text:
# A Julia Implementation of a streaming MedianFilter using Heaps with the same syntax as Matlabs medfilt1
module MedianFilter
# Dependencies
using DataStructures
import Base.isless, Base.isequal

export medfilt1

"""
    medfilt1(x::Array{<:Real}; n::Int = 3, padding::String = "zeropad", dim::Int = -1)

Apply a median filter to a signal vector x, similar to Matlabs medfilt1. Using Heap based calculation of the median to increase performance for larger windows n.

# Args:

* 'x::Array{<:Real}': Array containing real values

#  Keywords:

* 'n::Int': Window length. The Median at point i is defined as median(x[i-n+1:i])
* 'padding::String': Specifies how to deal with Endpoints. The modes 'zeropad' and 'truncate' are available, with the first as default.
* 'dim::Int': Specifies the dimension to be filtered along. As default the first non singleton dimension is chosen.

# Return:

* 'Array{Float64,N}': Always type Float64 with the same length as the input x

# Examples

```jldoctest
julia> medfilt1(collect(1:10))
10-element Array{Float64,1}:
 1.0
 2.0
 3.0
 4.0
 5.0
 6.0
 7.0
 8.0
 9.0
 9.0
```
"""
function medfilt1(x::Array{<:Real}; n::Int = 3, padding::String = "zeropad", dim::Int = -1)
    #Check Input
    @assert n > 0 "The window n needs to be greater 0. Is: n = $n"
    @assert padding in ["truncate", "zeropad"] "$padding is not a valid mode for padding. Use 'truncate' or 'zeropad'."
    @assert maximum(size(x)) > 1 "There is no non singleton dimension in the input."
    if dim != -1
        @assert dim > 0 "The dimension needs to be positive. Is: dim = $dim"
        @assert length(size(x)) >= dim "The input has only $(length(size(x))) dimensions but dim = $dim."
    else
        dim = find_first_nonsingleton(x)
    end

    if n == 1
        return x
    elseif padding == "zeropad"
        return mapslices(input -> medfilt1_worker_zeropad(input, n), x, dims = [dim])
    else
        return mapslices(input -> medfilt1_worker_truncate(input, n), x, dims = [dim])
    end

end

# Helper function to find the first non singleton dimension
function find_first_nonsingleton(x)
    dim = 1
    for i in size(x)
        if i > 1
            return dim
        end
        dim += 1
    end
    return dim
end

# The working horse, median Filter of an Vector with zero padding
function medfilt1_worker_zeropad(x::Vector{<:Real}, n::Int)
    if length(x) == 1
        return x
    end
    # adding zeros for padding
    x_pad = zeros(typeof(x[1]), (n -1))
    x_iter = Base.Iterators.Stateful([x_pad; x; x_pad])

    # initialization phase
    (result, left, right, delete_queue) = initialize_filter(x_iter)
    currentMedian = result[2]
    # loop through all elements
    while !isempty(x_iter)
        # delete Element if n elements are in the delete_queue
        if length(delete_queue) == n
            delete_element = dequeue!(delete_queue)
            if delete_element.in_left
                delete!(left, delete_element.handle)
            else
                delete!(right, delete_element.handle)
            end
        end

        # add new item
        # compare to current median
        currentElement = MedianElement(popfirst!(x_iter), false, -1)
        if currentElement.value >= currentMedian
            # add to left
            currentElement.handle = push!(right, currentElement)
        else
            # add to right
            currentElement.in_left = true
            currentElement.handle = push!(left, currentElement)
        end
        enqueue!(delete_queue, currentElement)

        # correct heaps if length differs more than one and get current median
        currentMedian = correct_heaps_return_median!(left, right)
        push!(result, currentMedian);
    end

    # cut result to the right length
    start = n + Int(ceil(n / 2)) - 1
    return result[start:start + length(x)-1]
end

# The working horse, median Filter of an Vector with smaller windows n at the endpoints
function medfilt1_worker_truncate(x::Vector{<:Real}, n::Int)
    if length(x) == 1
        return x
    end
    x_iter = Base.Iterators.Stateful(x)

    # initialization phase
    (result, left, right, delete_queue) = initialize_filter(x_iter)
    currentMedian = result[2]
    # loop through all elements an expand window (n_actual) to desired n
    n_actual = 2
    while length(delete_queue) > 1
        # delete Element if desired window size n was reached
        if n_actual >= n
            delete_element = dequeue!(delete_queue)
            if delete_element.in_left
                delete!(left, delete_element.handle)
            else
                delete!(right, delete_element.handle)
            end
        end

        # for each repetition increase n_actual
        n_actual += 1
        # add new item if there is something left in the iterator
        # compare to current median
        if !isempty(x_iter)
            currentElement = MedianElement(popfirst!(x_iter), false, -1)
            if currentElement.value >= currentMedian
                # add to left
                currentElement.handle = push!(right, currentElement)
            else
                # add to right
                currentElement.in_left = true
                currentElement.handle = push!(left, currentElement)
            end
            enqueue!(delete_queue, currentElement)
        end

        # correct heaps if length differs more than one and get current median
        currentMedian = correct_heaps_return_median!(left, right)
        push!(result, currentMedian);
    end

    # return desired excerpt
    start = Int(ceil((n) / 2))
    return result[start:start + length(x) - 1]
end

# Helper function that takes to Heaps, corrects them if length differs by more than 2 and returns the current median
function correct_heaps_return_median!(left, right)
    len_left = length(left)
    len_right = length(right)
    if abs(len_left - len_right) >= 2
        len_left > len_right ? transferTop!(left, right) : transferTop!(right, left)
        if abs(len_left - len_right) == 2
            currentMedian = 0.5 * top(left).value + 0.5 * top(right).value
            return currentMedian
        end
    end
    if len_left == len_right
        currentMedian = 0.5 * top(left).value + 0.5 * top(right).value
    elseif len_left > len_right
        currentMedian = top(left).value
    else
        currentMedian = top(right).value
    end
    return currentMedian
end

# Helper function that initializes all data structures and add first 2 elements
function initialize_filter(x_iter)
    # result vector
    result = Vector{Float64}()
    # >= current Median
    left = MutableBinaryMaxHeap{MedianElement}()
    # < current Median
    right = MutableBinaryMinHeap{MedianElement}()
    # A Queue that saves the Elements in the Heaps and pops the ones to delete next
    delete_queue = Queue{MedianElement}()

    # initialization phase
    # first value
    currentMedian = popfirst!(x_iter)
    push!(result, currentMedian);
    currentElement = MedianElement(currentMedian, true, -1)
    currentElement.handle  = push!(left, currentElement)
    enqueue!(delete_queue, currentElement)

    # second value
    currentElement = MedianElement(popfirst!(x_iter), false, -1)
    currentElement.handle = push!(right, currentElement)
    enqueue!(delete_queue, currentElement)
    if top(right).value < top(left).value
        swapTops!(left, right)
    end
    currentMedian = 0.5 * top(left).value + 0.5 * top(right).value
    push!(result, currentMedian);

    return (result, left, right, delete_queue)
end

# Data structure to store values and later identify them in a Heap to delete them
mutable struct MedianElement
    value
    in_left::Bool
    handle::Int
end

### Helper functions to get the Heap running on custom type MedianElement
function isless(x::MedianElement, y::MedianElement)
    return x.value < y.value
end
function isequal(x::MedianElement, y::MedianElement)
    return x.value == y.value
end
###

# Helper function to swap tops of the Heap
function swapTops!(left, right)
    temp_left = pop!(left)
    temp_left.in_left = false
    temp_right = pop!(right)
    temp_right.in_left = true
    temp_right.handle = push!(left, temp_right)
    temp_left.handle = push!(right, temp_left)
end

# Helper function to transfer the top to the smaller Heap
function transferTop!(larger, smaller)
    temp = pop!(larger)
    temp.in_left = !temp.in_left
    temp.handle = push!(smaller, temp)
end

end  # module MedianFilter


text:
module test_julia_modules1

using Test

module B
export c, d
f() = 2
end

@test [:B, :c, :d] == names(B)
@test test_julia_modules1 === parentmodule(B)
@test B === parentmodule(B.f)
@test (:test_julia_modules1, :B) == fullname(B)[end-1:end]
@test :B === nameof(B)
@test Base.isexported(B, :c)

@test Core === parentmodule(Int)
@test @isdefined B
@test !@isdefined Meh
@test isdefined(@__MODULE__, :B)

@test length(Base.loaded_modules) >= 33

end # module test_julia_modules1

module test_julia_modules2

using Test

module A
f(a::T) where T<:AbstractArray = 2
end # test_julia_modules2.A

module B
struct MyArray <: AbstractArray{Int,1}
end
Base.size(::MyArray) = 0
end # test_julia_modules2.B

using .A: f
using .B: MyArray

@test f(MyArray()) == 2

end # module test_julia_modules2

module test_julia_modules3

global a = 2

module A
global a = 3
end

using Test
@test a == 2

using .A
@test a == 2

end # module test_julia_modules3

module test_julia_modules4

using Test
m = Base.require(Base.PkgId(
                     Base.UUID((0x37e2e46d_f89d_539d,0xb4ee_838fcccc9c8e)), "LinearAlgebra"))
@test m isa Module

end # module test_julia_modules4

module test_julia_loaded_modules

pkg = Base.PkgId(Base.UUID("8dfed614-e22c-5e08-85e1-65c5234f0b40"), "Test")
loaded = Base.loaded_modules[pkg]
using .loaded

@test loaded === Test

end # module test_julia_loaded_modules

module test_julia_modules_macro

using Test

macro m()
    __module__
end

@test (@__MODULE__) === (@m) === test_julia_modules_macro

end # module test_julia_modules_macro

module test_julia_anonymous_module

using Test

m = Module()
@test m isa Module
@test nameof(m) === :anonymous
@test names(m, all=true, imported=true) == [:anonymous]

GC.@preserve m begin
    f = Base.eval(m, (x) -> 2x)
    @test f(3) == 6
end

end # module test_julia_anonymous_module

module test_julia_plus

using Test

module M
    + = 3
end

@test (+) !== (M.:+)

end # module test_julia_plus

module test_julia_baremodule

baremodule B
end

using Test

@test !isdefined(B, :Base)
@test isdefined(B, :Core)

baremodule A
+(x, y) = (x, y)
end

using .A: +

@test 1 + 2 == (1, 2)

end # module test_julia_baremodule

ext: jl


text:
using Test

function test_Group_interface(G::Group)
    @testset "Group interface" begin
        @testset "Iteration protocol" begin
            IS = Base.IteratorSize(typeof(G))
            if IS isa Base.HasLength || IS isa Base.HasShape
                @test isfinite(G) == true
                @test length(G) isa Int
                @test length(G) > 0

                @test eltype(G) <: GroupElement
                @test GroupsCore.elem_type(typeof(G)) == eltype(G)
                @test one(G) isa eltype(G)

                if GroupsCore.hasgens(G)
                    @test first(iterate(G)) isa eltype(G)
                    _, s = iterate(G)
                    @test first(iterate(G, s)) isa eltype(G)
                    @test isone(first(G))
                end
            else
                @test isfinite(G) == false
            end
        end

        @testset "Group generators" begin
            @test GroupsCore.hasgens(G) isa Bool

            if GroupsCore.hasgens(G)
                @test ngens(G) isa Int
                @test gens(G) isa AbstractVector{eltype(G)}
                @test length(gens(G)) == ngens(G)
                @test first(gens(G)) == gens(G, 1)
                @test last(gens(G)) == gens(G, ngens(G))
            else
                # TODO: throw something more specific
                @test_throws ErrorException gens(G)
                @test_throws ErrorException ngens(G)
            end
        end

        @testset "order, rand" begin
            if isfinite(G)
                @test order(Int16, G) isa Int16
                @test order(BigInt, G) isa BigInt
                @test order(G) >= 1
            else
                @test_throws GroupsCore.InfiniteOrder order(G)
            end

            @test rand(G) isa GroupElement
            @test rand(G, 2) isa AbstractVector{eltype(G)}
            g, h = rand(G, 2)
            @test parent(g) === parent(h) === G

            @test GroupsCore.rand_pseudo(G) isa eltype(G)
            @test GroupsCore.rand_pseudo(G, 2, 2) isa AbstractMatrix{eltype(G)}

            g, h = GroupsCore.rand_pseudo(G, 2)
            @test parent(g) === parent(h) === G
        end
    end
end

function test_GroupElement_interface(g::GEl, h::GEl) where {GEl<:GroupElement}

    @testset "GroupElement interface" begin

        @testset "Parent methods" begin
            @test parent(g) isa Group
            @test parent(g) === parent(h)
            G = parent(g)

            @test eltype(G) == typeof(g)
            @test GroupsCore.parent_type(typeof(g)) == typeof(G)
            @test GroupsCore.elem_type(GroupsCore.parent_type(g)) == typeof(g)

            @test one(g) isa eltype(G)

            @test one(G) == one(g) == one(h)

            if GroupsCore._is_deepcopiable(g)
                @test one(G) !== one(g)
            end

            @test isone(one(G))
        end

        @testset "Equality, deepcopy && hash" begin
            @test (g == h) isa Bool

            @test ==(g, h) isa Bool
            @test isequal(g, g)
            @test ==(h, h)

            if g != h
                @test !isequal(g, h)
            end

            @test deepcopy(g) isa typeof(g)
            @test deepcopy(g) == g
            if GroupsCore._is_deepcopiable(g)
                @test deepcopy(g) !== g
            end
            k = deepcopy(g)
            @test parent(k) === parent(g)
            @test hash(g) isa UInt
            @test hash(g) == hash(k)

            if isequal(g, h)
                @test hash(g) == hash(h)
            end
        end

        @testset "Group operations" begin
            old_g, old_h = deepcopy(g), deepcopy(h)

            # check that the default operations don't mutate their arguments
            @test inv(g) isa typeof(g)
            @test (g, h) == (old_g, old_h)

            @test g * h isa typeof(g)
            @test (g, h) == (old_g, old_h)

            @test g^2 == g * g
            @test (g, h) == (old_g, old_h)

            @test g^-3 == inv(g) * inv(g) * inv(g)
            @test (g, h) == (old_g, old_h)

            @test (g * h)^-1 == inv(h) * inv(g)
            @test (g, h) == (old_g, old_h)

            @test conj(g, h) == inv(h) * g * h
            @test (g, h) == (old_g, old_h)

            @test ^(g, h) == inv(h) * g * h
            @test (g, h) == (old_g, old_h)

            @test commutator(g, h) == g^-1 * h^-1 * g * h
            @test (g, h) == (old_g, old_h)

            @test commutator(g, h, g) == conj(inv(g), h) * conj(conj(g, h), g)
            @test (g, h) == (old_g, old_h)

            @test isone(g * inv(g)) && isone(inv(g) * g)
            @test (g, h) == (old_g, old_h)

            @test g / h == g * inv(h)
            @test (g, h) == (old_g, old_h)
        end

        @testset "Misc GroupElement methods" begin
            @test one(g) isa typeof(g)
            @test isone(g) isa Bool
            @test isone(one(g))

            @test isfiniteorder(g) isa Bool

            if isfiniteorder(g)
                @test order(Int16, g) isa Int16
                @test order(BigInt, g) isa BigInt
                @test order(g) >= 1
                if isfinite(parent(g))
                    @test iszero(order(parent(g)) % order(g))
                end
                if !isone(g) && !isone(g^2)
                    @test order(g) > 2
                end
                @test order(inv(g)) == order(g)
                @test order(one(g)) == 1
            else
                @test_throws GroupsCore.InfiniteOrder order(g)
            end

            @test similar(g) isa typeof(g)
        end

        one!, inv!, mul!, conj!, commutator!, div_left!, div_right! = (
            GroupsCore.one!,
            GroupsCore.inv!,
            GroupsCore.mul!,
            GroupsCore.conj!,
            GroupsCore.commutator!,
            GroupsCore.div_left!,
            GroupsCore.div_right!,
        )

        @testset "In-place operations" begin
            old_g, old_h = deepcopy(g), deepcopy(h)
            out = similar(g)

            @test isone(one!(g))
            g = deepcopy(old_g)

            @test inv!(out, g) == inv(old_g)
            @test g == old_g
            @test inv!(out, g) == inv(old_g)
            g = deepcopy(old_g)

            @testset "mul!" begin
                @test mul!(out, g, h) == old_g * old_h
                @test (g, h) == (old_g, old_h)

                @test mul!(out, g, h) == old_g * old_h
                @test (g, h) == (old_g, old_h)

                @test mul!(g, g, h) == old_g * old_h
                @test h == old_h
                g = deepcopy(old_g)

                @test mul!(h, g, h) == old_g * old_h
                @test g == old_g
                h = deepcopy(old_h)

                @test mul!(g, g, g) == old_g * old_g
                g = deepcopy(old_g)
            end

            @testset "conj!" begin
                res = old_h^-1 * old_g * old_h
                @test conj!(out, g, h) == res
                @test (g, h) == (old_g, old_h)

                @test conj!(g, g, h) == res
                @test h == old_h
                g = deepcopy(old_g)

                @test conj!(h, g, h) == res
                @test g == old_g
                h = deepcopy(old_h)

                @test conj!(g, g, g) == old_g
                g = deepcopy(old_g)
            end

            @testset "commutator!" begin
                res = old_g^-1 * old_h^-1 * old_g * old_h

                @test commutator!(out, g, h) == res
                @test (g, h) == (old_g, old_h)

                @test commutator!(out, g, h) == res
                @test (g, h) == (old_g, old_h)

                @test commutator!(g, g, h) == res
                @test h == old_h
                g = deepcopy(old_g)

                @test commutator!(h, g, h) == res
                @test g == old_g
                h = deepcopy(old_h)
            end

            @testset "div_[left|right]!" begin
                res = g * h^-1
                @test div_right!(out, g, h) == res
                @test (g, h) == (old_g, old_h)

                @test div_right!(g, g, h) == res
                @test h == old_h
                g = deepcopy(old_g)

                @test div_right!(h, g, h) == res
                @test g == old_g
                h = deepcopy(old_h)

                @test div_right!(g, g, g) == one(g)
                g = deepcopy(old_g)

                res = h^-1 * g
                @test div_left!(out, g, h) == res
                @test (g, h) == (old_g, old_h)

                @test div_left!(g, g, h) == res
                @test h == old_h
                g = deepcopy(old_g)

                @test div_left!(h, g, h) == res
                @test g == old_g
                h = deepcopy(old_h)

                @test div_left!(g, g, g) == one(g)
                g = deepcopy(old_g)
            end
        end
    end
end


text:
using GridapGeosciences
using Test

@testset "GridapGeosciences" begin
  @time @testset "CubedSphereDiscreteModelsTests" begin include("CubedSphereDiscreteModelsTests.jl") end
  @time @testset "DarcyCubedSphereTests" begin include("DarcyCubedSphereTests.jl") end
  @time @testset "WeakDivPerpTests" begin include("WeakDivPerpTests.jl") end
  @time @testset "LaplaceBeltramiCubedSphereTests" begin include("LaplaceBeltramiCubedSphereTests.jl") end
  @time @testset "WaveEquationCubedSphereTests" begin include("WaveEquationCubedSphereTests.jl") end
  #@time @testset "Williamsom2ThetaMethodFullNewtonTests" begin include("Williamsom2ThetaMethodFullNewtonTests.jl") end
  @time @testset "Williamson2ShallowWaterExplicitTests" begin include("Williamson2ShallowWaterExplicitTests.jl") end
  @time @testset "Williamson2ShallowWaterRosenbrockTests" begin include("Williamson2ShallowWaterRosenbrockTests.jl") end
  @time @testset "Williamson2ShallowWaterIMEXTests" begin include("Williamson2ShallowWaterIMEXTests.jl") end
end


text:
#check the following:
#    -piezo isn't moving too fast
#    -sufficient time interval between camera exposure starts
#        -also emit warning if frame rate exceeds full chip max (100fps)
#    -sufficient time interval between exposure stop and next exposure start (based on jitter time, see MINIMUM_EXPOSURE_SEPARATION variable
#    -sufficient time between laser state changes

function check_piezos(sigs::AbstractVector{TS}; window_sz = 100) where TS<:ImagineSignal
    ps = getpositioners(getoutputs(sigs))
    vs = check_piezo(ps[1]; window_sz = window_sz)
    for i = 2:length(ps)
        check_piezo(ps[i]; window_sz = window_sz, val_state = vs)
    end
    return vs
end

#check piezo speed
#window_sz is the size of the sliding window (in samples) over which to measure cumulative distance moved
#If the piezo is commanded to move faster than the speed set by PIEZO_MAX_SPEED within the window, then the test fails
#current method used by Imagine(?):  integer value can never change by more than +/-1, and cumulative change in a 10-sample window stays below some max speed
#by contrast here we are using a simple cumulative change over 100 samples
#with 1M sample rate this equates to a max change of 8 raw bits per 100 samps
#TODO: make this more robust to oscillations that occur faster than a 100-sample period.  Currently we're just hoping that those don't happen.  The strict way to fix this
#would be to check all window sizes less than 100 (or at least a couple more) but that is expensive
function check_piezo(pos::ImagineSignal; window_sz = 100, val_state::ValidationState = ValidationState()) #val_state is for bootstrapping with already-checked sequences
    max_sp = PIEZO_MAX_SPEED[rig_name(pos)]
    window_dur = window_sz / samprate(pos)
    max_dist = max_sp * window_dur #length units
    max_dist_raw = ceil(Int, (max_dist/width(interval_world(pos))) * width(interval_raw(pos)))
    val_func = (samps, win_sz) -> check_max_speed(samps, max_dist_raw, win_sz)
    if length(pos) < window_sz #should this throw an error?
        @warn "Insufficient samples to check this signal with a window size setting of $window_sz"
    else
        try
            window_validate!(val_state, val_func, window_sz, pos)
        catch
            error("The command $(name(pos)) exceeds the maximum speed of $max_sp permitted for this piezo / rig")
        end
    end
    return val_state
end

function check_max_speed(raw_samps::Vector, raw_change::Int, in_n_samps::Int)
    @assert length(raw_samps) >= in_n_samps
    niter = length(raw_samps) - in_n_samps + 1
    first_window = raw_samps[1:in_n_samps]
    last_samp = first_window[end]
    diff_vec = diff(first_window)
    #abs_vec will be a circular buffer for efficiency
    abs_vec = abs.(diff_vec)
    circ_i = 1
    tot = sum(abs_vec)
    for i = 2:niter
        if tot > raw_change
            return false
        end
        next_samp = raw_samps[i+in_n_samps-1]
        #update diff_vec, abs_vec, tot, and circ_i
        tot -= abs_vec[circ_i]
        next_abs = abs(next_samp - last_samp)
        abs_vec[circ_i] = next_abs
        tot += next_abs
        circ_i = circ_i % (in_n_samps-1) + 1
        last_samp = next_samp
    end
    #check last sample
    if tot > raw_change
        return false
    end
    return true
end

function check_pulse_duration(start_is::Vector{Int}, stop_is::Vector{Int}, nsamps_on_tol::Int, nsamps_off_tol::Int)
    bad_i = 0
    if nsamps_on_tol > 0
        for i = 1:length(start_is)
            if stop_is[i] - start_is[i] + 1 < nsamps_on_tol
                bad_i = i
                return bad_i
            end
        end
    end
    return bad_i
end

function check_interpulse_duration(start_is::Vector{Int}, stop_is::Vector{Int}, nsamps_on_tol::Int, nsamps_off_tol::Int)
    bad_i = 0
    if nsamps_off_tol > 0
        for i = 1:(length(start_is) - 1)
            if start_is[i+1] - stop_is[i] - 1 < nsamps_off_tol
                bad_i = i
                return bad_i
            end
        end
    end
    return bad_i
end

#checks that the start-to-start time is within tolerance
#also keeps track of the shortest start-to-start time
function check_pulse_interval(start_is::Vector{Int}, nsamps_interval_tol::Int)
    shortest_width = typemax(Int)
    bad_i = 0
    for i = 1:(length(start_is)-1)
        cur_width = start_is[i+1] - start_is[i]
        shortest_width = min(cur_width, shortest_width)
        if cur_width < nsamps_interval_tol
            bad_i = i
            return bad_i
        end
    end
    return bad_i #shortest_width
end

function check_pulses(sig::ImagineSignal, on_time::HasTimeUnits, off_time::HasTimeUnits, pulse_rate::HasInverseTimeUnits)
    rig = rig_name(sig)
    strts = find_pulse_starts(sig)
    stps = find_pulse_stops(sig)
    min_on_samps = ceil(Int, on_time * uconvert(inv(unit(on_time)), samprate(sig)))
    min_off_samps = ceil(Int, off_time * uconvert(inv(unit(off_time)), samprate(sig)))
    rslt = check_pulse_duration(strts, stps, min_on_samps, min_off_samps)
    if rslt > 0
        error("Pulse #$rslt of signal $(name(sig)) has insufficient duration")
    end
    rslt = check_interpulse_duration(strts, stps, min_on_samps, min_off_samps)
    if rslt > 0
        error("The interval between pulse #$rslt end and pulse #$(rslt+1) start of signal $(name(sig)) is too small")
    end
    nsamps_start_to_start = ceil(Int, (1/pulse_rate) * uconvert(unit(pulse_rate), samprate(sig)))
    if !isinf(pulse_rate)
        rslt = check_pulse_interval(strts, nsamps_start_to_start)
        if rslt > 0
            error("The interval between pulse #$rslt start and pulse #$(rslt+1) start of signal $(name(sig)) is too small")
        end
    end
    return true
end

check_cameras(sigs::AbstractVector{TS}; vertical_lines=nothing) where TS<:ImagineSignal = map(getcameras(getoutputs(sigs))) do cam
    default_chip_sz = chip_size(rig_name(cam))
    chip_sz = vertical_lines === nothing ? default_chip_sz : (default_chip_sz[1], vertical_lines)
    return check_camera(cam; chip_sz=chip_sz)
end

#check framerate, interpulse, intrapulse
#   the 0->1->0 interval is greater than or equal to CAMERA_ON_TIME
#   the 1->0->1 interval is greater than or equal to CAMERA_OFF_TIME
#   the 0->1->0->1 interval is greater than or equal to 1 / the max framerate
#   max_framerate = max_framerate(rig, chip_size(rig)...)
function check_camera(cam::ImagineSignal; chip_sz = chip_size(rig_name(cam)))
    if isempty(cam)
        @warn "Signal $(name(cam)) is empty.  Skipping validation."
        return true
    end
    rig = rig_name(cam)
    min_on_dur = CAMERA_ON_TIME[rig]
    min_off_dur = CAMERA_OFF_TIME[rig]
    max_fr = max_framerate(rig, chip_sz...)
    check_pulses(cam, min_on_dur, min_off_dur, max_fr)
end

check_lasers(sigs::AbstractVector{TS}) where TS<:ImagineSignal = map(check_laser, getlasers(getoutputs(sigs)))
#check laser pulse duration and frequency
#Imagine currently checks every transition points(0->1) and calculate frequency between two consecutive transition points
#We check that:
#   the 0->1->0 interval is greater than or equal to LASER_ON_TIME
#   the 1->0->1 interval is greater than or equal to LASER_OFF_TIME
#   the 0->1->0->1 interval is unconstrained
function check_laser(las::ImagineSignal)
    if isempty(las)
        @warn "Signal $(name(las)) is empty.  Skipping validation."
        return true
    end
    rig = rig_name(las)
    min_on_dur = LASER_ON_TIME[rig]
    min_on_samps = min_on_dur * samprate(las)
    min_off_dur = LASER_OFF_TIME[rig]
    check_pulses(las, min_on_dur, min_off_dur, Inf*inv(Unitful.s))
end

function validate_singles(sigs::AbstractVector{TS}; vertical_lines=nothing) where TS<:ImagineSignal
    check_piezos(sigs)
    check_cameras(sigs; vertical_lines=vertical_lines)
    check_lasers(sigs)
end

function validate_all(sigs::AbstractVector{TS}; check_is_sufficient = true, vertical_lines=nothing) where TS<:ImagineSignal
    validate_group(sigs; check_is_sufficient = check_is_sufficient)
    validate_singles(sigs; vertical_lines=vertical_lines)
end

max_forks_repo_name: HolyLab/ImagineInterface


text:
# Controller
function controller!(mechanism, k;
    U=0.5,
    timestep=0.01)
    for joint in mechanism.joints
        nu = input_dimension(joint)
        u = (nu <= 5 && k ∈ (1:100)) * U * timestep * sones(nu)
        set_input!(joint, u)
    end
    return
end
nocontrol!(mechanism, k) = controller!(mechanism, k, U=0.0)

# Parameters
ϵ0 = 1e-12
timestep0 = 0.01
gravity0= 0.0
joint_types = [
    :Fixed,
    :Prismatic,
    :Planar,
    :FixedOrientation,
    :Revolute,
    :Cylindrical,
    :PlanarAxis,
    :FreeRevolute,
    :Orbital,
    :PrismaticOrbital,
    :PlanarOrbital,
    :FreeOrbital,
    :Spherical,
    :CylindricalFree,
    :PlanarFree
    ]

################################################################################
# Box
################################################################################
# single body
# initial linear and angular velocity
# no gravity
# no spring and damper
# no control
################################################################################
@testset "Box" begin
    mech = get_mechanism(:block,
        timestep=timestep0,
        gravity=gravity0,
        contact=false)

    v0 = [1,2,3.0]
    ω0 = [10,10,10.0]
    initialize!(mech, :block,
        v=v0, ω=ω0)

    storage = simulate!(mech, 5.0, nocontrol!,
        record=true,
        verbose=false,
        opts=SolverOptions(rtol=ϵ0, btol=ϵ0))
    # visualize(mech, storage, vis = vis)

    m0 = momentum(mech, storage)[5:end]
    mlin0 = [Vector(m - m0[1])[1:3] for m in m0]
    mang0= [Vector(m - m0[1])[4:6] for m in m0]
    # plot([(i-1) * timestep0 for i in 1:length(m0)], hcat(mlin0...)')
    # plot([(i-1) * timestep0 for i in 1:length(m0)], hcat(mang0...)')
    @test all(norm.(mlin0, Inf) .< 1.0e-7)
    @test all(norm.(mang0, Inf) .< 1.0e-7)
end

################################################################################
# SINGLE PENDULUM
################################################################################
# single body
# initial angular velocity
# no gravity
# no spring and damper
# no control
################################################################################
@testset "Pendulum" begin
    mech = get_mechanism(:pendulum,
        timestep=timestep0,
        gravity=gravity0)

    ϕ0 = 0.7
    ω0 = 5.0
    initialize!(mech, :pendulum,
        angle=ϕ0,
        angular_velocity=ω0)

    storage = simulate!(mech, 5.0, nocontrol!,
        record=true,
        verbose=false,
        opts=SolverOptions(rtol=ϵ0, btol=ϵ0))

    # visualize(mech, storage, vis = vis)

    m0 = momentum(mech, storage)[10:end]
    mlin0 = [Vector(m - m0[1])[1:3] for m in m0]
    mang0= [Vector(m - m0[1])[4:6] for m in m0]
    # plot([(i-1) * timestep0 for i in 1:length(m0)], hcat(mlin0...)')
    # plot([(i-1) * timestep0 for i in 1:length(m0)], hcat(mang0...)')
    # @test all(norm.(mlin0, Inf) .< 1e-11)
    @test all(norm.(mang0, Inf) .< 1.0e-7)
end

################################################################################
#  HUMANOID
################################################################################
# multiple bodies
# initial linear and angular velocity
# no gravity
# with spring and damper
# with control
################################################################################
@testset "Humanoid" begin
    spring0 = 1.0
    damper0 = 1.0
    mech = get_mechanism(:humanoid,
        timestep=timestep0,
        gravity=gravity0,
        spring=spring0,
        damper=damper0,
        contact_feet=false,
        contact_body=false)

    initialize!(mech, :humanoid)
    bodies = mech.bodies
    set_maximal_velocities!.(bodies, ω=1e-0rand(3))

    storage = simulate!(mech, 10.0, controller!,
        record=true,
        verbose=false,
        opts=SolverOptions(rtol=ϵ0, btol=ϵ0))

    # visualize(mech, storage, vis = vis)

    m0 = momentum(mech, storage)[1:end]
    mlin0 = [Vector(m - m0[1])[1:3] for m in m0]
    mang0= [Vector(m - m0[1])[4:6] for m in m0]
    # plot([(i-1) * timestep0 for i in 1:length(m0)], hcat(mlin0...)')
    # plot([(i-1) * timestep0 for i in 1:length(m0)], hcat(mang0...)')
    @test all(norm.(mlin0, Inf) .< 1.0e-8)
    @test all(norm.(mang0, Inf) .< 1.0e-8)
end

################################################################################
#  ATLAS
################################################################################
# multiple bodies
# initial linear and angular velocity
# no gravity
# with spring and damper
# with control
################################################################################
@testset "Atlas" begin
    spring0 = 10.0
    damper0 = 1.0
    mech = get_mechanism(:atlas,
        timestep=timestep0,
        gravity=gravity0,
        spring=spring0,
        damper=damper0,
        contact_feet=false,
        contact_body=false)

    initialize!(mech, :atlas)
    storage = simulate!(mech, 5.0, controller!,
        record=true,
        verbose=false,
        opts=SolverOptions(rtol=ϵ0, btol=ϵ0))

    # visualize(mech, storage, vis = vis)

    m0 = momentum(mech, storage)[1:end]
    mlin0 = [Vector(m - m0[1])[1:3] for m in m0]
    mang0= [Vector(m - m0[1])[4:6] for m in m0]
    # plot([(i-1) * timestep0 for i in 1:length(m0)], hcat(mlin0...)')
    # plot([(i-1) * timestep0 for i in 1:length(m0)], hcat(mang0...)')
    @test all(norm.(mlin0, Inf) .< 1.0e-8)
    @test all(norm.(mang0, Inf) .< 1.0e-8)
end

################################################################################
#  QUADRUPED
################################################################################
# multiple bodies
# initial linear and angular velocity
# no gravity
# with spring and damper
# with control
################################################################################
@testset "Quadruped" begin
    spring0 = 0.3
    damper0 = 0.1
    mech = get_mechanism(:quadruped,
        timestep=timestep0,
        gravity=gravity0,
        spring=spring0,
        damper=damper0,
        contact_feet=false,
        contact_body=false)

    initialize!(mech, :quadruped)
    storage = simulate!(mech, 5.0, controller!,
        record=true,
        verbose=false,
        opts=SolverOptions(rtol=ϵ0, btol=ϵ0))

    # visualize(mech, storage, vis = vis)

    m0 = momentum(mech, storage)[1:end]
    mlin0 = [Vector(m - m0[1])[1:3] for m in m0]
    mang0= [Vector(m - m0[1])[4:6] for m in m0]
    # plot([(i-1) * timestep0 for i in 1:length(m0)], hcat(mlin0...)')
    # plot([(i-1) * timestep0 for i in 1:length(m0)], hcat(mang0...)')
    @test all(norm.(mlin0, Inf) .< 1.0e-8)
    @test all(norm.(mang0, Inf) .< 1.0e-8)
end

################################################################################
# 5-lINK SNAKE
################################################################################
# multiple bodies
# initial linear and angular velocity
# no gravity
# with spring and damper
# with control
################################################################################
@testset "Snake" begin
    Nb0 = 5
    spring0 = 0.0 * 4e0
    damper0 = 0.0 * 2e+1

    mech = get_mechanism(:snake,
        timestep=timestep0,
        gravity=gravity0,
        num_bodies=Nb0,
        spring=spring0,
        damper=damper0,
        joint_type=:Revolute,
        contact=false,
        radius=0.05);

    v0 = 100.0 * [1, 2, 3] * timestep0
    ω0 = 100.0 * [1, 2, 3.0] * timestep0
    q10 = Dojo.RotX(0.5*π)

    initialize!(mech, :snake,
        base_orientation=q10,
        base_linear_velocity=v0,
        base_angular_velocity=ω0)
    storage = simulate!(mech, 1.50,
        record=true,
        verbose=false,
        opts=SolverOptions(rtol=ϵ0, btol=ϵ0))

    # visualize(mech, storage, vis = vis)

    m0 = momentum(mech, storage)[5:end]
    mlin0 = [Vector(m - m0[1])[1:3] for m in m0]
    mang0= [Vector(m - m0[1])[4:6] for m in m0]
    # plot([(i-1) * timestep0 for i in 1:length(m0)], hcat(mlin0...)')
    # plot([(i-1) * timestep0 for i in 1:length(m0)], hcat(mang0...)')
    @test all(norm.(mlin0, Inf) .< 1.0e-8)
    @test all(norm.(mang0, Inf) .< 1.0e-8)

    for joint_type in joint_types
        mech = get_mechanism(:snake,
            timestep=timestep0,
            gravity=gravity0,
            num_bodies=Nb0,
            spring=spring0,
            damper=damper0,
            joint_type=joint_type,
            contact=false,
            radius=0.05)

        v0 = 10.0 * [1, 2, 3] * timestep0
        ω0 = 10.0 * [1, 2, 3.0] * timestep0
        q10 = Dojo.RotX(0.5 * π)

        initialize!(mech, :snake,
            base_orientation=q10,
            base_linear_velocity=v0,
            base_angular_velocity=ω0)
        storage = simulate!(mech, 1.50, controller!,
            record=true,
            verbose=false,
            opts=SolverOptions(rtol=ϵ0, btol=ϵ0))

        # visualize(mech, storage, vis = vis)

        m0 = momentum(mech, storage)[5:end]
        # @show mlin0 = [Vector(m - m0[1])[1:3] for m in m0]
        # @show mang0= [Vector(m - m0[1])[4:6] for m in m0]
        # plot([(i-1) * timestep0 for i in 1:length(m0)], hcat(mlin0...)')
        # plt = plot()
        # plot!([(i-1) * timestep0 for i in 1:length(m0)], hcat(mang0...)')
        # display(plt)
        @test all(norm.(mlin0, Inf) .< 1.0e-8)
        @test all(norm.(mang0, Inf) .< 1.0e-8)
    end
end

################################################################################
# 5-lINK TWISTER
################################################################################
# multiple bodies
# initial linear and angular velocity
# no gravity
# with spring and damper
################################################################################
@testset "Twister" begin
    Nb0 = 5
    spring0 = 1.0 * 4e0
    damper0 = 1.0 * 2e+1

    mech = get_mechanism(:twister,
        timestep=timestep0,
        gravity=gravity0,
        num_bodies=Nb0,
        spring=spring0,
        damper=damper0,
        joint_type=:FixedOrientation,
        contact=false,
        radius=0.05);

    v0 = 100.0 * [1, 2, 3] * timestep0
    ω0 = 100.0 * [1, 2, 3.0] * timestep0
    q10 = Dojo.RotX(0.5*π)

    initialize!(mech, :twister,
        base_orientation=q10,
        base_linear_velocity=v0,
        base_angular_velocity=ω0)
    storage = simulate!(mech, 2.5,
        record=true,
        verbose=false,
        opts=SolverOptions(rtol=ϵ0, btol=ϵ0))

    # visualize(mech, storage, vis = vis)

    m0 = momentum(mech, storage)[5:end]
    mlin0 = [Vector(m - m0[1])[1:3] for m in m0]
    mang0= [Vector(m - m0[1])[4:6] for m in m0]
    # plot([(i-1) * timestep0 for i in 1:length(m0)], hcat(mlin0...)')
    # plot([(i-1) * timestep0 for i in 1:length(m0)], hcat(mang0...)')
    @test all(norm.(mlin0, Inf) .< 1.0e-8)
    @test all(norm.(mang0, Inf) .< 1.0e-8)

    for joint_type in joint_types
        mech = get_mechanism(:twister,
            timestep=timestep0,
            gravity=gravity0,
            num_bodies=Nb0,
            spring=spring0,
            damper=damper0,
            joint_type=joint_type,
            contact=false,
            radius=0.05)

        v0 = 10.0 * [1, 2, 3] * timestep0
        ω0 = 10.0 * [1, 2, 3.0] * timestep0
        q10 = Dojo.RotX(0.5*π)
        initialize!(mech, :twister,
            base_orientation=q10,
            base_linear_velocity=v0,
            base_angular_velocity=ω0)
        storage = simulate!(mech, 1.50,
            record=true,
            verbose=false,
            opts=SolverOptions(rtol=ϵ0, btol=ϵ0))
        # visualize(mech, storage, vis = vis)

        m0 = momentum(mech, storage)[5:end]
        mlin0 = [Vector(m - m0[1])[1:3] for m in m0]
        mang0= [Vector(m - m0[1])[4:6] for m in m0]

        @test all(norm.(mlin0, Inf) .< 1.0e-8)
        @test all(norm.(mang0, Inf) .< 1.0e-8)
    end
end


text:
function permutest(disp ::Disp,n_perm = 10000)
    """
    This function permutes the residuals returned by 'dispersion' to provide global and pairwise
    P-values. Currently, no corrections are made for multiple testing. MultipleTesting.jl provides
    this functionality, if desired

    The function accepts:
        disp: The named tuple returned by 'dispersion'
        n_perm: the number of permutations

    This function returns a named tuple containing:
        P = Global P-value
        pairwise_P = pairwise P-value
        F = Global F-Statistic
        pairwise_F = pairwise F-statistics

    https://github.com/juliangehring/MultipleTesting.jl
    """
    # read in values from `disp`
    F = disp.F
    R = disp.residuals
    F_pairs = disp.pairwise_F
    group = disp.group
    levels = disp.levels
    level_names = string.(levels)
    inds = [group .== level for level in levels]
    n = length(levels)

    # pre-calculate values for ANOVA function
    r = vcat(R...)
    X = vcat(R...)
    N = length(X)
    nj = Tuple(length.(R))
    X̅ = mean(X)
    k = length(R)
    # generate empty containers for F values
    perm = Vector{Float64}(undef,n_perm)
    #rs = copy(R)#Vector{Vector}(undef,n)
    # run permutation
    @inbounds for p in 1:n_perm
        shuffle!(r)
        rs = [view(r,inds[i]) for i in 1:n]

        perm[p] =f(rs,N,nj,X̅,k)
    end
    ppairs =zeros(k,k)

    Threads.@threads for i in 1:k-1
        for j in i+1:k
            r_1 = r[group .== levels[i]]
            r_2 = r[group .== levels[j]]
            pstat =permutest(r_1,r_2, F_pairs[j,i],n_perm)
            ppairs[j,i] = pstat
        end
    end
    # calculate P-values and return P and F values
    P = sum(perm .> F)/n_perm
    P_pairs = NamedArray(ppairs,( level_names,level_names), ("group","group"))
    F_pairs = NamedArray(LowerTriangular(F_pairs),( level_names,level_names ), ("group","group"))

    return (P =P,pairwise_P = P_pairs,F = F, pairwise_F =F_pairs)
end

function permutest(r_1,r_2 ,F,n_perm = 1000)
    r = vcat(r_1,r_2)
    nj = (length(r_1),length(r_2))
    n = 2
    N = sum(nj)
    # pre-calculate values for ANOVA function
    X̅ = mean(r)
    k = 2
    inds = (1:nj[1],1:nj[2])
    # generate empty containers for F values
    perm = Vector{Float64}(undef,n_perm)
    # run permutation
    @inbounds for p in 1:n_perm
        shuffle!(r)
        rs = (view(r,inds[1]) ,view(r,inds[2]) )
        perm[p] =f(rs,N,nj,X̅,k)
    end
    P = sum(perm .> F)/n_perm
    return P
end


text:
using Makie
using SimplexQuad

vertices = [
    0.0 0.0
    1.0 0.0
    0.0 1.0
]
# Need at least 4 vertices for plotting?
vertices1 = [vertices; vertices[1:1, :]]

for npoints in [1, 2, 3, 4, 10, 100]
    X, W = simplexquad(npoints, vertices)

    scene = Scene()
    lines!(scene, vertices1; color = :black, linewidth = 2)
    scatter!(scene, X[:, 1], X[:, 2]; color = :red, markersize = 0.3 * sqrt.(W))
    text!(scene, "N=$npoints"; textsize = 0.15, position = (0.5, 0.9))
    scale!(scene, 1, 1)

    Makie.save("figures/gauß-points-$npoints.png", scene;
               resolution = (300, 300))
end


text:
using Formatting

function chop_suffix_from_path(path, suffix)
    if endswith(path, suffix)
        return path[1:(end-length(suffix))]
    end
    return path
end

function test_sleep(fd, ntests)
    for _ in 1:ntests
        elapsed_time = @elapsed sleep(100)
        println("Sleep(100): $elapsed_time")
        write(fd, string("A\t", fmt(".6f", elapsed_time), "\t1\n"))
    end
    write(fd, "\n")
end

function test_for_sleep(fd, ntests)
    niter = 100
    for _ in 1:ntests
        elapsed_time = @elapsed for _ in 1:niter
                        sleep(1)
        end
        println("For 100 times sleep(1): $elapsed_time")
        write(fd, string("B\t", fmt(".6f", elapsed_time), "\t$niter\n"))
    end
    write(fd, "\n")
end

function test_add(fd, ntests, nincr)
    for _ in 1:ntests
        initial_value, val = 0, 0
        elapsed_time = @elapsed for _ in 1:nincr
                        val += 1.0
        end
        println("For $nincr add 1 to $initial_value: $elapsed_time")
        write(fd, string("C\t", fmt(".6f", elapsed_time), "\t$nincr\n"))
    end
    write(fd, "\n")
end

function test_add_if(fd, ntests, nincr)
    for _ in 1:ntests
        initial_value, val = 0, 0
        elapsed_time = @elapsed for i in 1:nincr
                        if i % 2 == 0
                            val += 1.0
                        end
        end
        println("For $nincr add 1 to $initial_value: $elapsed_time")
        write(fd, string("D\t", fmt(".6f", elapsed_time), "\t$nincr\n"))
    end
    write(fd, "\n")
end

function test_subtract(fd, ntests, nincr)
    for _ in 1:ntests
        initial_value, val = 1000000, 1000000
        elapsed_time = @elapsed for _ in 1:nincr
                        val -= 1.0
        end
        println("For $nincr subtract 1 from $initial_value: $elapsed_time")
        write(fd, string("E\t", fmt(".6f", elapsed_time), "\t$nincr\n"))
    end
    write(fd, "\n")
end

function test_subtract_if(fd, ntests, nincr)
    for _ in 1:ntests
        initial_value, val = 1000000, 1000000
        elapsed_time = @elapsed for i in 1:nincr
                        if i % 2 == 0
                            val -= 1.0
                        end
        end
        println("For $nincr subtract 1 from $initial_value: $elapsed_time")
        write(fd, string("F\t", fmt(".6f", elapsed_time), "\t$nincr\n"))
    end
    write(fd, "\n")
end

function test_multiply(fd, ntests, nincr)
    for _ in 1:ntests
        initial_value, val = 1, 1
        elapsed_time = @elapsed for _ in 1:nincr
                        val *= 1.01
        end
        println("For $nincr multiply $initial_value by 1.01: $elapsed_time")
        write(fd, string("G\t", fmt(".6f", elapsed_time), "\t$nincr\n"))
    end
    write(fd, "\n")
end

function test_multiply_if(fd, ntests, nincr)
    for _ in 1:ntests
        initial_value, val = 1, 1
        elapsed_time = @elapsed for i in 1:nincr
                        if i % 2 == 0
                            val *= 1.01
                        end
        end
        println("For $nincr multiply $initial_value by 1.01: $elapsed_time")
        write(fd, string("H\t", fmt(".6f", elapsed_time), "\t$nincr\n"))
    end
    write(fd, "\n")
end

function test_divide(fd, ntests, nincr)
    for _ in 1:ntests
        initial_value, val = 1000000, 1000000
        elapsed_time = @elapsed for _ in 1:nincr
                        val /= 1.01
        end
        println("For $nincr divide $initial_value by 1.01: $elapsed_time")
        write(fd, string("I\t", fmt(".6f", elapsed_time), "\t$nincr\n"))
    end
    write(fd, "\n")
end

function test_divide_if(fd, ntests, nincr)
    for _ in 1:ntests
        initial_value, val = 1000000, 1000000
        elapsed_time = @elapsed for i in 1:nincr
                        if i % 2 == 0
                            val /= 1.01
                        end
        end
        println("For $nincr divide $initial_value by 1.01: $elapsed_time")
        write(fd, string("J\t", fmt(".6f", elapsed_time), "\t$nincr\n"))
    end
    write(fd, "\n")
end

max_stars_repo_name: mjaglarz/BA-thesis

max_issues_repo_name: mjaglarz/BA-thesis

max_forks_repo_name: mjaglarz/BA-thesis


text:
#from PC-SAFT paper https://doi.org/10.1021/ie0003887
function fugacity_coeff_impl(mt::MultiVT,model::HelmholtzModel,v,t,x)
    Z = compressibility_factor_impl(mt,model,v,t,x)
    ar = z -> αR_impl(mt,model,one(v)/v,t,z)
    dadx = ForwardDiff.gradient(ar,x)
    sumdadx = dot(x,dadx)
    lnZ = log(Z)
    c0 = ar(x) + (Z-1) - lnZ - sumdadx
    return map(xi->xi+c0,dadx)
end

function flash_impl(mt::MultiPT,model,p,t,z0;threaded=true,moles=one(eltype(z0)))
    p = float(p)
    t = float(t)
    _0 = zero(p*t*first(z0))
    _1 = one(_0)
    logϕ0 = similar(z0)
    if threaded
        res_z0 = Threads.@spawn update_logϕ!!(model,logϕ0,$p,$t,z0,nothing,:undefined)
    end
    kmodel = WilsonK(model)
    k0 = kvalues_impl(mt,kmodel,p,t,z0)
    k = copy(k0)
    #@show k
    g_0 = flash_eval(k,z0,_0)
    g_1 = flash_eval(k,z0,_1)
    if g_0 <= 0  #bubble point assumption
        #println("bubble point")
        β = _0
        β0 = _0
        xil = copy(z0)
        xiv = flash_vapor(k,z0,_0)
    elseif g_1 >= 0 #dew point assumption
        #println("dew point")
        β = _1
        β0 = _1
        xil = flash_liquid(k,z0,_1)
        xiv = copy(z0)
    else #two phase assumption
        #println("flash point")
        β = flash_vfrac(RR(),k,z0)
        xil = flash_liquid(k,z0,β)
        xiv = flash_vapor(k,z0,β)
    end
    logϕl = similar(xil)
    logϕv = similar(xiv)
    xil = normalizefrac!!(xil)
    xiv = normalizefrac!!(xiv)

    vl,logϕl,vv,logϕv = update_logϕ_vl!!(model,p,t,z0,xil,xiv,logϕl,logϕv,nothing,nothing,threaded)
    lnk = logϕl - logϕv
    @! k .= exp.(lnk)

    stable_phase = false
    for i = 1:5
        res = refine_k!!(model,p,t,z0,xil,xiv,vl,vv,k,logϕl,logϕv,lnk,threaded)
        xil,xiv,vl,vv,k,logϕl,logϕv,lnk,β = res

        if !(_0 <= β <= _1)
            stable_phase = true
            break
        end
    end
    #if the process was threaded, the value is fetched here
    if !threaded
        vz0,logϕ0 = update_logϕ!!(model,logϕ0,p,t,z0,nothing,:undefined)
    else
        vz0,logϕ0 = fetch(res_z0)
    end

    if stable_phase == false
        tpdl = flash_tpd(xil, logϕl,z0,logϕ0)
        tpdv = flash_tpd(xiv, logϕv,z0,logϕ0)
        ΔG = (1 - β) * tpdl + β * tpdv
        if all(>(0),(ΔG,tpdv,tpdl))
            stable_phase = true # more detailed stability phase analisis
        end
    end

    if stable_phase == true
        mul_v = copy(xiv)
        mul_l = copy(xiv)
        initial_data = (model,p,t,z0)
        xil,xiv,vl,vv = trial_compositions!!(initial_data,xil,xiv,vl,vv,logϕl,logϕv,logϕ0,mul_v,mul_l,true;k0=k0)
        count = 1
        max_count = 10
        while count < max_count
            xil,xiv,vl,vv,tpdl,tpdv = trial_compositions!!(initial_data,xil,xiv,vl,vv,logϕl,logϕv,logϕ0,mul_v,mul_l,true)
            converged,stable_phase = trial_convergence(z0,xil,xiv,mul_l,mul_v,tpdl,tpdv,stable_phase)
            if converged == true
                break
            else
                count += 1
            end
        end
        if stable_phase == false
            for i in 1:3
                xil,xiv,vl,vv,tpdl,tpdv = trial_compositions!!(initial_data,xil,xiv,vl,vv,logϕl,logϕv,logϕ0,mul_v,mul_l,true)
            end
        end
    end

    if stable_phase == false
        vle_res = VLEResults(
            model = model
            ,vl = vl #liquid mol volume
            ,vv = vv #gas mol volume
            ,tv = t #temperature
            ,pv = p #pressure
            ,tl = t #temperature
            ,pl = p #pressure
            ,z0 = z0 #initial composition
            ,lnϕl = logϕl #efective log fugacity coef - liquid
            ,lnϕv = logϕv #efective log fugacity coef - gas
            ,nl = copy(xil) #moles of liquid
            ,nv = copy(xiv) #moles of gas
            ,xil = xil #mol fraction of liquid
            ,xiv = xiv) #mol fraction of gas
            #@show vle_res.xiv
            #@show vle_res.xil
            #@show vle_res.nv
            #@show vle_res.nl
            fg! = generate_ptflash_objective(vle_res,threaded=threaded)
            nv0 = zeros(length(xiv))
            nv0 .+= xiv
            nv0 .*= β
            opts = Optim.Options(g_tol = 1e-12,
                iterations = 10,
                store_trace = true,
                show_trace = true)
            optim_obj = OnceDifferentiable(only_fg!(fg!), nv0)
         try
            optimize(optim_obj,nv0, Optim.BFGS(linesearch = Optim.LineSearches.MoreThuente()))
         catch err
            @show vle_res
         finally
            return vtx_states(vle_res,moles)
         end
    end

    return state(t=t,v=vz0,xn = z0)
end

function trial_compositions!!(initial_data,xil,xiv,vl,vv,logϕl,logϕv,logϕ0,mul_v,mul_l,threaded=true;k0=nothing)
    model,p,t,z0 = initial_data
    if k0 !== nothing

        @! xiv .= z0 .* k0
        @! xil .= z0 ./ k0

        xil = normalizefrac!!(xil)
        xiv = normalizefrac!!(xiv)
        vl = nothing
        vv = nothing
        tpdv = convert(typeof(p),Inf)
        tpdl = convert(typeof(p),Inf)
    else
        vl,logϕl,vv,logϕv = update_logϕ_vl!!(model,p,t,z0,xil,xiv,logϕl,logϕv,nothing,nothing,threaded)
        mul_v = exp.(logϕ0 .- logϕv)
        mul_l = exp.(logϕ0 .- logϕl)
        @! xiv += z0 .* mul_v - xiv
        @! xil += z0 .* mul_l - xil
        xil = normalizefrac!!(xil)
        xiv = normalizefrac!!(xiv)
        tpdl = flash_tpd(xil, logϕl,z0,logϕ0)
        tpdv = flash_tpd(xiv, logϕv,z0,logϕ0)
        @show tpdl
        @show tpdv

    end
    return xil,xiv,vl,vv,tpdl,tpdv
end

function trial_convergence(z0,xil,xiv,mul_l,mul_v,tpdl,tpdv,stable_phase)
    new_stable_phase = stable_phase
    _ones = ones(length(mul_l))

    #this checks for:
    # a: the composition aproaches z0
    # b: variation between iterations is very small

    if mul_v ≈ _ones
        converged = true
     elseif mul_l ≈ _ones
         converged = true
     elseif xiv ≈ z0
        converged = true
    elseif xil ≈ z0
        converged = true
    elseif tpdl < 0
        converged = true
        new_stable_phase = false
    elseif tpdv < 0
        converged = true
        new_stable_phase = false
    else
        converged = false
    end
    return converged,new_stable_phase
end

function refine_k!!(model,p,t,z0,xil,xiv,vl,vv,k,logϕl,logϕv,lnk,threaded)
    vl,logϕl,vv,logϕv = update_logϕ_vl!!(model,p,t,z0,xil,xiv,logϕl,logϕv,vl,vv,threaded)
    #@show xil,xiv,vl,vv
    @! lnk .= logϕl .- logϕv
    @! k .= exp.(lnk)
   k = remove_maximums!!(k)
    β = flash_vfrac(RR(),k,z0)
    #@show β
    xil = flash_liquid!!(xil,k,z0,β)
    xiv = flash_vapor!!(xiv,k,z0,β)

    xil = normalizefrac!!(xil)
    xiv = normalizefrac!!(xiv)
    return xil,xiv,vl,vv,k,logϕl,logϕv,lnk,β

end

#update both liquid and gas fugacity coefficients, considering threading
function update_logϕ_vl!!(model,p,t,z0,xil,xiv,logϕl,logϕv,vl=nothing,vv=nothing,threaded=true)
    if threaded
        res_l = Threads.@spawn update_logϕ!!(model,logϕl,$p,$t,xil,$vl,:liquid)
        res_v = Threads.@spawn update_logϕ!!(model,logϕv,$p,$t,xiv,$vv,:vapor)
         vl,logϕl = fetch(res_l)
         vv,logϕv = fetch(res_v)
     else
         vl,logϕl = update_logϕ!!(model,logϕl,p,t,xil,vl,:liquid)
         vv,logϕv = update_logϕ!!(model,logϕv,p,t,xiv,vv,:vapor)
     end
     #@show logϕl
     #@show logϕv
     #@! lnk .= logϕl .- logϕv
     #@! k .= exp.(lnk)
     #@show k
     return vl,logϕl,vv,logϕv

end
#with values of p,t,x; returns a new logϕ,v
function update_logϕ!!(model,logϕ,p,t,x,v=nothing,phase=:liquid)
    x = normalizefrac!!(x)
    _vtx = QuickStates.vtx()
    _ptx = QuickStates.ptx()
    _v =  v_zero(_ptx,model,p,t,x,v;phase=phase)
    @! logϕ = fugacity_coeff_impl(_vtx,model,_v,t,x)
    return _v,logϕ
end

function flash_tpd(xi, lnϕ,x0,lnϕ0)
    f(xiᵢ,lnϕᵢ,x0ᵢ,lnϕ0ᵢ) = xiᵢ*(log(xiᵢ) + lnϕᵢ-log(x0ᵢ)- lnϕ0ᵢ)
    return mapreduce(f,+,xi, lnϕ,x0,lnϕ0)
end

Base.@kwdef mutable struct VLEResults{MODEL,S,V}
    model::MODEL
    vl::S #liquid mol volume
    vv::S #gas mol volume
    tv::S #temperature, gas
    pv::S #pressure, gas
    tl::S #temperature,liquid
    pl::S #pressure,liquid
    z0::V #initial composition
    lnϕl::V #efective log fugacity coef - liquid
    lnϕv::V#efective log fugacity coef - gas
    nl::V #moles of liquid
    nv::V #moles of gas
    xil::V #mol fraction of liquid
    xiv::V #mol fraction of gas
end

function vtx_states(obj::VLEResults,moles=1.0)
    xv = obj.xiv
    xl = obj.xil
    z0 = obj.z0
    β = (z0[1] - xl[1])/(xv[1]-xl[1])
    nv = β*moles
    nl = (1-β)*moles
    stᵥ = state(mol_v=obj.vv,t=obj.tv,xn=obj.xiv,moles=nv,phase=:gas)
    stₗ = state(mol_v=obj.vl,t=obj.tl,xn=obj.xil,moles=nl,phase=:liquid)
    return (stₗ,stᵥ)
end

function vtn_states(obj::VLEResults)
    stᵥ = state(mol_v=obj.vv,t=obj.tv,n=obj.nv,phase=:gas)
    stₗ = state(mol_v=obj.vl,t=obj.tl,n=obj.nl,phase=:liquid)
    return (stₗ,stᵥ)
end

function _flash_obj(xlᵢ,xvᵢ,lnϕlᵢ,lnϕvᵢ,β)
    _0 = zero(β)
    l = ifelse(xlᵢ>_0,(1-β)*xlᵢ*(log(xlᵢ)+lnϕlᵢ),_0)
    v = ifelse(xvᵢ>_0,β*xvᵢ*(log(xvᵢ)+lnϕvᵢ),_0)
    return l+v
end

function generate_ptflash_objective(obj::VLEResults;threaded = true)
    #TODO: implement threading
    function flash_objective!(F, G, nv)
        t = obj.tv
        p = obj.pv
        β = sum(nv)
        obj.nv = nv
        #@show nv
        obj.nl =obj.z0 - obj.nv
        #corr = minimum(obj.nl)
        #@show obj.nv
        #obj.nl = obj.z0 - obj.nv - corr
        #@show obj.nl
        obj.nl = remove_minimums!!(obj.nl)
        obj.nv = remove_minimums!!(obj.nv)
        obj.xil = normalizefrac!!(obj.nl)
        obj.xiv = normalizefrac!!(obj.nv)

        #@show obj.xil
        obj.vl ,obj.lnϕl = update_logϕ!!(obj.model,obj.lnϕl,p,t,obj.xil,obj.vl,:liquid)
        obj.vv ,obj.lnϕv = update_logϕ!!(obj.model,obj.lnϕv,p,t,obj.xiv,obj.vv,:gas)
        #@show obj.lnϕl
        #@show obj.lnϕv
        if !(G === nothing)
            #@show G
            G .= zero(eltype(G))
            for i in 1:length(G)
                if obj.xiv[i] <= eps(obj.xiv[i])
                    G[i] += obj.lnϕv[i]  - log(obj.xil[i])
                elseif obj.xil[i] <= eps(obj.xil[i])
                    G[i] += log(obj.xiv[i])   - obj.lnϕl[i]
                else
                    G[i] += obj.lnϕv[i] + log(obj.xiv[i])   - obj.lnϕl[i] - log(obj.xil[i])
                #G .+= obj.lnϕv .+ log.(obj.xiv)   .- obj.lnϕl .- log.(obj.xil)
                end
            #@show G
            end
        end
        if !(F === nothing)
            fobjective = (xlᵢ,xvᵢ,lnϕlᵢ,lnϕvᵢ)->_flash_obj(xlᵢ,xvᵢ,lnϕlᵢ,lnϕvᵢ,β)
            F = mapreduce(fobjective,+,obj.xil,obj.xiv,obj.lnϕl,obj.lnϕv)
            #@show F
            return F
        end
    end
    return flash_objective!
end

function equilibria(mt::MultiPT,model::HelmholtzModel,st::ThermodynamicState)
    t = temperature(FromState(),st)
    opts = options(FromState(),st)
    thr = get(opts,:threaded,true)
    p = pressure(FromState(),st)
    z0 = mol_fraction(FromState(),st,nothing,molecular_weight(model))
    _moles = moles(FromState(),st,u"mol",molecular_weight(model))
    return flash_impl(mt,model,p,t,z0,threaded=thr,moles=_moles)
end

#v = xiv*β
#l = z0 - v

#dfugv = I/β - 1/vt + dlnfugv/vt
#dfugl = I/(1-β) - 1/lt + dlnfugl/lt

#=

water-o2 (GERG2008)
logϕl = [-0.7705794332079066, 14.528869934800472]
logϕv = [-0.009360929108700455, 0.0031516545028221]
k = [0.46709692059495184, 2.034413807942005e6]

logϕl = [-0.7705794332079066, 14.528869934800472]
logϕv = [-0.009360929108700452, 0.0031516545028221037]
k = [0.46709692059495184, 2.034413807942005e6]

logϕl = [-0.7705816360739632, 14.539131948004911]
logϕv = [-0.007992505036874167, 0.0017762880364641798]
k = [0.4664571435209862, 2.0582273486877447e6]

logϕl = [-0.770581636063973, 14.539132013140865]
logϕv = [-0.007980392574881982, 0.0017656849290196598]
k = [0.4664514936154415, 2.0582493064751942e6]

logϕl = [-0.7705816360640014, 14.53913201314133]
logϕv = [-0.007980285727848792, 0.0017655915176444346]
k = [0.4664514437764727, 2.0582494987400582e6]

logϕl = [-0.7705816360640014, 14.53913201314133]
logϕv = [-0.007980284785338597, 0.0017655906936611884]
k = [0.4664514433368374, 2.05824950043602e6]

water-o2 (VdW)
logϕl = [2.297372915993078, 7.406207150097631]
logϕv = [-0.005105694210269807, -6.038786149844145e-5]
k = [9.998935228783932, 1646.2702241435495]

logϕl = [2.297372915993078, 7.406207150097631]
logϕv = [-0.005105694210269807, -6.038786149844145e-5]
k = [9.998935228783932, 1646.2702241435495]

logϕl = [2.297372728601342, 7.406704158470847]
logϕv = [-0.003828630603820233, -0.0005244515961122089]
k = [9.986172231292917, 1647.8531690781929]

logϕl = [2.297411298402489, 7.389040025797274]
logϕv = [-0.005202637502164627, 6.585518575033793e-5]
k = [10.000288433265684, 1618.0454227403882]

logϕl = [2.297412692566457, 7.388698667718173]
logϕv = [-0.005203259181312256, 6.678678750513239e-5]
k = [10.000308592298598, 1617.4916772651154]
=#


text:
export has_cuda, has_cuda_gpu, usable_cuda_gpus

"""
    has_cuda()::Bool

Check whether the local system provides an installation of the CUDA driver and toolkit.
Use this function if your code loads packages that require CUDA, such as CuArrays.jl.

Note that CUDA-dependent packages might still fail to load if the installation is broken,
so it's recommended to guard against that and print a warning to inform the user:

```
using CUDAapi
if has_cuda()
    try
        using CuArrays
    catch ex
        @warn "CUDA is installed, but CuArrays.jl fails to load" exception=(ex,catch_backtrace())
    end
end
```
"""
function has_cuda()
    toolkit_dirs = find_toolkit()

    # check for the CUDA driver library
    libcuda = find_cuda_library("cuda", toolkit_dirs)
    return libcuda !== nothing
end

"""
    has_cuda_gpu()::Bool

Check whether the local system provides an installation of the CUDA driver and toolkit, and
if it contains a CUDA-capable GPU. See [`has_cuda`](@ref) for more details.

Note that this function initializes the CUDA API in order to check for the number of GPUs.
"""
function has_cuda_gpu()
    toolkit_dirs = find_toolkit()

    # find the CUDA driver library
    libcuda = find_cuda_library("cuda", toolkit_dirs)
    if libcuda === nothing
        return false
    end
    lib = Libdl.dlopen(libcuda)

    # initialize the API
    cuInit = Libdl.dlsym(lib, :cuInit)
    status = ccall(cuInit, Cint, (Cuint,), 0)
    if status != 0
        @warn "CUDA is installed, but fails to load (error code $status)"
        return false
    end

    # count the GPUs
    cuDeviceGetCount = Libdl.dlsym(lib, :cuDeviceGetCount)
    count_ref = Ref{Cint}()
    status = ccall(cuDeviceGetCount, Cint, (Ptr{Cint},), count_ref)
    @assert status == 0
    return count_ref[] > 0
end

"""
    usable_cuda_gpus(; suppress_output=false)::Int

Returns the number of CUDA GPUs that are available for use on the local system.

Note that this function initializes the CUDA API in order to check for the number of GPUs.
"""
function usable_cuda_gpus(; suppress_output=false)
    toolkit_dirs = find_toolkit()

    # find the CUDA driver library
    libcuda = find_cuda_library("cuda", toolkit_dirs)
    if libcuda === nothing
        return 0
    end
    lib = Libdl.dlopen(libcuda)

    # initialize the API
    cuInit = Libdl.dlsym(lib, :cuInit)
    status = ccall(cuInit, Cint, (Cuint,), 0)
    if status != 0
        !suppress_output && @warn "CUDA is installed, but fails to load (error code $status)"
        return -1
    end

    # count the GPUs
    cuDeviceGetCount = Libdl.dlsym(lib, :cuDeviceGetCount)
    count_ref = Ref{Cint}()
    status = ccall(cuDeviceGetCount, Cint, (Ptr{Cint},), count_ref)
    @assert status == 0
    return Int(count_ref[])
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -9.734460520574642 <= q <= 2.116911729869532)

                     add_NL_constraint(m, :(sigmoid(-0.8046949419530325 + -0.2278494590955038*sigmoid(0.1770285029406451 + -0.5067338360802731*sigmoid(-0.7998422425711689 + 0.5883237869246916*sigmoid(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5])) + 0.47235910506966317*sigmoid(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5]))) + -0.09889031570146667*sigmoid(-0.7622628423532047 + 0.6518912481623378*sigmoid(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5])) + 0.4006903724184485*sigmoid(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5])))) + -0.18111143847071443*sigmoid(-0.12446319977672715 + 0.7018706953017659*sigmoid(-0.7998422425711689 + 0.5883237869246916*sigmoid(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5])) + 0.47235910506966317*sigmoid(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5]))) + 0.5275006439747312*sigmoid(-0.7622628423532047 + 0.6518912481623378*sigmoid(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5])) + 0.4006903724184485*sigmoid(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5]))))) + sigmoid(-0.4025778812055263 + 0.08853333939955377*sigmoid(0.1770285029406451 + -0.5067338360802731*sigmoid(-0.7998422425711689 + 0.5883237869246916*sigmoid(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5])) + 0.47235910506966317*sigmoid(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5]))) + -0.09889031570146667*sigmoid(-0.7622628423532047 + 0.6518912481623378*sigmoid(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5])) + 0.4006903724184485*sigmoid(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5])))) + -0.4269192713141794*sigmoid(-0.12446319977672715 + 0.7018706953017659*sigmoid(-0.7998422425711689 + 0.5883237869246916*sigmoid(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5])) + 0.47235910506966317*sigmoid(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5]))) + 0.5275006439747312*sigmoid(-0.7622628423532047 + 0.6518912481623378*sigmoid(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5])) + 0.4006903724184485*sigmoid(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
#===================================================================================================
  Kernel Affinity
===================================================================================================#

doc"KernelAffinity(κ;a,c) = a⋅κ + c"
immutable KernelAffinity{T<:AbstractFloat} <: KernelOperation{T}
    a::HyperParameter{T}
    c::HyperParameter{T}
    kappa::Kernel{T}
    KernelAffinity(a::Variable{T}, c::Variable{T}, κ::Kernel{T}) = new(
        HyperParameter(a, leftbounded(zero(T), :open)),
        HyperParameter(c, leftbounded(zero(T), :closed)),
        κ
    )
end
function KernelAffinity{T<:AbstractFloat}(a::Argument{T}, c::Argument{T}, κ::Kernel{T})
    KernelAffinity{T}(Variable(a), Variable(c), κ)
end

function ==(ψ1::KernelAffinity, ψ2::KernelAffinity)
    (ψ1.a == ψ2.a) && (ψ1.c == ψ2.c) && (ψ1.kappa == ψ2.kappa)
end

ismercer(ψ::KernelAffinity) = ismercer(ψ.kappa)
isnegdef(ψ::KernelAffinity) = isnegdef(ψ.kappa)

attainszero(ψ::KernelAffinity)     = attainszero(ψ.kappa)
attainspositive(ψ::KernelAffinity) = attainspositive(ψ.kappa)
attainsnegative(ψ::KernelAffinity) = attainsnegative(ψ.kappa)

function description_string(κ::KernelAffinity, showtype::Bool = true)
    obj_str = string("KernelAffinity", showtype ? string("{", eltype(κ), "}") : "")
    kernel_str = description_string(κ.kappa)
    string(obj_str, "(a=$(κ.a.value),c=$(κ.c.value),", kernel_str, ")")
end

function convert{T<:AbstractFloat}(::Type{KernelAffinity{T}}, ψ::KernelAffinity)
    KernelAffinity(convert(T, ψ.a.value), convert(T, ψ.c.value), convert(Kernel{T}, ψ.kappa))
end

@inline phi{T<:AbstractFloat}(ψ::KernelAffinity{T}, z::T) = ψ.a*z + ψ.c

# Operations

+{T<:AbstractFloat}(κ::Kernel{T}, c::Real) = KernelAffinity(one(T), convert(T, c), κ)
+(c::Real, κ::Kernel) = +(κ, c)

*{T<:AbstractFloat}(κ::Kernel{T}, a::Real) = KernelAffinity(convert(T, a), zero(T), κ)
*(a::Real, κ::Kernel) = *(κ, a)

function +{T<:AbstractFloat}(κ::KernelAffinity{T}, c::Real)
    KernelAffinity(κ.a.value, κ.c + convert(T,c), κ.kappa)
end
+(c::Real, κ::KernelAffinity) = +(κ, c)

function *{T<:AbstractFloat}(κ::KernelAffinity{T}, a::Real)
    a = convert(T, a)
    KernelAffinity(a * κ.a, a * κ.c, κ.kappa)
end
*(a::Real, κ::KernelAffinity) = *(κ, a)

function ^{T<:AbstractFloat}(ψ::KernelAffinity{T}, d::Integer)
    KernelComposition(PolynomialClass(ψ.a.value, ψ.c.value, convert(T,d)), ψ.kappa)
end

function ^{T<:AbstractFloat}(ψ::KernelAffinity{T}, γ::AbstractFloat)
    KernelComposition(PowerClass(ψ.a.value, ψ.c.value, convert(T,γ)), ψ.kappa)
end

function exp{T<:AbstractFloat}(ψ::KernelAffinity{T})
    KernelComposition(ExponentiatedClass(ψ.a.value, ψ.c.value), ψ.kappa)
end

function tanh{T<:AbstractFloat}(ψ::KernelAffinity{T})
    KernelComposition(SigmoidClass(ψ.a.value, ψ.c.value), ψ.kappa)
end

#===================================================================================================
  Kernel Product and Sum
===================================================================================================#

# Kernel Product

immutable KernelProduct{T<:AbstractFloat} <: KernelOperation{T}
    a::HyperParameter{T}
    kappa1::Kernel{T}
    kappa2::Kernel{T}
    function KernelProduct(a::Variable{T}, κ1::Kernel{T}, κ2::Kernel)
        if !(ismercer(κ1) && ismercer(κ2))
            error("Kernels must be Mercer for closure under multiplication.")
        end
        new(HyperParameter(a, leftbounded(zero(T), :open)), κ1, κ2)
    end
end
function KernelProduct{T<:AbstractFloat}(a::Argument{T}, κ1::Kernel{T}, κ2::Kernel{T})
    KernelProduct{T}(Variable(a), κ1, κ2)
end

# Kernel Sum

immutable KernelSum{T<:AbstractFloat} <: KernelOperation{T}
    c::HyperParameter{T}
    kappa1::Kernel{T}
    kappa2::Kernel{T}
    function KernelSum(c::Variable{T}, κ1::Kernel{T}, κ2::Kernel{T})
        if !(ismercer(κ1) && ismercer(κ2)) && !(isnegdef(κ1) && isnegdef(κ2))
            error("All kernels must be Mercer or negative definite for closure under addition")
        end
        new(HyperParameter(c, leftbounded(zero(T), :closed)), κ1, κ2)
    end
end
function KernelSum{T<:AbstractFloat}(c::Argument{T}, κ1::Kernel{T}, κ2::Kernel{T})
    KernelSum{T}(Variable(c), κ1, κ2)
end

# Common Functions

for (kernel_obj, kernel_op, identity, scalar) in (
        (:KernelProduct, :*, :1, :a),
        (:KernelSum,     :+, :0, :c)
    )
    other_identity = identity == :1 ? :0 : :1
    scalar_str = string(scalar)
    @eval begin
        function description_string(ψ::$kernel_obj, showtype::Bool = true)
            constant_str = string($scalar_str,"=", ψ.$scalar.value)
            kernel1_str = "kappa1=" * description_string(ψ.kappa1, false)
            kernel2_str = "kappa2=" * description_string(ψ.kappa2, false)
            obj_str = string($kernel_obj.name.name, showtype ? string("{", eltype(ψ), "}") : "")
            string(obj_str, "(", constant_str, ",", kernel1_str, ",", kernel2_str, ")")
        end

        function ==(ψ1::$kernel_obj, ψ2::$kernel_obj)
            (ψ1.$scalar == ψ2.$scalar) && (ψ1.kappa1 == ψ2.kappa1) && (ψ1.kappa2 == ψ2.kappa2)
        end

        function convert{T<:AbstractFloat}(::Type{($kernel_obj){T}}, ψ::$kernel_obj)
            $kernel_obj(Variable(convert(T, ψ.$scalar.value), ψ.$scalar.isfixed),
                           convert(Kernel{T}, ψ.kappa1), convert(Kernel{T}, ψ.kappa2))
        end

        ismercer(ψ::$kernel_obj) = ismercer(ψ.kappa1) && ismercer(ψ.kappa2)
        isnegdef(ψ::$kernel_obj) = isnegdef(ψ.kappa1) && isnegdef(ψ.kappa2)
    end
end

for (kernel_obj, kernel_op, identity, scalar, op2_identity, op2_scalar) in (
        (:KernelProduct, :*, :1, :a, :0, :c),
        (:KernelSum,     :+, :0, :c, :1, :a)
    )
    @eval begin

        function $kernel_op{T}($scalar::Real, ψ::$kernel_obj{T})
            $kernel_obj($kernel_op(convert(T, $scalar), ψ.$scalar.value), ψ.kappa1, ψ.kappa2)
        end
        $kernel_op(ψ::$kernel_obj, $scalar::Real) = $kernel_op($scalar, ψ)

        $kernel_op{T}(κ1::Kernel{T}, κ2::Kernel{T}) = $kernel_obj(convert(T, $identity), κ1, κ2)

        function ($kernel_op){T}(κ1::KernelAffinity{T}, κ2::KernelAffinity{T})
            if κ1.$op2_scalar == $op2_identity && κ2.$op2_scalar == $op2_identity
                $kernel_obj($kernel_op(κ1.$scalar.value, κ2.$scalar.value), κ1.kappa, κ2.kappa)
            else
                $kernel_obj(convert(T, $identity), κ1, κ2)
            end
        end

        function ($kernel_op){T}(κ1::KernelAffinity{T}, κ2::StandardKernel{T})
            if κ1.$op2_scalar == $op2_identity
                $kernel_obj(κ1.$scalar.value, κ1.kappa, κ2)
            else
                $kernel_obj(convert(T, $identity), κ1, κ2)
            end
        end

        function ($kernel_op){T}(κ1::StandardKernel{T}, κ2::KernelAffinity{T})
            if κ2.$op2_scalar == $op2_identity
                $kernel_obj(κ2.$scalar.value, κ1, κ2.kappa)
            else
                $kernel_obj(convert(T, $identity), κ1, κ2)
            end
        end
    end
end


text:
### Generic rand methods

"""
    rand([rng::AbstractRNG,] s::Sampleable)

Generate one sample for `s`.

    rand([rng::AbstractRNG,] s::Sampleable, n::Int)

Generate `n` samples from `s`. The form of the returned object depends on the variate form of `s`:

- When `s` is univariate, it returns a vector of length `n`.
- When `s` is multivariate, it returns a matrix with `n` columns.
- When `s` is matrix-variate, it returns an array, where each element is a sample matrix.

    rand([rng::AbstractRNG,] s::Sampleable, dim1::Int, dim2::Int...)
    rand([rng::AbstractRNG,] s::Sampleable, dims::Dims)

Generate an array of samples from `s` whose shape is determined by the given
dimensions.
"""
rand(s::Sampleable, dims::Int...) = rand(GLOBAL_RNG, s, dims...)
rand(s::Sampleable, dims::Dims) = rand(GLOBAL_RNG, s, dims)
rand(rng::AbstractRNG, s::Sampleable, dim1::Int, moredims::Int...) =
    rand(rng, s, (dim1, moredims...))

# default fallback (redefined for univariate distributions)
function rand(rng::AbstractRNG, s::Sampleable{<:ArrayLikeVariate})
    return @inbounds rand!(rng, s, Array{eltype(s)}(undef, size(s)))
end

# multiple samples
function rand(rng::AbstractRNG, s::Sampleable{Univariate}, dims::Dims)
    out = Array{eltype(s)}(undef, dims)
    return @inbounds rand!(rng, sampler(s), out)
end
function rand(
    rng::AbstractRNG, s::Sampleable{<:ArrayLikeVariate}, dims::Dims,
)
    sz = size(s)
    ax = map(Base.OneTo, dims)
    out = [Array{eltype(s)}(undef, sz) for _ in Iterators.product(ax...)]
    return @inbounds rand!(rng, sampler(s), out, false)
end

# these are workarounds for sampleables that incorrectly base `eltype` on the parameters
function rand(rng::AbstractRNG, s::Sampleable{<:ArrayLikeVariate,Continuous})
    return @inbounds rand!(rng, sampler(s), Array{float(eltype(s))}(undef, size(s)))
end
function rand(rng::AbstractRNG, s::Sampleable{Univariate,Continuous}, dims::Dims)
    out = Array{float(eltype(s))}(undef, dims)
    return @inbounds rand!(rng, sampler(s), out)
end
function rand(
    rng::AbstractRNG, s::Sampleable{<:ArrayLikeVariate,Continuous}, dims::Dims,
)
    sz = size(s)
    ax = map(Base.OneTo, dims)
    out = [Array{float(eltype(s))}(undef, sz) for _ in Iterators.product(ax...)]
    return @inbounds rand!(rng, sampler(s), out, false)
end

"""
    rand!([rng::AbstractRNG,] s::Sampleable, A::AbstractArray)

Generate one or multiple samples from `s` to a pre-allocated array `A`. `A` should be in the
form as specified above. The rules are summarized as below:

- When `s` is univariate, `A` can be an array of arbitrary shape. Each element of `A` will
  be overriden by one sample.
- When `s` is multivariate, `A` can be a vector to store one sample, or a matrix with each
  column for a sample.
- When `s` is matrix-variate, `A` can be a matrix to store one sample, or an array of
  matrices with each element for a sample matrix.
"""
function rand! end
Base.@propagate_inbounds rand!(s::Sampleable, X::AbstractArray) = rand!(GLOBAL_RNG, s, X)
Base.@propagate_inbounds function rand!(rng::AbstractRNG, s::Sampleable, X::AbstractArray)
    return _rand!(rng, s, X)
end

# default definitions for arraylike variates
@inline function rand!(
    rng::AbstractRNG,
    s::Sampleable{ArrayLikeVariate{N}},
    x::AbstractArray{<:Real,N},
) where {N}
    @boundscheck begin
        size(x) == size(s) || throw(DimensionMismatch("inconsistent array dimensions"))
    end
    return _rand!(rng, s, x)
end

@inline function rand!(
    rng::AbstractRNG,
    s::Sampleable{ArrayLikeVariate{N}},
    x::AbstractArray{<:Real,M},
) where {N,M}
    @boundscheck begin
        M > N ||
            throw(DimensionMismatch(
                "number of dimensions of `x` ($M) must be greater than number of dimensions of `s` ($N)"
            ))
        ntuple(i -> size(x, i), Val(N)) == size(s) ||
            throw(DimensionMismatch("inconsistent array dimensions"))
    end
    # the function barrier fixes performance issues if `sampler(s)` is type unstable
    return _rand!(rng, sampler(s), x)
end

function _rand!(
    rng::AbstractRNG,
    s::Sampleable{<:ArrayLikeVariate},
    x::AbstractArray{<:Real},
)
    @inbounds for xi in eachvariate(x, variate_form(typeof(s)))
        rand!(rng, s, xi)
    end
    return x
end

Base.@propagate_inbounds function rand!(
    rng::AbstractRNG,
    s::Sampleable{ArrayLikeVariate{N}},
    x::AbstractArray{<:AbstractArray{<:Real,N}},
) where {N}
    sz = size(s)
    allocate = !all(isassigned(x, i) && size(@inbounds x[i]) == sz for i in eachindex(x))
    return rand!(rng, s, x, allocate)
end

Base.@propagate_inbounds function rand!(
    s::Sampleable{ArrayLikeVariate{N}},
    x::AbstractArray{<:AbstractArray{<:Real,N}},
    allocate::Bool,
) where {N}
    return rand!(GLOBAL_RNG, s, x, allocate)
end
@inline function rand!(
    rng::AbstractRNG,
    s::Sampleable{ArrayLikeVariate{N}},
    x::AbstractArray{<:AbstractArray{<:Real,N}},
    allocate::Bool,
) where {N}
    @boundscheck begin
        if !allocate
            sz = size(s)
            all(size(xi) == sz for xi in x) ||
                throw(DimensionMismatch("inconsistent array dimensions"))
        end
    end
    # the function barrier fixes performance issues if `sampler(s)` is type unstable
    return _rand!(rng, sampler(s), x, allocate)
end

function _rand!(
    rng::AbstractRNG,
    s::Sampleable{ArrayLikeVariate{N}},
    x::AbstractArray{<:AbstractArray{<:Real,N}},
    allocate::Bool,
) where {N}
    if allocate
        @inbounds for i in eachindex(x)
            x[i] = rand(rng, s)
        end
    else
        @inbounds for xi in x
            rand!(rng, s, xi)
        end
    end
    return x
end

"""
    sampler(d::Distribution) -> Sampleable
    sampler(s::Sampleable) -> s

Samplers can often rely on pre-computed quantities (that are not parameters
themselves) to improve efficiency. If such a sampler exists, it can be provided
with this `sampler` method, which would be used for batch sampling.
The general fallback is `sampler(d::Distribution) = d`.
"""
sampler(s::Sampleable) = s

# Random API
Random.Sampler(::Type{<:AbstractRNG}, s::Sampleable, ::Val{1}) = s
Random.Sampler(::Type{<:AbstractRNG}, s::Sampleable, ::Val{Inf}) = sampler(s)


text:
#===============================================================================

    rep_integer.jl - LR, Jan 2021

    Representation of states is done via an Integer.

===============================================================================#
const SpinState = UInt64 # That's a single species state.
const FullState = UInt128 # That's the full state describing the entire system.

# That's the shift which is applied to fuse the two spins to a full state.
const _bshift = 2^60

function f_to_s(fs::FullState)::Tuple{SpinState,SpinState}
    """ Takes a state that represents the full state in the Hilbert space
        and maps it to the constituent spin states.
    """
    return SpinState(fs%_bshift), SpinState(fs÷_bshift)
end

function s_to_f(su::SpinState, sd::SpinState)::FullState
    """ Maps two spin states to a full state.
    """
    return FullState(sd)*_bshift + FullState(su)
end

function create(s::SpinState, i::Integer)::Union{Nothing,SpinState}
    """ Destroys a particle at position i, false if there's a particle already.
    """
    shift = 1 << (i-1)
    if s & shift == 0
        return s ⊻ shift
    else
        return nothing
    end
end

function annihilate(s::SpinState, i::Integer)::Union{Nothing,SpinState}
    """ Destroys a particle at position i, false if there's no particle to destroy.
    """
    shift = 1 << (i-1)
    if s & shift > 0
        return s ⊻ shift
    else
        return nothing
    end
end

function count_occupancies(s::SpinState, src::Integer, dst::Integer)
    # return sum([(s >>> k) & 1 for k=0:63]) --> slower!
    c = 0
    for k=src-1:dst-1
        c += (s >>> k) & 1
    end
    return c
end

function sign(s::SpinState, src::Integer, dst::Integer)
    """ This is the sum of densities between the anihilation and creation operators
        (including both endpoints) which gives the entire fermionic sign factor.
        The contributions from up and down states simply add.
    """
    if  src > dst
        return (-1)^count_occupancies(s, dst, src)
    end
    return (-1)^count_occupancies(s, src, dst)
end

# Indexing.
function max_orbital(s::SpinState)
    """ Returns the index of the maximally occupied orbital in a SpinState.

        Note:
            - For integers, this is essentially the bit length.
    """
    # sizeof = number of bytes.
    # leading_zeros = built in for leading zeros in representation.
    return sizeof(s) * 8 - leading_zeros(s)
end
function Base.getindex(s::SpinState, i::Unsigned)::SpinState
    """ Note: Indexing works with the lowest bits first (i.e., the ones from the
        right). It doesn't matter how it is used, only the implementation here
        needs to be consistent.

        Only positive integers are allowed. Negatives would only give zeros.
    """
    return s >>> (i-1) & 1
end
Base.getindex(s::SpinState, i::Integer)::SpinState = s[UInt(i)]
Base.firstindex(s::SpinState) =  1
Base.lastindex(s::SpinState) = max_orbital(s)
Base.length(s::SpinState) = max_orbital(s)


text: abstract type AbstractPolicy end


text:
using Yao, Yao.Blocks
using LinearAlgebra

circuit(n) = chain(
    n,
    put(1=>X),
    repeat(H, 2:n),
    control(2, 1=>X),
    control(4, 3=>X),
    control(3, 1=>X),
    control(4, 3=>X),
    repeat(H, 1:n),
)

cc = circuit(4)

number_of_bits = nqubits(cc)

## root block
typeof(cc)
# PutBlock{4,1,XGate{Complex{Float64}},Complex{Float64}}
typeof(cc[1]).types
# svec(XGate{Complex{Float64}}, Tuple{Int64})
typeof(cc[1]).name
# PutBlock
typeof(cc[1]).parameters
# svec(4, 1, XGate{Complex{Float64}}, Complex{Float64})

#=
- [Yao.jl/BlockBasics.jl at 6af3ba2160181e592d0a0b66f14af489f4959c57 · QuantumBFS/Yao.jl](https://github.com/QuantumBFS/Yao.jl/blob/6af3ba2160181e592d0a0b66f14af489f4959c57/docs/src/tutorial/BlockBasics.jl#L75-L84)
=#

using InteractiveUtils: subtypes
function subtypetree(t, level=1, indent=4)
   level == 1 && println(t)
   for s in subtypes(t)
     println(join(fill(" ", level * indent)) * string(s))
     subtypetree(s, level+1, indent)
   end
end

subtypetree(Yao.Blocks.AbstractBlock);
#=
AbstractBlock
    AbstractMeasure
        Measure
        MeasureAndRemove
        MeasureAndReset
    FunctionBlock
    MatrixBlock
        AbstractContainer
            Concentrator
            ControlBlock
            PutBlock
            RepeatedBlock
            TagBlock
                AbstractDiff
                    BPDiff
                    QDiff
                AbstractScale
                    Scale
                    StaticScale
                CachedBlock
                Daggered
        CompositeBlock
            AddBlock
            ChainBlock
            KronBlock
            PauliString
            Roller
        PrimitiveBlock
            ConstantGate
                CNOTGate
                HGate
                I2Gate
                P0Gate
                P1Gate
                PdGate
                PuGate
                SGate
                SWAPGate
                SdagGate
                TGate
                TdagGate
                ToffoliGate
                XGate
                YGate
                ZGate
            GeneralMatrixGate
            MathBlock
            PhaseGate
            ReflectBlock
            RotationGate
            ShiftGate
            Swap
            TimeEvolution
    Sequential
=#


text:
function string_compression(str :: ASCIIString)
  """
  Update the string builder depending on
  the given character and its count
  """
  function update(prev :: Char, count :: Int)
    push!(str_builder, prev)
    append!(str_builder, map(ss->ss[1], split(string(count), "")))
  end

  @assert length(str) > 1 "[string_compression] cannot compress empty string"

  str_builder = Vector{Char}()
  prev :: Char = str[1]
  count :: Int = 0

  for ch in str
    if ch != prev         # if find a new character
      update(prev, count) # update the string builder
      prev, count = ch, 1 # reset
    else                  # if the character is the same as previous
      count += 1          # count
    end
  end
  update(prev, count)     # update the last type of character

  compressed_str = ASCIIString(str_builder)

  return length(str) > length(compressed_str) ?  compressed_str : str
end

# Time  : O(n)
# Space : O(n)

using Base: Test

@test string_compression("aabcccccaaa") == "a2b1c5a3"

max_stars_repo_name: Kiandr/crackingcodinginterview

max_issues_repo_name: Kiandr/crackingcodinginterview

max_forks_repo_name: Kiandr/crackingcodinginterview


text:
"""
    LKJCholesky(d::Int, η::Real, uplo='L')

The `LKJCholesky` distribution of size ``d`` with shape parameter ``\\eta`` is a
distribution over `LinearAlgebra.Cholesky` factorisations of ``d\\times d`` real correlation
matrices (positive-definite matrices with ones on the diagonal).

Variates or samples of the distribution are `LinearAlgebra.Cholesky` objects, as might
be returned by `F = LinearAlgebra.cholesky(R)`, so that `Matrix(F) ≈ R` is a variate or
sample of [`LKJ`](@ref).

Sampling `LKJCholesky` is faster than sampling `LKJ`, and often having the correlation
matrix in factorized form makes subsequent computations cheaper as well.

!!! note
    `LinearAlgebra.Cholesky` stores either the upper or lower Cholesky factor, related by
    `F.U == F.L'`. Both can be accessed with `F.U` and `F.L`, but if the factor
    not stored is requested, then a copy is made. The `uplo` parameter specifies whether the
    upper (`'U'`) or lower (`'L'`) Cholesky factor is stored when randomly generating
    samples. Set `uplo` to `'U'` if the upper factor is desired to avoid allocating a copy
    when calling `F.U`.

See [`LKJ`](@ref) for more details.

External links

* Lewandowski D, Kurowicka D, Joe H.
  Generating random correlation matrices based on vines and extended onion method,
  Journal of Multivariate Analysis (2009), 100(9): 1989-2001
  doi: [10.1016/j.jmva.2009.04.008](https://doi.org/10.1016/j.jmva.2009.04.008)
"""
struct LKJCholesky{T <: Real} <: Distribution{CholeskyVariate,Continuous}
    d::Int
    η::T
    uplo::Char
    logc0::T
end

#  -----------------------------------------------------------------------------
#  Constructors
#  -----------------------------------------------------------------------------

function LKJCholesky(d::Int, η::Real, _uplo::Union{Char,Symbol} = 'L'; check_args::Bool=true)
    @check_args(
        LKJCholesky,
        (d, d > 0, "matrix dimension must be positive"),
        (η, η > 0, "shape parameter must be positive"),
    )
    logc0 = lkj_logc0(d, η)
    uplo = _char_uplo(_uplo)
    T = Base.promote_eltype(η, logc0)
    return LKJCholesky(d, T(η), uplo, T(logc0))
end

# adapted from LinearAlgebra.char_uplo
function _char_uplo(uplo::Union{Symbol,Char})
    uplo ∈ (:U, 'U') && return 'U'
    uplo ∈ (:L, 'L') && return 'L'
    throw(ArgumentError("uplo argument must be either 'U' (upper) or 'L' (lower)"))
end

#  -----------------------------------------------------------------------------
#  REPL display
#  -----------------------------------------------------------------------------

Base.show(io::IO, d::LKJCholesky) = show(io, d, (:d, :η, :uplo))

#  -----------------------------------------------------------------------------
#  Conversion
#  -----------------------------------------------------------------------------

function Base.convert(::Type{LKJCholesky{T}}, d::LKJCholesky) where T <: Real
    return LKJCholesky{T}(d.d, T(d.η), d.uplo, T(d.logc0))
end
function convert(::Type{LKJCholesky{T}}, d::Integer, η::Real, uplo::Char, logc0::Real) where T <: Real
    return LKJCholesky{T}(Int(d), T(η), uplo, T(logc0))
end

#  -----------------------------------------------------------------------------
#  Properties
#  -----------------------------------------------------------------------------

Base.eltype(::Type{LKJCholesky{T}}) where {T} = T

function Base.size(d::LKJCholesky)
    p = d.d
    return (p, p)
end

function insupport(d::LKJCholesky, R::LinearAlgebra.Cholesky)
    p = d.d
    factors = R.factors
    (isreal(factors) && size(factors, 1) == p) || return false
    iinds, jinds = axes(factors)
    # check that the diagonal of U'*U or L*L' is all ones
    @inbounds if R.uplo === 'U'
        for (j, jind) in enumerate(jinds)
            col_iinds = view(iinds, 1:j)
            sum(abs2, view(factors, col_iinds, jind)) ≈ 1 || return false
        end
    else  # R.uplo === 'L'
        for (i, iind) in enumerate(iinds)
            row_jinds = view(jinds, 1:i)
            sum(abs2, view(factors, iind, row_jinds)) ≈ 1 || return false
        end
    end
    return true
end

function StatsBase.mode(d::LKJCholesky)
    factors = Matrix{eltype(d)}(LinearAlgebra.I, size(d))
    return LinearAlgebra.Cholesky(factors, d.uplo, 0)
end

StatsBase.params(d::LKJCholesky) = (d.d, d.η, d.uplo)

@inline partype(::LKJCholesky{T}) where {T <: Real} = T

#  -----------------------------------------------------------------------------
#  Evaluation
#  -----------------------------------------------------------------------------

function logkernel(d::LKJCholesky, R::LinearAlgebra.Cholesky)
    factors = R.factors
    p, η = params(d)
    c = p + 2(η - 1)
    p == 1 && return c * log(first(factors))
    # assuming D = diag(factors) with length(D) = p,
    # logp = sum(i -> (c - i) * log(D[i]), 2:p)
    logp = sum(Iterators.drop(enumerate(diagind(factors)), 1)) do (i, di)
        return (c - i) * log(factors[di])
    end
    return logp
end

function logpdf(d::LKJCholesky, R::LinearAlgebra.Cholesky)
    insupport(d, R) || throw(ArgumentError("provided point is not in the support"))
    return _logpdf(d, R)
end

_logpdf(d::LKJCholesky, R::LinearAlgebra.Cholesky) = logkernel(d, R) + d.logc0

pdf(d::LKJCholesky, R::LinearAlgebra.Cholesky) = exp(logpdf(d, R))

loglikelihood(d::LKJCholesky, R::LinearAlgebra.Cholesky) = logpdf(d, R)
function loglikelihood(d::LKJCholesky, Rs::AbstractArray{<:LinearAlgebra.Cholesky})
    return sum(R -> logpdf(d, R), Rs)
end

#  -----------------------------------------------------------------------------
#  Sampling
#  -----------------------------------------------------------------------------

function Base.rand(rng::AbstractRNG, d::LKJCholesky)
    factors = Matrix{eltype(d)}(undef, size(d))
    R = LinearAlgebra.Cholesky(factors, d.uplo, 0)
    return _lkj_cholesky_onion_sampler!(rng, d, R)
end
function Base.rand(rng::AbstractRNG, d::LKJCholesky, dims::Dims)
    p = d.d
    uplo = d.uplo
    T = eltype(d)
    TM = Matrix{T}
    Rs = Array{LinearAlgebra.Cholesky{T,TM}}(undef, dims)
    for i in eachindex(Rs)
        factors = TM(undef, p, p)
        Rs[i] = R = LinearAlgebra.Cholesky(factors, uplo, 0)
        _lkj_cholesky_onion_sampler!(rng, d, R)
    end
    return Rs
end

Random.rand!(d::LKJCholesky, R::LinearAlgebra.Cholesky) = Random.rand!(GLOBAL_RNG, d, R)
function Random.rand!(rng::AbstractRNG, d::LKJCholesky, R::LinearAlgebra.Cholesky)
    return _lkj_cholesky_onion_sampler!(rng, d, R)
end

function Random.rand!(
    rng::AbstractRNG,
    d::LKJCholesky,
    Rs::AbstractArray{<:LinearAlgebra.Cholesky{T,TM}},
    allocate::Bool,
) where {T,TM}
    p = d.d
    uplo = d.uplo
    if allocate
        for i in eachindex(Rs)
            Rs[i] = _lkj_cholesky_onion_sampler!(
                rng,
                d,
                LinearAlgebra.Cholesky(TM(undef, p, p), uplo, 0),
            )
        end
    else
        for i in eachindex(Rs)
            _lkj_cholesky_onion_sampler!(rng, d, Rs[i])
        end
    end
    return Rs
end
function Random.rand!(
    rng::AbstractRNG,
    d::LKJCholesky,
    Rs::AbstractArray{<:LinearAlgebra.Cholesky{<:Real}},
)
    allocate = any(!isassigned(Rs, i) for i in eachindex(Rs)) || any(R -> size(R, 1) != d.d, Rs)
    return Random.rand!(rng, d, Rs, allocate)
end

#
# onion method
#

function _lkj_cholesky_onion_sampler!(
    rng::AbstractRNG,
    d::LKJCholesky,
    R::LinearAlgebra.Cholesky,
)
    if R.uplo === 'U'
        _lkj_cholesky_onion_tri!(rng, R.factors, d.d, d.η, Val(:U))
    else
        _lkj_cholesky_onion_tri!(rng, R.factors, d.d, d.η, Val(:L))
    end
    return R
end

function _lkj_cholesky_onion_tri!(
    rng::AbstractRNG,
    A::AbstractMatrix,
    d::Int,
    η::Real,
    ::Val{uplo},
) where {uplo}
    # Section 3.2 in LKJ (2009 JMA)
    # reformulated to incrementally construct Cholesky factor as mentioned in Section 5
    # equivalent steps in algorithm in reference are marked.
    @assert size(A) == (d, d)
    A[1, 1] = 1
    d > 1 || return R
    β = η + (d - 2)//2
    #  1. Initialization
    w0 = 2 * rand(rng, Beta(β, β)) - 1
    @inbounds if uplo === :L
        A[2, 1] = w0
    else
        A[1, 2] = w0
    end
    @inbounds A[2, 2] = sqrt(1 - w0^2)
    #  2. Loop, each iteration k adds row/column k+1
    for k in 2:(d - 1)
        #  (a)
        β -= 1//2
        #  (b)
        y = rand(rng, Beta(k//2, β))
        #  (c)-(e)
        # w is directionally uniform vector of length √y
        @inbounds w = @views uplo === :L ? A[k + 1, 1:k] : A[1:k, k + 1]
        Random.randn!(rng, w)
        rmul!(w, sqrt(y) / norm(w))
        # normalize so new row/column has unit norm
        @inbounds A[k + 1, k + 1] = sqrt(1 - y)
    end
    #  3.
    return A
end


text:
stripc0{T<:String}(a::T) = replace(a, r"[\x00-\x1f\x7f]", "")
stripc0x{T<:String}(a::T) = replace(a, r"[^\x20-\x7e]", "")

a = "a\n\tb\u2102d\u2147f"

println("Original String:\n    ", a)
println("\nWith C0 control characters removed:\n    ", stripc0(a))
println("\nWith C0 and extended characters removed:\n    ", stripc0x(a))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode