text:
# MvNormalCDF
# Copyright © 2019-2021 Vladimir Arnautov aka PharmCat <mail@pharmcat.net>, Andrew Gough

#=
This function uses an algorithm given in the paper
"Numerical Computation of Multivariate Normal Probabilities", in
 J. of Computational and Graphical Stat., 1(1992), pp. 141-149, by
Alan Genz, WSU Math, PO Box 643113, Pullman, WA 99164-3113
Email : alangenz@wsu.edu
The primary references for the numerical integration are
"On a Number-Theoretical Integration Method"
H. Niederreiter, Aequationes Mathematicae, 8(1972), pp. 304-11, and
"Randomization of Number Theoretic Methods for Multiple Integration"
R. Cranley and T.N.L. Patterson, SIAM J Numer Anal, 13(1976), pp. 904-14.
Re-coded in Julia from the MATLAB function qsimvnv(m,r,a,b)
Alan Genz is the author the MATLAB qsimvnv() function.
Alan Genz software website: http://archive.is/jdeRh
Source code to MATLAB qsimvnv() function: http://archive.is/h5L37
% QSIMVNV(m,r,a,b) and _chlrdr(r,a,b)
%
% Copyright (C) 2013, Alan Genz,  All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided the following conditions are met:
%   1. Redistributions of source code must retain the above copyright
%      notice, this list of conditions and the following disclaimer.
%   2. Redistributions in binary form must reproduce the above copyright
%      notice, this list of conditions and the following disclaimer in
%      the documentation and/or other materials provided with the
%      distribution.
%   3. The contributor name(s) may not be used to endorse or promote
%      products derived from this software without specific prior
%      written permission.
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
% "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
% FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
% COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
% INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
% BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
% OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
% ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
% TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF USE
% OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%
=#

copy_oftype(A::AbstractArray{T}, ::Type{T}) where {T} = copy(A)
copy_oftype(A::AbstractArray{T,N}, ::Type{S}) where {T,N,S} = convert(AbstractArray{S,N}, A)

"""
    mvnormcdf(dist::MvNormal, a, b; m::Integer = 1000*size(dist.Σ,1), rng = RandomDevice())

Computes the Multivariate Normal probability integral using a quasi-Monte-Carlo
algorithm with m points for multivariate normal distributions ([MvNormal](https://juliastats.org/Distributions.jl/stable/multivariate/#Distributions.MvNormal))

Probability p is output with error estimate e.
# Arguments
- `dist::MvNormal`: multivariate normal distributions from Distributions.jl
- `a::AbstractVector`: lower integration limit column vector
- `b::AbstractVector`: upper integration limit column vector
- `m::Integer`:        number of integration points (default 1000*dimension)
- `rng`: random number generator

# Reference
- Genz, A. (1992). Numerical computation of multivariate normal probabilities. Journal of Computational and Graphical Statistics, 1, 141--150
- Genz, A. (1993). Comparison of methods for the computation of multivariate normal probabilities. Computing Science and Statistics, 25, 400--405
"""
function mvnormcdf(dist::MvNormal, a, b; m::Integer = 1000*size(dist.Σ,1), rng = RandomDevice())
    mvnormcdf(dist.μ, dist.Σ, a, b, m = m, rng = rng)
end

"""
    mvnormcdf(μ::AbstractVector{<:Real}, Σ::AbstractMatrix{<:Real}, a::AbstractVector{<:Real}, b::AbstractVector{<:Real}; m::Integer = 1000*size(Σ,1), rng = RandomDevice())

Computes the Multivariate Normal probability integral using a quasi-Monte-Carlo
algorithm with m points for positive definite covariance matrix Σ, mean [0,...], with lower
integration limit vector a and upper integration limit vector b.
```math
\\Phi_k(\\mathbf{a},\\mathbf{b},\\mathbf{\\Sigma} ) = \\frac{1}{\\sqrt{\\left | \\mathbf{\\Sigma}  \\right |{(2\\pi )^k}}}\\int_{a_1}^{b_1}\\int_{a_2}^{b_2}\\begin{align*}
 &...\\end{align*} \\int_{a_k}^{b_k}e^{^{-\\frac{1}{2}}\\mathbf{x}^t{\\mathbf{\\Sigma }}^{-1}\\boldsymbol{\\mathbf{x}}}dx_k...dx_1
```
Probability p is output with error estimate e.
# Arguments
- `μ::AbstractVector`: vector of means
- `Σ::AbstractMatrix`: positive-definite covariance matrix of MVN distribution
- `a::AbstractVector`: lower integration limit column vector
- `b::AbstractVector`: upper integration limit column vector
- `m::Integer`:        number of integration points (default 1000*dimension)
- `rng`: random number generator

# Example
```julia
Σ = [4 3 2 1; 3 5 -1 1; 2 -1 4 2; 1 1 2 5]
μ = zeros(4)
a = [-Inf; -Inf; -Inf; -Inf]
b = [1; 2; 3; 4]
m = 5000
(p,e) = mvnormcdf(μ, Σ, a, b; m=m)
#(0.605219554009911, 0.0015718064928452481)
```
Results will vary slightly from run-to-run due to the quasi-Monte-Carlo
algorithm.
"""
function mvnormcdf(Σ::AbstractMatrix{<:Real}, a::AbstractVector{<:Real}, b::AbstractVector{<:Real}; m::Integer = 1000*size(Σ,1), rng = RandomDevice())
    mvnormcdf(Zeros(size(Σ, 1)), Σ, a, b, m = m, rng = rng)
end

"""
     mvnormcdf(Σ::AbstractMatrix{<:Real}, a::AbstractVector{<:Real}, b::AbstractVector{<:Real}; m::Integer = 1000*size(Σ,1), rng = RandomDevice())

Non-central MVN distributions (with non-zero mean) can use this function by adjusting
the integration limits. Subtract the mean vector, μ, from each
integration vector.
# Example
```julia
Σ = [4 2; 2 3]
μ = [1; 2]
a = [-Inf; -Inf]
b = [2; 2]
(p,e) = mvnormcdf(Σ, a-μ, b-μ)
#(0.4306346895870772, 0.00015776288569406053)
```
"""
function mvnormcdf(μ, Σ::AbstractMatrix{<:Real}, a::AbstractVector{<:Real}, b::AbstractVector{<:Real}; m::Integer = 1000*size(Σ,1), rng = RandomDevice())
    # check for proper dimensions
    n =size(Σ, 1)
    nc=size(Σ, 2) 	# assume square Cov matrix nxn
    # check dimension > 1
    n >= 2   || throw(ErrorException("dimension of Σ must be 2 or greater. Σ dimension: $(size(Σ))"))
    n == nc  || throw(DimensionMismatch("Σ matrix must be square. Σ dimension: $(size(Σ))"))
    # check dimensions of lower vector, upper vector, and cov matrix match
    (n == length(a) == length(b)) || throw(DimensionMismatch("iconsistent argument dimensions. Sizes: Σ $(size(Σ))  a $(size(a))  b $(size(b))"))

    # check that lower integration limit a < upper integration limit b for all elements
    all(a .<= b) || throw(ArgumentError("lower integration limit a must be <= upper integration limit b"))
    # check that Σ is positive definate; if not, print warning
    #isposdef(Σ) || @warn "covariance matrix Σ fails positive definite check"
    # check if Σ, a, or b contains NaNs
    if any(isnan.(Σ)) || any(isnan.(a)) || any(isnan.(b))
    	p = NaN
    	e = NaN
    	return (p,e)
    end
    # check if a==b
    if a == b
    	p = 0.0
    	e = 0.0
    	return (p,e)
    end
    # check if a = -Inf & b = +Inf
    if all(a .== -Inf) && all(b .== Inf)
    	p = 1.0
    	e = 0.0
    	return (p,e)
    end

##################################################################
#
# Special cases: positive Orthant probabilities for 2- and
# 3-dimesional Σ have exact solutions. Integration range [0,∞]
#
##################################################################
    if all(a .== zero(eltype(a))) && all(b .== Inf) && n <= 3
    	#Σstd = sqrt.(diag(Σ))
        Σstd  = Vector{Float64}(undef, n)
        @inbounds for i in 1:n
            Σstd[i] = sqrt(Σ[i, i])
        end
    	Rcorr = cov2cor(Σ, Σstd)
    	if n == 2
    		p = 1/4 + asin(Rcorr[1,2])/(2π)
    		e = eps()
    	elseif n == 3
    		p = 1/8 + (asin(Rcorr[1,2]) + asin(Rcorr[2,3]) + asin(Rcorr[1,3]))/(4π)
    		e = eps()
    	end
    return (p,e)
    end
    a = a .- μ
    b = b .- μ
    qsimvnv!(copy_oftype(Σ, Float64), a, b, m, rng)
end

"""
Re-coded in Julia from the MATLAB function qsimvnv(m,r,a,b)
Alan Genz is the author the MATLAB qsimvnv() function.

! Mutate Σ, a, b.
"""
function qsimvnv!(Σ::AbstractMatrix, a::AbstractVector{<:Real}, b::AbstractVector{<:Real}, m::Integer, rng)
##################################################################
#
# get lower cholesky matrix and (potentially) re-ordered integration vectors
#
##################################################################
	(ch, as, bs) = _chlrdr!(Σ, a, b) # ch =lower cholesky; as=lower vec; bs=upper vec
##################################################################
#
# quasi-Monte Carlo integration of MVN integral
#
##################################################################
### setup initial values
	ai = as[1]
	bi = bs[1]
	ct = ch[1,1]
	unitnorm = Normal() # unit normal distribution
# rng = RandomDevice()
# if ai is -infinity, explicity set c=0
# implicitly, the algorith classifies anythign > 9 std. deviations as infinity
	if ai > -9*ct
		if ai < 9*ct
			c1 = cdf(unitnorm, ai/ct)
		else
			c1 = 1.0
		end
	else
		c1 = 0.0
	end
# if bi is +infinity, explicity set d=0
	if bi > -9*ct
		if bi < 9*ct
			d1 = cdf(unitnorm, bi/ct)
		else
			d1 = 1.0
		end
	else
		d1 = 0.0
	end
	n   = size(Σ,1) 	# assume square Cov matrix nxn
	cxi = c1			# initial cxi; genz uses ci but it conflicts with Lin. Alg. ci variable
	dci = d1-cxi		# initial dcxi
	p   = 0.0			# probablity = 0
	e   = 0.0			# error = 0
	# Richtmyer generators
	ps  = sqrt.(primes(Int(floor(5*n*log(n+1)/4)))) # Richtmyer generators
	q   = ps[1:n-1,1]
	ns  = 12
	nv  = Int(max(floor(m/ns),1))
    ##
	Jnv    = ones(1, nv)
	cfill  = transpose(fill(cxi, nv)) 	# evaulate at nv quasirandom points row vec
	dpfill = transpose(fill(dci, nv))
    y      = zeros(n - 1, nv)				# n-1 rows, nv columns, preset to zero
#=
Randomization loop for ns samples
	j is the number of samples to integrate over,
	but each with a vector nv in length
	i is the number of dimensions, or integrals to comptue
=#
     c  = similar(cfill)
     dc = similar(dpfill)
     pv = similar(dpfill)
     d  = similar(Jnv)
     tv = zeros(nv)'

	for j in 1:ns					# loop for ns samples
		copyto!(c, cfill)
		copyto!(dc, dpfill)
		copyto!(pv, dpfill)

		for i in 2:n
            xr = rand(rng)
            @inbounds for cnt = 1:length(c)
                x = abs(2 * mod(cnt * q[i-1] + xr, 1) - 1)
                y[i-1, cnt] = quantile(unitnorm, c[cnt] + x * dc[cnt])
            end
            s = mul!(tv , view(ch, i, 1:i-1)' , view(y, 1:i-1, :))
			ct = ch[i,i] # ch is cholesky matrix
            copyto!(c, Jnv)										# preset to 1 (>9 sd, +∞)
			copyto!(d, Jnv)										# preset to 1 (>9 sd, +∞)
            asi = as[i]
            bsi = bs[i]
            @inbounds for cnt in 1:length(c)
                aicnt = asi - s[cnt]
                bicnt = bsi - s[cnt]
                if isless(aicnt, -9*ct)
                    c[cnt] = 0.0
                elseif isless(abs(aicnt),9*ct)
                    c[cnt] = cdf(unitnorm, aicnt/ct)
                end
                if isless(bicnt, -9*ct)
                    d[cnt] = 0.0
                elseif isless(abs(bicnt),9*ct)
                    d[cnt] = cdf(unitnorm, bicnt/ct)
                end
            end
			@. dc = d - c
			@. pv = pv * dc
		end # for i
	    dm = (mean(pv)-p)/j
	    p += dm
	    e = (j-2)*e/j+dm^2
    end # for j
    e = 3*sqrt(e) 	# error estimate is 3 times standard error with ns samples
    return (p,e)  	# return probability value and error estimate
end # function qsimvnv

"""
Computes permuted lower Cholesky factor c for R which may be singular,
  also permuting integration limit vectors a and b.

! Mutate Σ, a, b

# Arguments
	Σ		matrix			Matrix for which to compute lower Cholesky matrix,
                            this is a covariance matrix
	a		vector			column vector for the lower integration limit
							algorithm may permutate this vector to improve integration
							accuracy
	b		vector			column vector for the upper integration limit
							algorithm may pertmutate this vector to improve integration
							accuracy
Output
			tuple		An a tuple with 3 returned arrays:
							1 - lower Cholesky root of r
							2 - lower integration limit (perhaps permutated)
							3 - upper integration limit (perhaps permutated)
# Examples
r = [1 0.25 0.2; 0.25 1 0.333333333; 0.2 0.333333333 1]
a = [-1; -4; -2]
b = [1; 4; 2]
(c, ap, bp) = _chlrdr(r,a,b)
result:
Lower cholesky root:
c = [ 1.00  0.0000  0.0000,
      0.20  0.9798  0.0000,
      0.25  0.2892  0.9241 ]
Permutated upper input vector:
ap = [-1, -2, -4]
Permutated lower input vector:
bp = [1, 2, 4]
"""
#############
function _chlrdr!(Σ::AbstractMatrix, a::AbstractVector, b::AbstractVector)
    # define constants
    ep = 1e-10 # singularity tolerance
    ϵ  = eps()
    # unit normal distribution
    unitnorm = Normal()
    n   = size(Σ, 1) # covariance matrix n x n square
    ckk = zero(Float64)
    dem = zero(Float64)
    am  = zero(Float64)
    bm  = zero(Float64)
    im  = zero(Float64)
    c   = Σ
    ap  = a
    bp  = b
    d   = Vector{Float64}(undef, n)
    @inbounds for i in 1:n
        d[i] = sqrt(c[i, i])
    end
    @inbounds for i in 1:n
        if d[i] > 0
            c[:,i] /= d[i]
            c[i,:] /= d[i]
            ap[i]  /= d[i]     # ap n x 1 vector
            bp[i]  /= d[i]     # bp n x 1 vector
        end
    end
    y = zeros(Float64, n) # n x 1 zero vector to start
    @inbounds for k in 1:n
        im  = k
        ckk = zero(Float64)
        dem = one(Float64)
        s   = zero(Float64)
        @inbounds for i in k:n
            if c[i,i] > ϵ  # machine error
                cii = sqrt(max(c[i,i], 0))
                if i > 1  && k > 1
                    #=
                    s = ch(i,1:i-1)*y(1:i-1);
                    seems this it multiplication of row-vector on column-vector
                    so this is dot product of two vectors in Julia,
                    not element whise multiplication
                    =#
                    #s=(c[i,1:(k-1)] .* y[1:(k-1)])[1]
                    s = dot(c[i, 1:(k-1)], y[1:(k-1)])
                end
                ai=(ap[i] - s)/cii
                bi=(bp[i] - s)/cii
                de = cdf(unitnorm, bi) - cdf(unitnorm, ai)
                if de <= dem
                    ckk = cii
                    dem = de
                    am  = ai
                    bm  = bi
                    im  = i
                end
            end # if c[i,i]> ϵ
        end # for i=
        if im > k
            c[im, im] = c[k, k]
            ap[im] , ap[k] = ap[k] , ap[im]
            bp[im] , bp[k] = bp[k] , bp[im]
            if k > 1
                c[im, 1:(k-1)], c[k,1:(k-1)] = c[k, 1:(k-1)], c[im, 1:(k-1)]
            end
            if im < n
                c[(im+1):n, im], c[(im+1):n, k] = c[(im+1):n, k], c[(im+1):n, im]
            end
            if k <= (n-1) && im <= n
                c[(k+1):(im-1), k], c[im,(k+1):(im-1)] = transpose(c[im,(k+1):(im-1)]), transpose(c[(k+1):(im-1), k])
            end
        end # if im>k
        if k < n
            c[k:k, (k+1):n] .= zero(Float64)
        end
        if ckk > k*ep
            c[k, k] = ckk
            for i in k+1:n
                c[i, k] /= ckk
                #c[i:i, (k+1):i] -= c[i,k]*transpose(c[(k+1):i,k])
                axpy!(-c[i,k], view(c, k+1:i, k), view(c, i:i, k+1:i))
            end
            if abs(dem) > ep
                y[k] = (exp(-am^2/2)-exp(-bm^2/2))/(sqrt2π*dem)
            else
                if am < -10
                    y[k] = bm
                elseif bm > 10
                    y[k] = am
                else
                    y[k] = (am + bm)/2
                end
            end # if abs
        else
            c[k:n, k] .= zero(Float64)
            y[k] = zero(Float64)
        end # if ckk>ep*k
    end # for k=
    return (c, ap, bp)
end

max_stars_repo_name: PharmCat/MvNormalCDF

max_issues_repo_name: PharmCat/MvNormalCDF

max_forks_repo_name: PharmCat/MvNormalCDF


text:
export OutgoingMessageEvent,
       MessageEvent,
       MessageAckEvent

@slackoutevent(OutgoingMessageEvent, "message",
    begin
        text::UTF8String
        channel::ChannelId
    end)

@slackevent(MessageEvent, "message",
    begin
        text::UTF8String
        channel::ChannelId
        user::UserId
        ts::EventTimestamp
    end)

@slackevent(MessageAckEvent, "message_ack",
    begin
        text::UTF8String
        channel::Nullable{UTF8String}
        ok::Bool
        ts::EventTimestamp
    end)

ext: jl


text:
# generated
module en_US # LOLTools.DataDragon
include(normpath(@__DIR__, "champion.jl"))
include(normpath(@__DIR__, "item.jl"))
include(normpath(@__DIR__, "summoner.jl"))
end # module LOLTools.DataDragon.en_US


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# ╔═╡ c63c2fbc-f1dd-11ea-3081-cd1f11630f23
using Pkg

# ╔═╡ bdbb7f1e-f1dd-11ea-1219-91d906455e4a
begin
	Pkg.add("Suppressor")
	using Suppressor
end

# ╔═╡ ca0e570a-f1dd-11ea-1655-29faf543c20a
Pkg.activate(mktempdir())

# ╔═╡ 41e94dd4-f1d9-11ea-2d2f-4b521ddb2284
md"""

# Looking under the hood

This segment introduces the tools available in Julia to **find and inspect** Julia code.

"""

# ╔═╡ 8983cd52-f1d9-11ea-1be0-c9d7d64e8942
md"""
## Finding the implementation of functions
"""

# ╔═╡ b00ae624-f1d9-11ea-003a-a149946054a9
a = 1.0 + 2.0im

# ╔═╡ c3433912-f1d9-11ea-1198-095259b8c36a
b = -2.0 + 2im

# ╔═╡ a81c0d26-f1d9-11ea-3319-19c78805c432
@which angle(a)

# ╔═╡ e779062e-f1d9-11ea-11c2-678f14dcd4a0
md"""**Tip:** infix operations like `+` are also functions.

`1 + 2` is the same as `+(1,2)`
"""

# ╔═╡ 2a66f8cc-f1da-11ea-3e26-cbe02936adb5
@which a * 2

# ╔═╡ 1b72c758-f1da-11ea-2a63-e164d7e0c63b
@which a * b

# ╔═╡ ef1f280e-fcd6-11ea-17c2-873d329119e7
md"Count the multiplies and adds in each case."

# ╔═╡ 56398582-f1da-11ea-2ae7-6f55de7e9e82
md"""
**Tip:** If you're using the default Julia REPL, you can type `@edit a * 2` and go to this line of code within your editor.

## Specialized implementations of functions

By now you may have guessed: there are many implementations of the same function and they are chosen based on the types of its arguments.

These implementations are called "**methods**".

You will see that it's very common for functions to have many methods, in fact the `*` function has a total of 364 methods!! Try to see them with `methods(*)`!
"""

# ╔═╡ 3e455938-f1da-11ea-1dc3-778cef1b6189
methods(*)

# ╔═╡ 73bef402-f1df-11ea-36d5-13e346ee1d59
md"""

Further, any other function which uses `*` will be able to make use of the correct implementation based on the types used!

For example, `prod` function which finds the product of an array uses `*` within itself.

"""

# ╔═╡ de7f0034-f1df-11ea-36e8-1d9720219a30
prod([a,b,a,b])

# ╔═╡ f0325696-f1df-11ea-2fc0-236d41409f16
prod([1,2,3,4])

# ╔═╡ dc718976-fce6-11ea-31fe-478aec584c33
prod( [rand(2,2) for i=1:4] )

# ╔═╡ f8fdbc82-fcd4-11ea-2c4f-0bd4db86b30e
md"This is a big deal, and it's called *Generic Programming* "

# ╔═╡ ba55539c-f1db-11ea-1744-41d9fc14d417
md"""

## Looking under the hood at code specialization

Julia takes a number of steps to run your code on a computer.

This involves transforming the requested Julia code into various forms of code and eventually machine code.

These are the forms your code will take in order:

- Julia code
- `@code_lowered`: Lowered Julia code ([syntactic sugar](https://en.wikipedia.org/wiki/Syntactic_sugar#:~:text=In%20computer%20science%2C%20syntactic%20sugar,style%20that%20some%20may%20prefer.)  is removed)
- `@code_typed`: Typed julia code (types are propagated depending on types of the arguments to the function)
- `@code_llvm`: Code compiled for the compiler infrastructure
- `@code_native`: Assembly code compiled for your specific machine!

"""

# ╔═╡ 7ab12188-fcd7-11ea-3d1b-db3bcb8c93e2
 #  c = 2
  c = 3.0 + 4.0im
 # c = [1 2 3]

# ╔═╡ ff22427c-fcd6-11ea-01a4-0d8fd0be3633
@code_lowered a * c

# ╔═╡ e6054270-f1dd-11ea-2c15-3f63ac42feca
md"""

**Notes on `@code_lowered` output**:

You can see that it has turned the function into a straight-line code.

- Each line does only 1 function call.
- The result of each line is assigned to a variable `%N`
"""

# ╔═╡ 32bc8b2e-f1dd-11ea-1ac2-5b2e6f2c213f
@code_typed a * c

# ╔═╡ 282d2bfe-f1de-11ea-3acc-9b4d284b99b3
md"""

**Notes on `@code_typed` output**:

Compared to `@code_lowered`, this code is more low level: it has chosen specific functions which will emit floating point instructions.
"""

# ╔═╡ 5e5f78c6-f1de-11ea-0156-25f96ef2b812
md"""
**Notes on the `@code_llvm` output**

Notice that here not only the variable names, but the field names are also gone!
"""

# ╔═╡ aa907166-f1de-11ea-00cc-a935c275eebd
md"""
**Notes on the `@code_native` output**

This is specific to your computer.

## Conclusion

Julia  allows you to look all the way to the exact instructions which will run on your CPU. Isn't that cool?

We hope you also got a feeling for the kind of specialization of code that Julia performs in order to get the fast performance for your code!

In a language without types this would not be possible! Or at least very hard!
"""

# ╔═╡ 4c131066-f1dd-11ea-3e4f-4122c5feb20a
function with_terminal(f)
	local spam_out, spam_err
	@color_output false begin
		spam_out = @capture_out begin
			spam_err = @capture_err begin
				f()
			end
		end
	end
	spam_out, spam_err

	HTML("""

$(Markdown.htmlesc(spam_out))

	""")
end

# ╔═╡ 3c9f179c-f1dd-11ea-1c63-59429d7fa23a
with_terminal() do
	@code_llvm debuginfo=:none a * c
end

# ╔═╡ dd322a8c-f1dd-11ea-245e-e54922c0f355
with_terminal() do
	@code_native a * c
end

# ╔═╡ 293c8888-fcd7-11ea-34f6-4b994aa6baf4
with_terminal() do
	@code_native a * c
end

# ╔═╡ Cell order:
# ╠═c63c2fbc-f1dd-11ea-3081-cd1f11630f23
# ╠═ca0e570a-f1dd-11ea-1655-29faf543c20a
# ╠═bdbb7f1e-f1dd-11ea-1219-91d906455e4a
# ╟─41e94dd4-f1d9-11ea-2d2f-4b521ddb2284
# ╟─8983cd52-f1d9-11ea-1be0-c9d7d64e8942
# ╠═b00ae624-f1d9-11ea-003a-a149946054a9
# ╠═c3433912-f1d9-11ea-1198-095259b8c36a
# ╠═a81c0d26-f1d9-11ea-3319-19c78805c432
# ╟─e779062e-f1d9-11ea-11c2-678f14dcd4a0
# ╠═2a66f8cc-f1da-11ea-3e26-cbe02936adb5
# ╠═1b72c758-f1da-11ea-2a63-e164d7e0c63b
# ╟─ef1f280e-fcd6-11ea-17c2-873d329119e7
# ╟─56398582-f1da-11ea-2ae7-6f55de7e9e82
# ╠═3e455938-f1da-11ea-1dc3-778cef1b6189
# ╟─73bef402-f1df-11ea-36d5-13e346ee1d59
# ╠═de7f0034-f1df-11ea-36e8-1d9720219a30
# ╠═f0325696-f1df-11ea-2fc0-236d41409f16
# ╠═dc718976-fce6-11ea-31fe-478aec584c33
# ╟─f8fdbc82-fcd4-11ea-2c4f-0bd4db86b30e
# ╟─ba55539c-f1db-11ea-1744-41d9fc14d417
# ╠═7ab12188-fcd7-11ea-3d1b-db3bcb8c93e2
# ╠═ff22427c-fcd6-11ea-01a4-0d8fd0be3633
# ╟─e6054270-f1dd-11ea-2c15-3f63ac42feca
# ╠═32bc8b2e-f1dd-11ea-1ac2-5b2e6f2c213f
# ╟─282d2bfe-f1de-11ea-3acc-9b4d284b99b3
# ╠═3c9f179c-f1dd-11ea-1c63-59429d7fa23a
# ╟─5e5f78c6-f1de-11ea-0156-25f96ef2b812
# ╠═dd322a8c-f1dd-11ea-245e-e54922c0f355
# ╠═293c8888-fcd7-11ea-34f6-4b994aa6baf4
# ╟─aa907166-f1de-11ea-00cc-a935c275eebd
# ╟─4c131066-f1dd-11ea-3e4f-4122c5feb20a

max_stars_repo_name: mathijsvdv/ComputationalThinking

max_issues_repo_name: mathijsvdv/ComputationalThinking

max_forks_repo_name: mathijsvdv/ComputationalThinking


text:
"""
Either{S,T} is a wrapper object.

# Type Parameters:
- `S`: left/right indicator.  Either `:L` or `:R`.
- `T`: specific type of either (:General, :Result)

"""
struct Either{S,T}
    value
end

# Either Monad (General type)

const LeftEither  = Either{:L, T} where T
const RightEither = Either{:R, T} where T

"""
    left(x)

Create a left-object with type LeftEither.
"""
left(x) = Either{:L, :General}(x)
left(x::LeftEither) = x
left(x::RightEither) = error("Can't make a Left out of a Right!")

"""
    right(x)
Create a right-object with type RightEither.
"""
right(x) = Either{:R, :General}(x)
right(x::RightEither) = x
right(x::LeftEither) = error("Can't make a Right out of a Left!")

left_right_value_doc = """
    left_value(x::Either{:L})
    right_value(x::Either{:R})

Extract the left (or right) underlying value from the either object `x`.
"""

"$left_right_value_doc"
left_value(x::LeftEither) where T = x.value

"$left_right_value_doc"
right_value(x::RightEither) where T = x.value

"""
    fmap(f::Function, x::Either)

Map function `f` over the unwrapped value of `x`.
"""
fmap(f::Function, x::LeftEither) = left(f(left_value(x)))
fmap(f::Function, x::RightEither) = right(f(right_value(x)))

"Return true if `x` is a left-object."
is_left(x::Either{S,T}) where {S,T} = S === :L

"Return true if `x` is a right-object."
is_right(x::Either{S,T}) where {S,T} = S === :R

Base.show(io::IO, x::LeftEither) =
    print(io, "MonadEither_Left(", left_value(x), ")")

Base.show(io::IO, x::RightEither) =
    print(io, "MonadEither_Right(", right_value(x), ")")

# Result Monad

const ResultEither = Either{:R, :Result}
const ErrorEither = Either{:L, :Result}

"""
    result(x)

Create a Result monad, which is either a left-object or right-object.
Both `Exception`, `Nothing`, and `None` are considered left.
Everything else is right.
"""
result(x) = ResultEither(x)
result(x::Exception) = ErrorEither(x)

fmap(f::Function, x::ResultEither) = result(f(right_value(x)))
fmap(f::Function, x::ErrorEither) = x

Base.show(io::IO, x::ResultEither) where {S} =
    print(io, "MonadResult_Value(", right_value(x), ")")

Base.show(io::IO, x::ErrorEither) where {S} =
    print(io, "MonadResult_Error(", left_value(x), ")")


text:
# First basic plot function
generatecolors() = [:red, :blue, :green, :yellow, :white, :cyan, :magenta]

# Lineplots TODO: Add LearnBase.ObsDim support so that one can either plot samples or variables (so far only variable plotting supported)
lineplot(xidx::Int=1; color=:white, kwargs...) = FunctionCell(lineplot, (xidx,), "Line Plot (xidx=$xidx)"; color=color, kwargs...)
lineplot(x::T where T<:CellData, xidx::Int=1; color=:white, kwargs...) = lineplot(getobs(_variable_(x,xidx)); color=color, kwargs...)
lineplot(x::T where T<:AbstractArray, xidx::Int=1; color=:white, kwargs...) = UnicodePlots.lineplot(getobs(_variable_(x,xidx)); color=color, kwargs...)
lineplot(f::Function, xidx::Int=1; color=:white, kwargs...) = FunctionCell(x->UnicodePlots.lineplot(f, getobs(_variable_(x,xidx)); color=color, kwargs...)
									   , (), "Line Plot (f=$f, xidx=$xidx)")
lineplot(x::T where T<:CellData, f::Function, xidx::Int=1; color=:white, kwargs...) = UnicodePlots.lineplot(f, getobs(_variable_(x,xidx)); color=color, kwargs...)

#=
# Cumbersome for vector of functions; TODO: UnicodePlots.lineplot! should be used by superimposing the plots of vf[i]( getobs(_variable_(x,xidx)) )
lineplot(vf::Vector{Function}, idx::Int=1; kwargs...) = FunctionCell(x->UnicodePlots.lineplot([y->f(y) for f in vf], minimum(x), maximum(x); kwargs...)
								, (), "Line Plot (vf=$(vf), idx=$idx")
lineplot{T <: CellData}(x::T, vf::Vector{Function}, idx::Int=1; kwargs...) = begin
	data = getobs(_variable_(x,idx))
	UnicodePlots.lineplot([, data, minimum(data), maximum(data); kwargs...)
end
=#

# Scatterplots
scatterplot(xidx::Int=1, yidx::Int=2; kwargs...) = j4pr.FunctionCell(scatterplot, (xidx, yidx), "Scatter Plot (xidx=$xidx, yidx=$yidx)"; kwargs...)
scatterplot(dc::T where T<:CellDataU, xidx::Int=1, yidx::Int=2; color=:white, kwargs...) = scatterplot(getx!(dc), xidx, yidx; color=color, kwargs...)
scatterplot(dc::T where T<:CellDataL, xidx::Int=1, yidx::Int=2; kwargs...) = scatterplot(strip(dc), xidx, yidx; kwargs...)

scatterplot(x::T where T<:AbstractMatrix, xidx::Int=1, yidx::Int=2; color=:white, kwargs...) = begin

	@assert nvars(x) > 1 "[scatterplot] Dimensionality of the data has to be larger than 1."

	x1 = getobs(_variable_(x,xidx))
	x2 = getobs(_variable_(x,yidx))
	maxx1 = maximum(x1)
	minx1 = minimum(x1)
	maxx2 = maximum(x2)
	minx2 = minimum(x2)

	# Plot
	p = UnicodePlots.scatterplot(x1, x2, xlim=[minx1,maxx1], ylim=[minx2, maxx2]; color=color, kwargs...)
end

scatterplot(x::Tuple{T,S} where T<:AbstractArray where S<:AbstractVector, xidx::Int=1, yidx::Int=2; kwargs...) = begin

	@assert nvars(x[1]) > 1 "[scatterplot] Dimensionality of the data has to be larger than 1."

	# Read needed data
	data = getobs(x[1])
	labels = getobs(x[2])
	classes = sort(unique(labels))

	# Plot 1-st class
	x1 = getobs(_variable_(data,xidx))
	x2 = getobs(_variable_(data,yidx))

	colors = generatecolors()
	cmask = labels .== classes[1]
	x1c = x1[cmask]
	x2c = x2[cmask]

	maxx1 = maximum(x1)
	minx1 = minimum(x1)
	maxx2 = maximum(x2)
	minx2 = minimum(x2)

	p = UnicodePlots.scatterplot(x1c, x2c, color = colors[1], xlim=[minx1,maxx1], ylim=[minx2, maxx2]; kwargs...)

	# Plot the rest of the classes
	for i in 2:length(classes)
		cmask = labels .== classes[i]
		x1c = getobs(_variable_(data,xidx))[cmask]
		x2c = getobs(_variable_(data,yidx))[cmask]
		UnicodePlots.scatterplot!(p, x1c, x2c, color = colors[(i-1)%length(colors)+1])
 	end
	p
end

# 2-D density plots
densityplot2d(xidx::Int=1, yidx::Int=2; kwargs...) = j4pr.FunctionCell(densityplot2d, (xidx, yidx), "Density Plot 2D (xidx=$xidx, yidx=$yidx)"; kwargs...)
densityplot2d(dc::T where T<:CellDataU, xidx::Int=1, yidx::Int=2; color=:white, kwargs...) = densityplot2d(getx!(dc), xidx, yidx; color=color, kwargs...)
densityplot2d(dc::T where T<:CellDataL, xidx::Int=1, yidx::Int=2; kwargs...) = densityplot2d(strip(dc), xidx, yidx; kwargs...)

densityplot2d(x::T where T<:AbstractMatrix, xidx::Int=1, yidx::Int=2; color=:white, kwargs...) = begin

	@assert nvars(x) > 1 "[scatterplot] Dimensionality of the data has to be larger than 1."

	x1 = getobs(_variable_(x,xidx))
	x2 = getobs(_variable_(x,yidx))
	maxx1 = maximum(x1)
	minx1 = minimum(x1)
	maxx2 = maximum(x2)
	minx2 = minimum(x2)

	# Plot
	p = UnicodePlots.densityplot(x1, x2, xlim=[minx1,maxx1], ylim=[minx2, maxx2]; color=color, kwargs...)
end

densityplot2d(x::Tuple{T,S} where T<:AbstractArray where S<:AbstractVector, xidx::Int=1, yidx::Int=2; kwargs...) = begin

	@assert nvars(x[1]) > 1 "[scatterplot] Dimensionality of the data has to be larger than 1."

	# Read needed data
	data = getobs(x[1])
	labels = getobs(x[2])
	classes = sort(unique(labels))

	# Plot 1-st class
	x1 = getobs(_variable_(data,xidx))
	x2 = getobs(_variable_(data,yidx))

	colors = generatecolors()
	cmask = labels .== classes[1]
	x1c = x1[cmask]
	x2c = x2[cmask]

	maxx1 = maximum(x1)
	minx1 = minimum(x1)
	maxx2 = maximum(x2)
	minx2 = minimum(x2)

	p = UnicodePlots.densityplot(x1c, x2c, color = colors[1], xlim=[minx1,maxx1], ylim=[minx2, maxx2]; kwargs...)

	# Plot the rest of the classes
	for i in 2:length(classes)
		cmask = labels .== classes[i]
		x1c = getobs(_variable_(data,xidx))[cmask]
		x2c = getobs(_variable_(data,yidx))[cmask]
		UnicodePlots.densityplot!(p, x1c, x2c, color = colors[(i-1)%length(colors)+1])
 	end
	p

end

# Histogram ... of sorts, using lineplot
densityplot1d(xidx::Int=1; nbins=100, kwargs...) = j4pr.FunctionCell(densityplot1d, (xidx,), "Density Plot (xidx=$xidx, nbins=$nbins)"; nbins=nbins, kwargs...)
densityplot1d(dc::T where T<:CellDataU, xidx::Int=1; nbins=100, color=:white, kwargs...) = densityplot1d(getx!(dc), xidx; nbins=nbins, color=color, kwargs...)
densityplot1d(dc::T where T<:CellDataL, xidx::Int=1; nbins=100, kwargs...) = densityplot1d(strip(dc), xidx; nbins=nbins, kwargs...)
densityplot1d(x::T where T<:AbstractMatrix, xidx::Int=1; nbins=100, color=:white, kwargs...) = densityplot1d(getobs(_variable_(x,xidx)), xidx; nbins=nbins, color=color, kwargs...)

densityplot1d(x::T where T<:AbstractVector, xidx::Int=1; nbins=100, color=:white, kwargs...) = begin

	h = fit(Histogram, x, closed=:left, nbins=nbins)
	r = h.edges[1]
	x1 = r[1:length(r)-1].+0.5*step(r)
	x2 = h.weights

	p = UnicodePlots.lineplot(x1,x2; color=color, kwargs...)
end

densityplot1d(x::Tuple{T,S} where T<:AbstractArray where S<:AbstractVector, xidx::Int=1; nbins=nbins, kwargs...) = begin

	# Read needed data
	data = getobs(_variable_(x[1],xidx))
	labels = getobs(x[2])
	classes = sort(unique(labels))

	colors = generatecolors()
	cmask = labels .== classes[1]
	maxx = maximum(data)
	minx = minimum(data)
	miny = 0
	maxy = 0

	# search first to find maximum height for plot
	for c in classes
		h = fit(Histogram, data[labels.==c], nbins=nbins, closed=:left)
		maxy = max(maxy,maximum(h.weights))
	end

	# Plot 1-st class
	h = fit(Histogram, data[cmask], nbins=nbins, closed=:left)
	r = h.edges[1]
	x = r[1:length(r)-1].+0.5*step(r)
	y = h.weights
	maxy = max(maxy,maximum(y))
	p = UnicodePlots.lineplot(x, y, color=colors[1], xlim=[minx,maxx], ylim=[miny, maxy]; kwargs...)

	# Plot the rest of the classes
	for i in 2:length(classes)
		cmask = labels .== classes[i]
		h = fit(Histogram, data[cmask], nbins=nbins, closed=:left)
		r = h.edges[1]
		x = r[1:length(r)-1].+0.5*step(r)
		y = h.weights
		maxy = max(maxy,maximum(y))
		UnicodePlots.lineplot!(p, x, y, color = colors[(i-1)%length(colors)+1], xlim=[minx,maxx],ylim=[miny,maxy] ; kwargs...)
 	end
	p
end

# ROC curve plots
if v"0.6" <= VERSION < v"0.7-"

using j4pr.ROC: ComplexOP, AbstractPerfMetric, TPr, FPr, TNr, FNr

rocplot(x::CellFunT{<:Model{<:ComplexOP}}, xmetric::AbstractPerfMetric=FPr(),
		ymetric::AbstractPerfMetric=TPr(); color=:white, kwargs...) =
	rocplot(getx!(x), xmetric, ymetric; color=color, kwargs...)

rocplot(model::Model{<:ComplexOP}, xmetric::AbstractPerfMetric=FPr(),
		ymetric::AbstractPerfMetric=TPr(); color=:white, kwargs...) =
	rocplot(model.data, xmetric, ymetric; color=color, kwargs...)

rocplot(op::ComplexOP, xmetric::AbstractPerfMetric=FPr(), ymetric::AbstractPerfMetric=TPr(); color=:white, kwargs...)=
begin

	_get_data_(::TPr, op::ComplexOP) = 1-op.rocdata[:,3] 	# TPr
	_get_data_(::FNr, op::ComplexOP) = op.rocdata[:,3]  	# FNr
	_get_data_(::FPr, op::ComplexOP) = op.rocdata[:,2]   	# FPr
	_get_data_(::TNr, op::ComplexOP) = 1-op.rocdata[:,2]  	# TNr

	# Get data from op
	x = _get_data_(xmetric, op)
	y = _get_data_(ymetric, op)

	p = UnicodePlots.lineplot(x, y; color=color, kwargs...)
end
end


text:
# Sample script for plotting the phase space density near a given spatial location.
#
# Hongyang Zhou, hyzhou@umich.edu

using Vlasiator, PyPlot
using Vlasiator: RE # Earth radius [m]

file = "bulk.0001347.vlsv"

meta = load(file)

coordinates = [0.0, 0.0, 0.0]

vdfslice(meta, coordinates; verbose=true)

# Show the spatial distribution of cells with saved VDF
cellswithVDF = getcellwithvdf(meta)
locations = [getcellcoordinates(meta, cid) for cid in cellswithVDF]

xcell = zeros(size(locations))
ycell = similar(xcell)

for i in eachindex(locations)
   xcell[i] = locations[i][1] / RE
   ycell[i] = locations[i][2] / RE
end

figure()
pcolormesh(meta, "vg_pressure", colorscale=Linear)
ax = plt.gca()
ax.scatter(xcell, ycell, marker="+", color="w")


text:
using Parameters

"""

Rank a vector.

$(SIGNATURES)

## Positional arguments
```julia
* `x::Vector` # Vector of the draws
* `nt_args::NamedTuple` # NamedTuple specifying n_draws and n_chains in x
```

Not exported.

"""
function rank_vector(x::Vector{Float64}, nt_args::NamedTuple)
    @unpack n_draws, n_chains = nt_args
    xc = zeros(Int, n_draws, n_chains)
    ranked = Int.(ordinalrank(x))
    for j in 1:n_chains
        start_draw = (j-1) * n_draws + 1
        last_draw = j * n_draws
        xc[:, j] = ranked[start_draw:last_draw]
    end
    xc
end


text:
export switch_map_to

"""
    switch_map_to(inner_observable)

Creates a `switch_map_to` operator, which projects each source value to the same Observable which is
flattened multiple times with `switch_map` in the output Observable.

# Arguments
- `inner_observable`: an Observable to replace each value from the source Observable.

# Producing

Stream of type `<: Subscribable{R}` where R refers to the eltype of `inner_observable`

# Examples
```jldoctest
using Rocket

source = from([ 0, 0, 0 ])
subscribe!(source |> switch_map_to(from([ 1, 2, 3 ])), logger())
;

# output

[LogActor] Data: 1
[LogActor] Data: 2
[LogActor] Data: 3
[LogActor] Data: 1
[LogActor] Data: 2
[LogActor] Data: 3
[LogActor] Data: 1
[LogActor] Data: 2
[LogActor] Data: 3
[LogActor] Completed
```

See also: [`switch_map`](@ref), [`AbstractOperator`](@ref), [`RightTypedOperator`](@ref), [`ProxyObservable`](@ref), [`logger`](@ref)
"""
switch_map_to(source::S) where S = as_switch_map_to(as_subscribable(S), source)

as_switch_map_to(::InvalidSubscribableTrait,  source)             = throw(InvalidSubscribableTraitUsageError(source))
as_switch_map_to(::SimpleSubscribableTrait{R}, source)    where R = switch_map(R, (_) -> source)
as_switch_map_to(::ScheduledSubscribableTrait{R}, source) where R = switch_map(R, (_) -> source)


text:
# using LinearAlgebra, SparseArrays

function svec(U)
    # Get the size
    n = size(U, 1)
    N = n*(n+1)÷2
    # Preallocate
    res = zeros(Float64, N)
    index = 1
    # Iterate in column-major order
    for i = 1:n
        for j = 1:i
            if i == j
                res[index] = U[j,i]
            else
                res[index] = sqrt(2.0)*U[j,i]
            end
            index += 1
        end
    end
    return res
end

function svec_lower(U)
    # Get the size
    n = size(U, 1)
    N = n*(n+1)÷2
    # Preallocate
    res = zeros(Float64, N)
    index = 1
    # Iterate in column-major order
    for i = 1:n
        for j = 1:i
            if i == j
                res[index] = U[j,i]
            else
                res[index] = sqrt(2.0)*U[i,j]
            end
            index += 1
        end
    end
    return res
end
function smat(u)
    # Get the size
    N = length(u)
    n = (isqrt(1+8*N)-1)÷2
    # Preallocate
    res = zeros(Float64, n, n)
    index = 1
    for j = 1:n
        for i = 1:j
            if i==j
                res[i,j] = u[index]
            else
                res[i,j] = u[index]/sqrt(2.0)
                res[j,i] = u[index]/sqrt(2.0)
            end
            index += 1
        end
    end
    return res
end

function readSDPA(filename)
    # Open file
    io = open(filename, "r")
    # Read the whole file into a string
    content = read(io, String)
    # Close file
    close(io)
    # Parse the string
    lines = split(content, "\n")
    nlines = length(lines)
    # Number of constraints
    m = parse(Int, lines[1])

    # Number of blocks
    nblocks = parse(Int, lines[2])

    # Number of variables
    blocks = split(lines[3], " ")
    blocksizes = zeros(Int, nblocks)
    blockstarts = zeros(Int, nblocks)
    blockstarts[1] = 1
    n = 0
    for i in 1:nblocks
        n += parse(Int, blocks[i])
        blocksizes[i] = parse(Int, blocks[i])
    end

    for i in 2:nblocks
        blockstarts[i] = blockstarts[i-1] + blocksizes[i]
    end

    # Allocate arrays
    b = zeros(Float64, m)
    bvals = split(lines[4], " ")
    for j in 1:m
        b[j] = -parse(Float64, bvals[j])
    end
    C = zeros(Float64, n, n)
    A = zeros(Float64, n, m*n)

    for l = 5:nlines-1 # \n at the end
        line = split(lines[l], " ")
        cstr_num = parse(Int, line[1])
        block_num = parse(Int, line[2])
        i = parse(Int, line[3])
        j = parse(Int, line[4])
        val = parse(Float64, line[5])
        if cstr_num == 0
            i_index = blockstarts[block_num] + i -1
            j_index = blockstarts[block_num] + j -1
            C[i_index, j_index] = -val
            C[j_index, i_index] = -val
        else
            i_index = blockstarts[block_num] + i -1
            j_index = blockstarts[block_num] + j -1
            A[i_index, (cstr_num-1)*n + j_index] = -val
            A[j_index, (cstr_num-1)*n + i_index] = -val
        end
    end

    return A, b, C
end

function smat!(U, u)
    # Get the size
    n = size(U, 1)
    # Preallocate
    index = 1
    @views for j = 1:n
        for i = 1:j
            if i==j
                U[i,j] = u[index]
            else
                U[i,j] = u[index]/sqrt(2.0)
                U[j,i] = u[index]/sqrt(2.0)
            end
            index += 1
        end
    end
end

function fast_smat!(U, u)
    # Copies a vectorized symmetric matrix u into a full matrix U.

    # Get the size of the matrix and precompute constants.
    n = size(U, 1)
    inv_sqr = 1.0/sqrt(2.0)

    # Computes the index of column in u.
    @inline col_start(col::Int) = (col-1)*col÷2+1

    # Iterate through columns of U.
    @inbounds for k = 1:n
        # Fill the part of column k in the upper triangle of U.
        c_start = col_start(k)
        for j = c_start:c_start+k-2
            U[j-c_start+1, k] = inv_sqr*u[j]
        end
        # Fill the diagonal entry.
        U[k, k] = u[col_start(k)+k-1]
        # Fill the remaining part in the lower triangle.
        for j = k+1:n
            U[j, k] = inv_sqr*u[col_start(j)+k-1]
        end
    end
end
#####################################################################
############## In-place symmetric matrices functions #####################
##########################################################################

function svec!(u, U)

    N = size(U, 1)
    index = 1
    # Iterate in column-major order
    for i = 1:N
        for j = 1:i
            if i == j
                @inbounds u[index] = U[j,i]
            else
                @inbounds u[index] = sqrt(2.0)*U[j,i]
            end
            index += 1
        end
    end
end

function svec_add!(u, U)
    N = size(U, 1)
    n = N*(N+1)÷2
    index = 1
    # Iterate in column-major order
    for i = 1:N
        for j = 1:i
            if i == j
                u[index] += U[j,i]
            else
                u[index] += sqrt(2.0)*U[j,i]
            end
            index += 1
        end
    end
end

function copy_mat!(C, R)
    n = size(R, 1)
    for i = 1:n
        for j = 1:i
            C[j,i] = R[j,i]
        end
    end
    triu!(C)
end

function add_transpose!(U)
    n = size(U, 1)
    for i = 1:n
        for j = 1:n
            @inbounds U[i,j] += U[j,i]
            @inbounds U[i,j] *= 0.5
        end
    end
end

function ind2cart(i)
    m = (isqrt(8*i+1)-1)÷2
    t1 = m*(m+1)÷2
    t2 = (m*m + 3*m +2)÷2
    if t1 == i
        t = t1
    else
        t = max(t1,t2)
    end
    y = isqrt(2*t)
    x = y - (t-i)
    # println(x,", ",y)
    return x, y
end

function sparse_A_mul(C, a, B, colptr, rowval, nzval)
    n = size(B, 1)
    NNZ = length(a.nzind)

    fill!(colptr, zero(BlasInt))
    fill!(rowval, zero(BlasInt))
    fill!(nzval, zero(Float64))

    for k = 1:length(a.nzind)
        i, j = ind2cart(a.nzind[k])
        colptr[j+1] += 1
        rowval[k] = i
        nzval[k] = i == j ? a.nzval[k] : a.nzval[k]/sqrt(2.0)
    end

    colptr[1] = 1
    for k = 1:n
        colptr[k+1] += colptr[k]
    end

    mat = SparseMatrixCSC{Float64, BlasInt}(n, n, colptr, rowval, nzval)
    # display(colptr)
    MKLSparse.BLAS.cscmm!('N', 1.0, "SUNF", mat, B, 0.0, C)
end

max_stars_repo_name: TomaszWasTaken/MySDPSolver

max_issues_repo_name: TomaszWasTaken/MySDPSolver

max_forks_repo_name: TomaszWasTaken/MySDPSolver


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

include("../../src/Sparso.jl")
include("../../src/simple-show.jl")
using Sparso

set_options(SA_ENABLE, SA_VERBOSE, SA_USE_SPMP, SA_REORDER)

function foo(A, x)
    A*x
end

m = 10
A = sprand(m, m, 0.1)
x = repmat([1/m], m)
@acc foo(A, x)

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
module Auth

using StructTypes, Dates, JWTs, HTTP
using ..DataModel

const JWT_AUTH_KEYS = Ref{JWKSet}()

function init(authkeysfile)
    JWT_AUTH_KEYS[] = JWKSet(authkeysfile)
    refresh!(JWT_AUTH_KEYS[])
    return
end

# generating a new token for a user
const DATE_FORMAT = DateFormat("e, dd u yyyy HH:MM:SS G\\MT") # Wed, 21 Oct 2015 07:28:00 GMT
const JWT_TOKEN_COOKIE_NAME = "X-P2Pai-Jwt-Token"

function addtoken!(resp::HTTP.Response, user::User)
    exp = Dates.now(Dates.UTC) + Dates.Hour(12)
    payload = Dict("iss" => "MusicAlbums.jl", "exp" => Dates.datetime2unix(exp), "sub" => "managing albums", "aud" => user.username, "uid" => user.id)
    jwt = JWT(; payload = payload)
    keyid = first(first(JWT_AUTH_KEYS[].keys))
    sign!(jwt, JWT_AUTH_KEYS[], keyid)
    HTTP.setheader(resp, "Set-Cookie" => "$JWT_TOKEN_COOKIE_NAME=$(join([jwt.header, jwt.payload, jwt.signature], '.')); Expires=$(Dates.format(exp, DATE_FORMAT))")
    HTTP.setheader(resp, JWT_TOKEN_COOKIE_NAME => join([jwt.header, jwt.payload, jwt.signature], '.'))
    return resp
end

struct Unauthenticated <: Exception end

# parsing token from HTTP request
function User(req::HTTP.Request)
    if HTTP.hasheader(req, "Cookie")
        cookies = filter(x -> x.name == JWT_TOKEN_COOKIE_NAME, HTTP.cookies(req))
        if !isempty(cookies) && !isempty(cookies[1].value)
            jwt = JWT(; jwt = cookies[1].value)
            verified = false
            for kid in JWT_AUTH_KEYS[].keys
                validate!(jwt, JWT_AUTH_KEYS[], kid[1])
                verified |= isverified(jwt)
            end
            if verified
                parts = claims(jwt)
                return User(parts["uid"], parts["aud"])
            end
        end
    elseif HTTP.hasheader(req, JWT_TOKEN_COOKIE_NAME)
        jwt = JWT(; jwt = String(HTTP.header(req, JWT_TOKEN_COOKIE_NAME)))
        verified = false
        for kid in JWT_AUTH_KEYS[].keys
            validate!(jwt, JWT_AUTH_KEYS[], kid[1])
            verified |= isverified(jwt)
        end
        if verified
            parts = claims(jwt)
            return User(parts["uid"], parts["aud"])
        end
    end
    throw(Unauthenticated())
end

end # module


text:
module MLJExampleInterface

# Example.jl doesn't actually provide machine learning models, so we
# provide a module below with the same name to furnish us with simple
# constant probabilistic classification. Note that `Example.fit`
# ignores the training features `Xmatrix`.

module Example

function fit(Xmatrix::Matrix, yint::AbstractVector{<:Integer})
    classes = sort(unique(yint))
    counts = [count(==(c), yint) for c in classes]
    Θ = counts / sum(counts)
end

predict(Xnew::Matrix, Θ) = vcat(fill(Θ', size(Xnew, 1))...)

# julia> yint = rand([1,3,4], 100);

# julia> Θ = fit(rand(100, 3), yint)
# 3-element Vector{Float64}:
#  0.35
#  0.23
#  0.42

# julia> predict(rand(5, 3), Θ)
# 5×3 Matrix{Float64}:
#  0.35  0.23  0.42
#  0.35  0.23  0.42
#  0.35  0.23  0.42
#  0.35  0.23  0.42
#  0.35  0.23  0.42

end # of module

### CONTINUATION OF TEMPLATE

import .Example # substitute model-providing package name here (no dot)
import MLJModelInterface
import ScientificTypesBase

const PKG = "Example"          # substitute model-providing package name
const MMI = MLJModelInterface
const STB = ScientificTypesBase

"""
    CoolProbabilisticClassifier()

A cool classifier that predicts `UnivariateFinite` probability
distributions. These are distributions for a finite sample space whose
elements are *labeled*.

"""
MMI.@mlj_model mutable struct CoolProbabilisticClassifier <: MMI.Probabilistic
    dummy_hyperparameter1::Float64 = 1.0::(_ ≥ 0)
    dummy_hyperparameter2::Int = 1::(0 < _ ≤ 1)
    dummy_hyperparameter3
end

function MMI.fit(::CoolProbabilisticClassifier, verbosity, X, y)

    Xmatrix = MMI.matrix(X)

    yint = MMI.int(y)
    decode = MMI.decoder(y[1])                # for decoding int repr.
    classes_seen = decode(sort(unique(yint))) # ordered by int repr.

    Θ = Example.fit(Xmatrix, yint)            # probability vector
    fitresult = (Θ, classes_seen)
    report = (n_classes_seen = length(classes_seen),)
    cache = nothing

    return fitresult, cache, report

end

function MMI.predict(::CoolProbabilisticClassifier, fitresult, Xnew)
    Xmatrix = MMI.matrix(Xnew)

    Θ, classes_seen = fitresult
    prob_matrix = Example.predict(Xmatrix, Θ)

    # `classes_seen` is a categorical vector whose pool actually
    # includes *all* classes. The `UnivariateFinite` constructor
    # automatically assigns zero probability to the unseen classes.

    return MMI.UnivariateFinite(classes_seen, prob_matrix)
end

# for returning user-friendly form of the learned parameters:
function MMI.fitted_params(::CoolProbabilisticClassifier, fitresult)
    Θ, classes_seen = fitresult
    return (raw_probabilities = Θ, classes_seen_in_training = classes_seen)
end

## META DATA

MMI.metadata_pkg(CoolProbabilisticClassifier,
             name="$PKG",
             uuid="7876af07-990d-54b4-ab0e-23690620f79a",
             url="https://github.com/JuliaLang/Example.jl",
             is_pure_julia=true,
             license="MIT",
             is_wrapper=false
)

MMI.metadata_model(CoolProbabilisticClassifier,
               input_scitype = MMI.Table(STB.Continuous),
               target_scitype = AbstractVector{<:STB.Finite},# ie, a classifier
               docstring = "Really cool classifier",         # brief description
               path = "$PKG.CoolProbabilisiticClassifier"
               )

end # module


text:
using ForwardDiff
using Distances
using LinearAlgebra
using ReverseDiff
using AutoGrad
#using CuArrays
using Flux
using Flux.Tracker
using Flux: @epochs

LR = 0.1 # learning rate
EPOCHS = 100 # number of epochs
DIM = 2 # dimensionality of embeddings

function radius(x)
    return x[end]
end

function centerpoint(x)
    return x[1:end-1]
end

function loss1(c, d) # loss for normal form 1
    return max(0, euclidean(centerpoint(c), centerpoint(d)) + abs(radius(c))-abs(radius(d)))
end

function loss2(c1, c2, d) # loss for c1 and c2 SubClassOf: d
    dist = euclidean(centerpoint(c1), centerpoint(c2))
    if dist > abs(radius(c1) + radius(c2))
        return abs(radius(c2) + radius(c2)) # no solution, circles are separate
    elseif dist < abs(radius(c1)-radius(c2)) # no solution, one circle contained in other
        if radius(c1) < radius(c2)
            return loss1(c1, d)
        else
            return loss1(c2, d)
        end
    elseif dist == 0 && radius(c1) == radius(c2) ## circles are coincident
        return loss1(c1, d)
    else
        a = (radius(c1)^2 - radius(c2)^2 + dist^2)/(2 * dist)
        rad = sqrt(radius(c1)^2 - a^2)
        cent = centerpoint(c1) + a * (centerpoint(c2) - centerpoint(c1)) / dist
        return loss1(vcat(cent, rad),d)
    end
end

function loss3(c, d, r) # normal form 3
    return max(0, euclidean(centerpoint(c)-r, centerpoint(d)) + abs(radius(d)) - abs(radius(c)))
end

function loss4(c, d, r) # normal form 4
    return max(0, euclidean(centerpoint(c) + r, centerpoint(d)) + abs(radius(d)) - abs(radius(c)))
end

classes = Dict()
relations = Dict()
ccount = 1
rcount = 1

function filldict(x)
    global ccount
    if haskey(classes, x)
        return classes[x]
    else
        classes[x] = ccount
        ccount += 1
        return ccount - 1
    end
end

function rfilldict(x)
    global rcount
    if haskey(relations, x)
        return relations[x]
    else
        relations[x] = rcount
        rcount += 1
        return rcount - 1
    end
end

nf1 = Set()
nf2 = Set()
nf3 = Set()
nf4 = Set()

# Reading the file of normalized axioms (generated by Normalize.groovy)
open(ARGS[1]) do file
    for line in eachline(file)
        if (startswith(line, "SubClassOf")) # ignore subproperty axioms
            if occursin(r"ObjectIntersectionOf", line) # normal form 2
                m = match(r"(<.*>).*(<.*>).*(<.*>)", line) # 3 captures: C and D SubClassOf: E
                if m != nothing
                    c, d, e = map(filldict, m.captures)
                    push!(nf2, (c,d,e))
                end
            elseif occursin(r"SubClassOf.<[^\s]*>[\s]*<[^\s]*>.$", line) # normal form 1
                m = match(r"(<.*>).*(<.*>)", line) # 2 captures: C SubClassOf: D
                c,d = map(filldict, m.captures)
                push!(nf1, (c,d))
            elseif occursin(r"SubClassOf.ObjectSomeValuesFrom", line)  # normal form 4
                m = match(r"(<.*>).*(<.*>).*(<.*>)", line) # 3 captures: R some C SubClassOf: D
                r,c,d = m.captures
                c = filldict(c)
                d = filldict(d)
                r = rfilldict(r)
                push!(nf4, (c,d,r))
            elseif occursin(r"SubClassOf.<[^\s]*>[\s]*ObjectSomeValuesFrom", line)  # normal form 3
                m = match(r"(<.*>).*(<.*>).*(<.*>)", line) # 3 captures: C SubClassOf: R some D
                c,r,d = m.captures
                c = filldict(c)
                d = filldict(d)
                r = rfilldict(r)
                push!(nf3, (c,d,r))
            end
        end
    end
end

# there are now rcount-1 relations and ccount-1 classes; generate random embeddings; offload to GPU; make them trainable
#cvec = param(cu(rand(ccount-1, DIM + 1)))
#rvec = param(cu(rand(rcount-1, DIM)))
cvec = param((rand(ccount-1, DIM + 1)))
rvec = param((rand(rcount-1, DIM)))

# setting up reverse maps for faster lookups
rclasses = Dict(value => key for (key, value) in classes)
rrelations = Dict(value => key for (key, value) in relations)

# set to array, just for convenience
nf1arr =  collect(nf1)
nf2arr =  collect(nf2)
nf3arr =  collect(nf3)
nf4arr =  collect(nf4)

# loss functions that work with arrays
function loss1(x::Int, y::Int)
    return loss1(cvec[x,:], cvec[y,:])
end
function loss2(x::Int, y::Int, z::Int)
    return loss2(cvec[x,:], cvec[y,:], cvec[z,:])
end
function loss3(x::Int, y::Int, r::Int)
    return loss3(cvec[x,:], cvec[y,:], rvec[r,:])
end
function loss4(x::Int, y::Int, r::Int)
    return loss4(cvec[x,:], cvec[y,:], rvec[r,:])
end

opt = Momentum([cvec, rvec], LR, decay=0.0)

for k in 1:EPOCHS
    println("Epoch: $k")
    Flux.train!(loss1, nf1, opt)
    Flux.train!(loss2, nf2, opt)
    Flux.train!(loss3, nf3, opt)
    Flux.train!(loss4, nf4, opt)
    if k % 1 == 0
        println("Loss 1: ",  sum([loss1(x,y) for (x,y) in nf1arr]))
        println("Loss 2: ",  sum([loss2(x,y,z) for (x,y,z) in nf2arr]))
        println("Loss 3: ",  sum([loss3(x,y,r) for (x,y,r) in nf3arr]))
        println("Loss 4: ",  sum([loss4(x,y,r) for (x,y,r) in nf4arr]))
    end
end

max_stars_repo_name: leechuck/julia-stuff

max_issues_repo_name: leechuck/julia-stuff

max_forks_repo_name: leechuck/julia-stuff


text:
# type conversion
cint(n) = convert(Cint, n)    # same as int32
csize_t(n) = convert(Csize_t, n)
cchar(n) = convert(Int8,n)
cuchar(n) = convert(Cuchar, n)
cshort(n) = convert(Cshort, n)
cushort(n) = convert(Cushort, n)

# Search directory
searchdir(path,key) = filter(x->contains(x,key), readdir(path))
swapext(f, new_ext) = "$(splitext(f)[1])$new_ext"

# C++ std::vector class thin-wrappings for Julia
cxx"""
template <typename T>
    std::vector<T> stdvector(int size, T x)
      {
        std::vector<T> cppvec(size, x);
        return (cppvec);
      }

template <typename T>
    std::vector<T> stdvectorSzt(std::size_t size, T x)
      {
        std::vector<T> cppvec(size, x);
        return (cppvec);
      }

template <typename T_>
    T_ at(std::vector<T_> cppvec, int index)
       {
           T_ value = cppvec[index];    // does not check for out of bounds (fast but dangerous)
           return(value);
       }

template <typename T_>
    T_ at_(std::vector<T_> cppvec, int index)
       {
          T_ value = cppvec.at(index); // checks for out of bounds (safe but slow)
          return(value);
       }

template <typename T_>
   void stdset(std::vector<T_>& cppvec, int index, T_ value)
       {
          cppvec[index] = value;
       }

template <typename T_>
   void stdset_(std::vector<T_>& cppvec, int index, T_ value)
       {
          cppvec.at(index) = value; // checks for out of bounds (safe but slow)
       }

"""

stdvec(size, value) = @cxxnew stdvector(size, value)
stdvecSzt(size, value) = @cxxnew stdvectorSzt(csize_t(size), value)
stdassign(ccpvec, size, value) = @cxx ccpvec->assign(size,value)
stddata(cppvec) = @cxx cppvec->data()     # Ptr to first elememt
stdempty!(cppvec) = @cxx cppvec->empty()   # check if it is empty
stdcapacity(cppvec) = Int(@cxx cppvec->capacity())
stdpush!(cppvec, value) = @cxx cppvec->push_back(value)
stdpop!(cppvec) = @cxx cppvec->pop_back()
stdsize(cppvec) = Int(@cxx cppvec->size())
stdresize!(cppvec, n::Int) = @cxx cppvec->resize(n)
stdshrink!(cppvec) = @cxx cppvec->shrink_to_fit()
stdswap!(cppvec1, cppvec2) = @cxx cppvec1->swap(cppvec2)

# make sure index is stdsize -1 (C++ indexing starts at 0)
function at(cppvec, index::Int)
   (index < 0 || index > (stdsize(cppvec)-1)) ? throw(ArgumentError("index is out of bounds")) : nothing
   @cxx at(cppvec, index)
end
function at_(cppvec, index::Int)
  # make sure index is stdsize -1 (C++ indexing starts at 0)
  (index < 0 || index > (stdsize(cppvec)-1)) ? throw(ArgumentError("index is out of bounds")) : nothing
   @cxx at_(cppvec, index)   # out of bounds check  (will lead to crash if out of bounds)
end

function set!(cppvec, index, value)
   (index < 0 || index > (stdsize(cppvec)-1)) ? throw(ArgumentError("index is out of bounds")) : nothing
   @cxx stdset(cppvec, index, value)
end

clear(cppvec) = @cxx cppvec->clear()

# Converting julia Array{Int64,1} to an std::vector
function tostdvec(jl_vector::Array{T,1}) where {T}
    # C++ must deduce type from template functions
    vec = stdvec(0,jl_vector[1])

    for i=1:length(jl_vector)
       stdpush!(vec, jl_vector[i])  # index -1 (C++ has 0-indexing)
    end
    return(vec)
end

# C++ string handling
# julia string => std::string
# std::string => julia string

cxx"""
std::string stdstring(const char* word)
      {
        std::string stdword = word;
        return(stdword);
      }
"""

stdstring(word::String) = @cxx stdstring(pointer(word))
juliastring(std_string) = bytestring(@cxx std_string->c_str())
stdstrdata(std_string) = @cxx std_string->data()     # Ptr to first elememt
stdstrempty(std_string) = @cxx std_string->empty()
stdstrsize(std_string) = @cxx std_string->size()
stdstrlength(std_string) = @cxx std_string->length()
stdstrclear(std_string) = @cxx std_string->clear()

#-------------------------------------------------------------------------------------------------------------------#
# Basic core utility functions

# Print Mat uchar array (very rough)
cxx"""
void printMat(const cv::Mat& mat)
{
    for(int i = 0; i < mat.rows; i++)
    {
        for(int j=0; j < mat.cols; j++)
        {
            std::cout << (int) mat.at<uchar>(i,j) << " "; break;
        }
    }
}
"""

printMat(img) = @cxx printMat(img)

#-------------------------------------------------------------------------------------------------------------------#
# Image processing (imgproc)

# imreplace - uses a rectangular ROI to copy a region from one image to another
cxx"""
cv::Mat imreplace(cv::Mat& src,cv::Mat& dst, cv::Rect roi)
{
    //initialize the ROI (location to copy)
    cv::Mat dstROI(dst(roi));
    //perform the copy
    src.copyTo(dstROI);
    return(dst);
}
"""
imreplace(src, dst, roi) = @cxx imreplace(src, dst, roi)

# Support for image convolution
cxx"""
// Function getSum returns total sum of all the elements of given matrix.

float getKernelSum(cv::Mat kernel)
{
   float sum = 0;
   for(typeof(kernel.begin<float>()) it = kernel.begin<float>(); it!= kernel.end<float>() ; it++)
   {
      sum+=*it;
   }
   return sum;
}

// Normalize the mask (kernel)
cv::Mat normalizeKernel (cv::Mat kernel, double ksum) {
    kernel = kernel/ksum;
    return(kernel);
    }
"""

# kernel = ones(5,5,CV_32F)
getKernelSum(kernel) = @cxx getKernelSum(kernel)
# ksum = getKernelSum(kernel)
normalizeKernel(kernel, ksum) = @cxx normalizeKernel(kernel, ksum)

# TO_DO:  Similar utility functions are likely required for other functions in imgproc

#-------------------------------------------------------------------------------------------------------------------#
# GUI interface (highui)

# Display an image in a window and close upon key press (and delay)
# For concurrent multiple display, set multi = "ON"
function imdisplay(img, windowName::String, flag=WINDOW_AUTOSIZE)
    namedWindow(windowName, flag)
    imshow(windowName, img)
end

function closeWindows(delay, key, windowName)
      (waitkey(delay) == key && windowName != "") ? destroyWindow(windowName) : destroyAllWindows()
end
# closeWindows(0,27,"")

function im2tile(imgArray, windowName::String, flag=WINDOW_AUTOSIZE)
    canvas = Mat()

    for i=1:length(imgArray)
        # check that images have same dims, format and channels
        if (i > 1)
            (cvtypeval(imgArray[i]) != cvtypeval(imgArray[i-1]) ||
             rows(imgArray[i]) != rows(imgArray[i-1]) ||
             cols(imgArray[i]) != cols(imgArray[i-1])) ?
          throw(ArgumentError("Images must have same dimensions and format")) : nothing
        end
        push(canvas, imgArray[i])
    end

    imdisplay(canvas, windowName, flag)
end

#-------------------------------------------------------------------------------------------------------------------#
# Video capture (video)

function videocam(device = CAP_ANY)
    cam = videoCapture(device)    # open Video device
    !isOpened(cam) ? throw(ArgumentError("Can not open camera!")) : nothing
    namedWindow("Video")
    frame = Mat()

    # Loop until user presses ESC or frame is empty
    while(true)
        if !(videoRead(cam, frame))
            throw(ArgumentError("Can not acquire video!"))
            break
        end

        imshow("Video", frame)

        if (waitkey(30) == 27)
            destroyAllWindows()
            break
       end
   end
  release(cam)
end

# Webstreaming
# src: full http link to the video source
# e.g., NASA TV: vid = "http://www.nasa.gov/multimedia/nasatv/NTV-Public-IPS.m3u8"
function webstream(src::String)
    cam = videoCapture(src)    # open Video device
    !isOpened(cam) ? throw(ArgumentError("Can not open stream!")) : nothing
    namedWindow("Web stream")
    frame = Mat()

    # Loop until user presses ESC or frame is empty
    while(true)
        if !(videoRead(cam, frame))
            throw(ArgumentError("Can not acquire video stream!"))
            break
        end

        imshow("Web stream", frame)

        if (waitkey(30) == 27)
            destroyAllWindows()
            break
       end
   end
  release(cam)
end

function videoWrite(cam, filename::String, fps::Float64, nframes = 0, frameSize=cvSize(0,0), fourcc=-1, isColor=true)
    !(isOpened(cam)) ? throw(ArgumentError("Video source is not available!")) : nothing

    # Set the video capture frame size based on camera input WIDTH and HEIGHT
    width = getVideoId(cam, CAP_PROP_FRAME_WIDTH)
    height = getVideoId(cam, CAP_PROP_FRAME_HEIGHT)

    frame = Mat(Int(height), Int(width), CV_8UC3)

    # Set the frameSize of output video frames
    (frameSize.data[1] == 0) ? frameSize = cvSize(Int(width), Int(height)) :
         throw(ArgumentError("Output frame dimension is wrong"))

    # Initialize and open video writer
    writer = videoWriter(filename, fourcc, fps, frameSize, isColor)
    openWriter(writer, filename, fourcc, fps, frameSize, isColor)
    !(isOpened(writer)) ? throw(ArgumentError("Can not open the video writer!")) : nothing

    # create a window for display
    namedWindow("Welcome!")

    # Loop until user presses ESC or frame is empty
    count = 0
    while(true)
        if !(videoRead(cam, frame))
            throw(ArgumentError("Can not acquire video!"))
            break
        end

        writeVideo(writer, frame)
        imshow("Welcome!", frame)

        (nframes > 0) ? count += 1 : nothing
        if ((waitkey(30) == 27) || (count == nframes))
            destroyAllWindows()
            break
       end
   end
   release(cam)
end


text:
"""
3D simulation of spacecraft

sphere domain with R = 5 and center = (2, 0, 0)

structured mesh: 65×37×41
- i: index along positive x direction
- j: index around the annulus perpendicular to x
- k: index from spacecraft to the outer surface

"""

using KitBase, CSV, DataFrames

cd(@__DIR__)

set = (ni = 65, nj = 37, nk = 41, xc = [2.0, 0.0, 0.0])

function point_coords(df)
    ni, nj, nk = set.ni, set.nj, set.nk
    x = reshape(collect(df.x), ni, nj, nk)
    y = reshape(collect(df.y), ni, nj, nk)
    z = reshape(collect(df.z), ni, nj, nk)
    points = cat(x, y, z, dims = 4)

    return points
end

df = CSV.read("../assets/mesh/spacecraft.csv", DataFrame)
points = point_coords(df)

#--- geometric illustration ---#
[
    @show point_distance(points[i, j, end, :], set.xc) for i in axes(points, 1),
    j in axes(points, 2)
]
[@show points[i, 1, end, :] for i in axes(points, 1)]
[@show points[30, j, end, :] for j in axes(points, 2)]

#--- physical space ---#
function physical_space(points)
    ni, nj, nk = size(points)[1:3] .- 1
    x0, y0, z0 = minimum(points[:, :, :, 1]),
    minimum(points[:, :, :, 2]),
    minimum(points[:, :, :, 3])
    x1, y1, z1 = maximum(points[:, :, :, 1]),
    maximum(points[:, :, :, 2]),
    maximum(points[:, :, :, 3])
end

physical_space(points)

ni, nj, nk = size(points)[1:3] .- 1
x0, y0, z0 =
    minimum(points[:, :, :, 1]), minimum(points[:, :, :, 2]), minimum(points[:, :, :, 3])
x1, y1, z1 =
    maximum(points[:, :, :, 1]), maximum(points[:, :, :, 2]), maximum(points[:, :, :, 3])
x = zeros(ni, nj, nk)
y = zero(x)
z = zero(x)
dx = zero(x)
dy = zero(x)
dz = zero(x)
for i = 1:ni, j = 1:nj, k = 1:nk
    x[i, j, k] = points[i, j, k, 1]
end


text:
## HSP family of heuristics ##
export HSPHeuristic, HAdd, HMax
export HSPRHeuristic, HAddR, HMaxR

"Precomputed domain information for HSP heuristic."
struct HSPCache
    domain::Domain # Preprocessed domain
    axioms::Vector{Clause} # Preprocessed axioms
    preconds::Dict{Symbol,Vector{Vector{Term}}} # Preconditions in DNF
    additions::Dict{Symbol,Vector{Term}} # Action add lists
end

"HSP family of delete-relaxation heuristics."
mutable struct HSPHeuristic <: Heuristic
    op::Function # Aggregator (e.g. maximum, sum) for fact costs
    cache::HSPCache # Precomputed domain information
    pre_key::UInt64 # Key to check if information needs to be precomputed again
    HSPHeuristic(op) = new(op)
end

Base.hash(heuristic::HSPHeuristic, h::UInt) =
    hash(heuristic.op, hash(HSPHeuristic, h))

function precompute!(h::HSPHeuristic,
                     domain::Domain, state::State, spec::Specification)
    # Check if cache has already been computed
    if is_precomputed(h, domain, state, spec) return h end
    h.pre_key = objectid(domain) # Precomputed data is unique to each domain
    domain = copy(domain) # Make a local copy of the domain
    # Preprocess axioms
    axioms = regularize_clauses(domain.axioms) # Regularize domain axioms
    axioms = [Clause(ax.head, [t for t in ax.body if t.name != :not])
              for ax in axioms] # Remove negative literals
    domain.axioms = Clause[] # Remove axioms so they do not affect execution
    # Preprocess actions
    preconds = Dict{Symbol,Vector{Vector{Term}}}()
    additions = Dict{Symbol,Vector{Term}}()
    for (act_name, act_def) in domain.actions
        # Convert preconditions to DNF without negated literals
        conds = get_preconditions(act_def; converter=to_dnf)
        conds = [c.args for c in conds.args]
        for c in conds filter!(t -> t.name != :not, c) end
        preconds[act_name] = conds
        # Extract additions from each effect
        diff = effect_diff(act_def.effect)
        additions[act_name] = diff.add
    end
    h.cache = HSPCache(domain, axioms, preconds, additions)
    return h
end

function is_precomputed(h::HSPHeuristic,
                        domain::Domain, state::State, spec::Specification)
    return (isdefined(h, :cache) && objectid(domain) == h.pre_key)
end

function compute(h::HSPHeuristic,
                 domain::Domain, state::State, spec::Specification)
    # Precompute if necessary
    if !is_precomputed(h, domain, state, spec)
        precompute!(h, domain, state, spec) end
    @unpack op, cache = h
    @unpack domain = cache
    @unpack types, facts = state
    goals = get_goal_terms(spec)
    # Initialize fact costs in a GraphPlan-style graph
    fact_costs = Dict{Term,Float64}(f => 0 for f in facts)
    while true
        facts = Set(keys(fact_costs))
        state = State(types, facts, Dict{Symbol,Any}())
        if is_goal(spec, domain, state)
            return op([0; [fact_costs[g] for g in goals]]) end
        # Compute costs of one-step derivations of domain axioms
        for ax in cache.axioms
            _, subst = resolve(ax.body, [Clause(f, []) for f in facts])
            for s in subst
                body = [substitute(t, s) for t in ax.body]
                cost = op([0; [get(fact_costs, f, 0) for f in body]])
                derived = substitute(ax.head, s)
                if cost < get(fact_costs, derived, Inf)
                    fact_costs[derived] = cost end
            end
        end
        # Compute costs of all effects of available actions
        actions = available(state, domain)
        for act in actions
            act_args = domain.actions[act.name].args
            subst = Subst(var => val for (var, val) in zip(act_args, act.args))
            # Look-up preconds and substitute vars
            preconds = cache.preconds[act.name]
            preconds = [[substitute(t, subst) for t in c] for c in preconds]
            # Compute cost of reaching each action
            cost = minimum([[op([0; [get(fact_costs, f, 0) for f in conj]])
                             for conj in preconds]; Inf])
            # Compute cost of reaching each added fact
            additions = [substitute(a, subst) for a in cache.additions[act.name]]
            cost = cost + 1 # TODO: Handle arbitrary action costs
            for fact in additions
                if cost < get(fact_costs, fact, Inf)
                    fact_costs[fact] = cost end
            end
        end
        # Terminate if there's no change to the number of facts
        if length(fact_costs) == length(facts) && keys(fact_costs) == facts
            return Inf end
    end
end

"HSP heuristic where a fact's cost is the maximum cost of its dependencies."
HMax(args...) = HSPHeuristic(maximum, args...)

"HSP heuristic where a fact's cost is the summed cost of its dependencies."
HAdd(args...) = HSPHeuristic(sum, args...)

"HSPr family of delete-relaxation heuristics for regression search."
mutable struct HSPRHeuristic <: Heuristic
    op::Function
    fact_costs::Dict{Term,Float64} # Est. cost of reaching each fact from goal
    pre_key::NTuple{3,UInt64} # Key to check if precomputation needs to re-run
    HSPRHeuristic(op) = new(op)
end

Base.hash(heuristic::HSPRHeuristic, h::UInt) =
    hash(heuristic.op, hash(HSPRHeuristic, h))

function precompute!(h::HSPRHeuristic,
                     domain::Domain, state::State, spec::Specification)
    # Check if data has already been computed
    if is_precomputed(h, domain, state, spec) return h end
    # Precomputed data is tied to all three inputs
    h.pre_key = (objectid(domain), hash(state), hash(spec))
    @unpack op = h
    @unpack types, facts = state
    goals = get_goal_terms(spec)
    # Preprocess domain and axioms
    domain = copy(domain)
    axioms = regularize_clauses(domain.axioms)
    axioms = [Clause(ax.head, [t for t in ax.body if t.name != :not])
              for ax in axioms]
    domain.axioms = Clause[]
    # Preprocess actions
    preconds = Dict{Symbol,Vector{Vector{Term}}}()
    additions = Dict{Symbol,Vector{Term}}()
    for (act_name, act_def) in domain.actions
        # Convert preconditions to DNF without negated literals
        conds = get_preconditions(act_def; converter=to_dnf)
        conds = [c.args for c in conds.args]
        for c in conds filter!(t -> t.name != :not, c) end
        preconds[act_name] = conds
        # Extract additions from each effect
        diff = effect_diff(act_def.effect)
        additions[act_name] = diff.add
    end
    # Initialize fact costs in a GraphPlan-style graph
    fact_costs = Dict{Term,Float64}(f => 0 for f in PDDL.get_facts(state))
    while true
        facts = Set(keys(fact_costs))
        state = State(types, facts, Dict{Symbol,Any}())
        # Compute costs of one-step derivations of domain axioms
        for ax in axioms
            _, subst = resolve(ax.body, [Clause(f, []) for f in facts])
            for s in subst
                body = [substitute(t, s) for t in ax.body]
                cost = op([0; [get(fact_costs, f, 0) for f in body]])
                derived = substitute(ax.head, s)
                if cost < get(fact_costs, derived, Inf)
                    fact_costs[derived] = cost end
            end
        end
        # Compute costs of all effects of available actions
        actions = available(state, domain)
        for act in actions
            act_args = domain.actions[act.name].args
            subst = Subst(var => val for (var, val) in zip(act_args, act.args))
            # Look-up preconds and substitute vars
            conds = preconds[act.name]
            conds = [[substitute(t, subst) for t in c] for c in conds]
            # Compute cost of reaching each action
            cost = minimum([[op([0; [get(fact_costs, f, 0) for f in conj]])
                             for conj in conds]; Inf])
            # Compute cost of reaching each added fact
            added = [substitute(a, subst) for a in additions[act.name]]
            cost = cost + 1 # TODO: Handle arbitrary action costs
            for fact in added
                if cost < get(fact_costs, fact, Inf)
                    fact_costs[fact] = cost end
            end
        end
        # Terminate when there's no change to the number of facts
        if length(fact_costs) == length(facts) && keys(fact_costs) == facts
            break end
    end
    h.fact_costs = fact_costs
    return h
end

function is_precomputed(h::HSPRHeuristic,
                        domain::Domain, state::State, spec::Specification)
    return (isdefined(h, :fact_costs) && objectid(domain) == h.pre_key[1]
            && hash(state) == h.pre_key[2] && hash(spec) == h.pre_key[3])
end

function compute(h::HSPRHeuristic,
                 domain::Domain, state::State, spec::Specification)
    # Compute cost of achieving all facts in current state
    return h.op([0; [get(h.fact_costs, f, 0) for f in PDDL.get_facts(state)]])
end

"HSPr heuristic where a fact's cost is the maximum cost of its dependencies."
HMaxR(args...) = HSPRHeuristic(maximum, args...)

"HSPr heuristic where a fact's cost is the summed cost of its dependencies."
HAddR(args...) = HSPRHeuristic(sum, args...)


text:
@testset "LAS" begin
    workdir = dirname(dirname(pathof(FileManager)))
    file = joinpath(workdir,"test","Files/PC.las")
    @show file
    PC = FileManager.PointCloud(rand(3,100),LasIO.FixedPointNumbers.N0f16.(rand(3,100)))
    FileManager.save_pointcloud(file, PC,  "TEST")
    PC_load = FileManager.source2pc(file)
    @test PC_load.n_points == PC.n_points
    @test PC_load.dimension == PC.dimension
end


text:
function divrem(a,b)
    q = []
    r = copy(a)
    while length(r) >= length(b)
        k = length(r) - length(b)
        c = vcat(b,zeros(Int64,k))
        push!(q,r[1])
        r = r-r[1]*c
        deleteat!(r,1)
    end
    return q,r
end


text:
module WinKnownPaths

using UUIDs

# https://docs.microsoft.com/en-us/windows/win32/api/shlobj_core/nf-shlobj_core-shgetknownfolderpath
# https://docs.microsoft.com/en-us/windows/win32/shell/knownfolderid

module FOLDERID
using UUIDs

const AccountPictures         = UUID("008ca0b1-55b4-4c56-b8a8-4de4b299d3be")
const AdminTools              = UUID("724EF170-A42D-4FEF-9F26-B60E846FBA4F")
const ApplicationShortcuts    = UUID("A3918781-E5F2-4890-B3D9-A7E54332328C")
const CameraRoll              = UUID("AB5FB87B-7CE2-4F83-915D-550846C9537B")
const CDBurning               = UUID("9E52AB10-F80D-49DF-ACB8-4330F5687855")
const CommonAdminTools        = UUID("D0384E7D-BAC3-4797-8F14-CBA229B392B5")
const CommonOEMLinks          = UUID("C1BAE2D0-10DF-4334-BEDD-7AA20B227A9D")
const CommonPrograms          = UUID("0139D44E-6AFE-49F2-8690-3DAFCAE6FFB8")
const CommonStartMenu         = UUID("A4115719-D62E-491D-AA7C-E74B8BE3B067")
const CommonStartup           = UUID("82A5EA35-D9CD-47C5-9629-E15D2F714E6E")
const CommonTemplates         = UUID("B94237E7-57AC-4347-9151-B08C6C32D1F7")
const Contacts                = UUID("56784854-C6CB-462b-8169-88E350ACB882")
const Cookies                 = UUID("2B0F765D-C0E9-4171-908E-08A611B84FF6")
const Desktop                 = UUID("B4BFCC3A-DB2C-424C-B029-7FE99A87C641")
const DeviceMetadataStore     = UUID("5CE4A5E9-E4EB-479D-B89F-130C02886155")
const Documents               = UUID("FDD39AD0-238F-46AF-ADB4-6C85480369C7")
const DocumentsLibrary        = UUID("7B0DB17D-9CD2-4A93-9733-46CC89022E7C")
const Downloads               = UUID("374DE290-123F-4565-9164-39C4925E467B")
const Favorites               = UUID("1777F761-68AD-4D8A-87BD-30B759FA33DD")
const Fonts                   = UUID("FD228CB7-AE11-4AE3-864C-16F3910AB8FE")
const GameTasks               = UUID("054FAE61-4DD8-4787-80B6-090220C4B700")
const History                 = UUID("D9DC8A3B-B784-432E-A781-5A1130A75963")
const ImplicitAppShortcuts    = UUID("BCB5256F-79F6-4CEE-B725-DC34E402FD46")
const InternetCache           = UUID("352481E8-33BE-4251-BA85-6007CAEDCF9D")
const Libraries               = UUID("1B3EA5DC-B587-4786-B4EF-BD1DC332AEAE")
const Links                   = UUID("bfb9d5e0-c6a9-404c-b2b2-ae6db6af4968")
const LocalAppData            = UUID("F1B32785-6FBA-4FCF-9D55-7B8E7F157091")
const LocalAppDataLow         = UUID("A520A1A4-1780-4FF6-BD18-167343C5AF16")
const LocalizedResourcesDir   = UUID("2A00375E-224C-49DE-B8D1-440DF7EF3DDC")
const Music                   = UUID("4BD8D571-6D19-48D3-BE97-422220080E43")
const MusicLibrary            = UUID("2112AB0A-C86A-4FFE-A368-0DE96E47012E")
const NetHood                 = UUID("C5ABBF53-E17F-4121-8900-86626FC2C973")
const OriginalImages          = UUID("2C36C0AA-5812-4b87-BFD0-4CD0DFB19B39")
const PhotoAlbums             = UUID("69D2CF90-FC33-4FB7-9A0C-EBB0F0FCB43C")
const PicturesLibrary         = UUID("A990AE9F-A03B-4E80-94BC-9912D7504104")
const Pictures                = UUID("33E28130-4E1E-4676-835A-98395C3BC3BB")
const Playlists               = UUID("DE92C1C7-837F-4F69-A3BB-86E631204A23")
const PrintHood               = UUID("9274BD8D-CFD1-41C3-B35E-B13F55A758F4")
const Profile                 = UUID("5E6C858F-0E22-4760-9AFE-EA3317B67173")
const ProgramData             = UUID("62AB5D82-FDC1-4DC3-A9DD-070D1D495D97")
const ProgramFiles            = UUID("905e63b6-c1bf-494e-b29c-65b732d3d21a")
const ProgramFilesX64         = UUID("6D809377-6AF0-444b-8957-A3773F02200E")
const ProgramFilesX86         = UUID("7C5A40EF-A0FB-4BFC-874A-C0F2E0B9FA8E")
const ProgramFilesCommon      = UUID("F7F1ED05-9F6D-47A2-AAAE-29D317C6F066")
const ProgramFilesCommonX64   = UUID("6365D5A7-0F0D-45E5-87F6-0DA56B6A4F7D")
const ProgramFilesCommonX86   = UUID("DE974D24-D9C6-4D3E-BF91-F4455120B917")
const Programs                = UUID("A77F5D77-2E2B-44C3-A6A2-ABA601054A51")
const Public                  = UUID("DFDF76A2-C82A-4D63-906A-5644AC457385")
const PublicDesktop           = UUID("C4AA340D-F20F-4863-AFEF-F87EF2E6BA25")
const PublicDocuments         = UUID("ED4824AF-DCE4-45A8-81E2-FC7965083634")
const PublicDownloads         = UUID("3D644C9B-1FB8-4f30-9B45-F670235F79C0")
const PublicGameTasks         = UUID("DEBF2536-E1A8-4c59-B6A2-414586476AEA")
const PublicLibraries         = UUID("48DAF80B-E6CF-4F4E-B800-0E69D84EE384")
const PublicMusic             = UUID("3214FAB5-9757-4298-BB61-92A9DEAA44FF")
const PublicPictures          = UUID("B6EBFB86-6907-413C-9AF7-4FC2ABF07CC5")
const PublicRingtones         = UUID("E555AB60-153B-4D17-9F04-A5FE99FC15EC")
const PublicUserTiles         = UUID("0482af6c-08f1-4c34-8c90-e17ec98b1e17")
const PublicVideos            = UUID("2400183A-6185-49FB-A2D8-4A392A602BA3")
const QuickLaunch             = UUID("52a4f021-7b75-48a9-9f6b-4b87a210bc8f")
const Recent                  = UUID("AE50C081-EBD2-438A-8655-8A092E34987A")
const RecordedTVLibrary       = UUID("1A6FDBA2-F42D-4358-A798-B74D745926C5")
const ResourceDir             = UUID("8AD10C31-2ADB-4296-A8F7-E4701232C972")
const Ringtones               = UUID("C870044B-F49E-4126-A9C3-B52A1FF411E8")
const RoamingAppData          = UUID("3EB685DB-65F9-4CF6-A03A-E3EF65729F3D")
const RoamedTileImages        = UUID("AAA8D5A5-F1D6-4259-BAA8-78E7EF60835E")
const RoamingTiles            = UUID("00BCFC5A-ED94-4e48-96A1-3F6217F21990")
const SampleMusic             = UUID("B250C668-F57D-4EE1-A63C-290EE7D1AA1F")
const SamplePictures          = UUID("C4900540-2379-4C75-844B-64E6FAF8716B")
const SamplePlaylists         = UUID("15CA69B3-30EE-49C1-ACE1-6B5EC372AFB5")
const SampleVideos            = UUID("859EAD94-2E85-48AD-A71A-0969CB56A6CD")
const SavedGames              = UUID("4C5C32FF-BB9D-43b0-B5B4-2D72E54EAAA4")
const SavedSearches           = UUID("7d1d3a04-debb-4115-95cf-2f29da2920da")
const Screenshots             = UUID("b7bede81-df94-4682-a7d8-57a52620b86f")
const SearchHistory           = UUID("0D4C3DB6-03A3-462F-A0E6-08924C41B5D4")
const SearchTemplates         = UUID("7E636BFE-DFA9-4D5E-B456-D7B39851D8A9")
const SendTo                  = UUID("8983036C-27C0-404B-8F08-102D10DCFD74")
const SidebarDefaultParts     = UUID("7B396E54-9EC5-4300-BE0A-2482EBAE1A26")
const SidebarParts            = UUID("A75D362E-50FC-4fb7-AC2C-A8BEAA314493")
const SkyDrive                = UUID("A52BBA46-E9E1-435f-B3D9-28DAA648C0F6")
const SkyDriveCameraRoll      = UUID("767E6811-49CB-4273-87C2-20F355E1085B")
const SkyDriveDocuments       = UUID("24D89E24-2F19-4534-9DDE-6A6671FBB8FE")
const SkyDrivePictures        = UUID("339719B5-8C47-4894-94C2-D8F77ADD44A6")
const StartMenu               = UUID("625B53C3-AB48-4EC1-BA1F-A1EF4146FC19")
const Startup                 = UUID("B97D20BB-F46A-4C97-BA10-5E3608430854")
const System                  = UUID("1AC14E77-02E7-4E5D-B744-2EB1AE5198B7")
const SystemX86               = UUID("D65231B0-B2F1-4857-A4CE-A8E7C6EA7D27")
const Templates               = UUID("A63293E8-664E-48DB-A079-DF759E0509F7")
const UserPinned              = UUID("9E3995AB-1F9C-4F13-B827-48B24B6C7174")
const UserProfiles            = UUID("0762D272-C50A-4BB0-A382-697DCD729B80")
const UserProgramFiles        = UUID("5CD7AEE2-2219-4A67-B85D-6C9CE15660CB")
const UserProgramFilesCommon  = UUID("BCBD3057-CA5C-4622-B42D-BC56DB0AE516")
const Videos                  = UUID("18989B1D-99B5-455B-841C-AB7C74E4DDFC")
const VideosLibrary           = UUID("491E922F-5643-4AF4-A7EB-4E7A138D8174")
const Windows                 = UUID("F38BF404-1D43-42F2-9305-67DE0B28FC23")
end

import .FOLDERID

# KNOWNFOLDERID is an alias for a GUID
struct GUID
    l1::UInt32
    w1::UInt16
    w2::UInt16
    b::NTuple{8,UInt8}
end
function Base.convert(::Type{GUID}, uuid::UUID)
    bytes = uuid.value
    data1 = UInt32((bytes >> 96) & 0xffffffff)
    data2 = UInt16((bytes >> 80) & 0xffff)
    data3 = UInt16((bytes >> 64) & 0xffff)
    data4 = ntuple(i -> UInt8((bytes >> (64 - i*8)) & 0xff), 8)
    return GUID(data1, data2, data3, data4)
end

const KF_FLAG_DEFAULT = 0x00000000

"""
    WinKnownPaths.path(folderid::UUID) -> AbstractString

    Retrieve the full path of a known Windows path identified by the folder's id
    specified in the module `FOLDERID`.

# Examples

```julia
julia> import WinKnownPaths: FOLDERID

julia> WinKnownPaths.path(FOLDERID.System)
"C:\\WINDOWS\\system32"

julia> WinKnownPaths.path(FOLDERID.Fonts)
"C:\\WINDOWS\\Fonts"
```

"""
function path(folderid::UUID)
    pathptr = Ref{Ptr{Cwchar_t}}()
    result = ccall((:SHGetKnownFolderPath,:shell32), stdcall, Cint,
                    (Ref{GUID}, UInt32, Ptr{Cvoid}, Ref{Ptr{Cwchar_t}}), folderid, KF_FLAG_DEFAULT, C_NULL, pathptr)
    Base.systemerror("SHGetKnownFolderPath", result != 0)
    pathbuf = unsafe_wrap(Vector{Cwchar_t}, pathptr[], ccall(:wcslen, Csize_t, (Ptr{Cwchar_t},), pathptr[]))
    path = transcode(String, pathbuf)
    ccall((:CoTaskMemFree, :ole32), Cvoid, (Ptr{Cvoid},), pathptr[])
    return path
end

end


text:
# N-Queens Solution with Genetic Algorithms in Julia
# Victhor S. Sartório

import Distributions
import StatsBase

# Generates a population containing $nchroms members of size $nqueens
# with a discrete uniform distribution with support in [1, $nqueens]
# in the form of a matrix, with each row being a member.
function genpop(nchroms::Int, nqueens::Int)
	# Check parameter validity
	if nchroms <= 9 || nqueens <= 3
		throw(DomainError())
	end

	# Set up the distribution
	U = Distributions.DiscreteUniform(1, nqueens)

	# Generate and return the matrix
	return reshape(rand(U, nchroms*nqueens), nchroms, nqueens)
end

# Calculates the fitness of each member of a population
function fitness(population::Array{Int,2})
	# Get the dimensions
	(nchroms, nqueens) = size(population)
	# Allocate return vector
	fitvec = zeros(Float64, nchroms)

	# For each member calculate fitness
	for k in 1:nchroms
		# afit measures the number of attacks, which is the opposite of fitness
		local afit::UInt = 0
		for i in 1:nqueens
			pivot = population[k,i]
			for j in i+1:nqueens
				current = population[k,j]
				if pivot == current
					afit += 1
				end
				if pivot+i-j == current || pivot-i+j == current
					afit += 1
				end
			end
		end
		# We invert it, so the less "antifitness" the higher the fitness, up to Inf
		fitvec[k] = 1/Float64(afit)
	end

	return fitvec
end

# Crossover
function crossover(population::Array{Int,2}, fitvec::Array{Float64,1}, mutrate::Float64)
	# Get the dimensions
	(nchroms, nqueens) = size(population)
	# Allocate the matrix for the new population
	newpop = zeros(Int, (nchroms, nqueens))

	# Extract random numbers we deterministically know will be needed
	# to save on function call overhead and possibly take advantage of
	# vectorization in the function implementations.
	sampler = StatsBase.sample(1:nchroms, StatsBase.WeightVec(fitvec), 2*nchroms)

	crosspt = rand(1:nqueens, nchroms)
	mutate = rand(MersenneTwister(), nchroms)

	# Generate new population
	for k in 1:nchroms
		dad = population[sampler[2*k-1],:]
		mom = population[sampler[2*k],:]
		pivot = crosspt[k]
		kid = [dad[1:pivot]; mom[pivot+1:nqueens]]

		if mutate[k] < mutrate
			(at, val) = rand(1:nqueens, 2)
			kid[at] = val
		end

		newpop[k,:] = kid
	end

	return newpop
end

# Run
function solvenq(n::Int, popsize::Int, mutrate::Float64)
	pop = genpop(popsize, n)
	fit = fitness(pop)

	while !(Inf in fit)
		pop = crossover(pop, fit, mutrate)
		fit = fitness(pop)
	end

	for k in 1:popsize
		if fit[k] == Inf
			return pop[k,:]
		end
	end
end

max_stars_repo_name: dimassores/GeneticQueens

max_issues_repo_name: dimassores/GeneticQueens

max_forks_repo_name: dimassores/GeneticQueens


text:
using CryptoGroups: Group, Hash, PRG, RO, PGroup, ROPRG
using XMLDict: parse_xml

using Base: @kwdef

@kwdef struct ProtocolSpec{G<:Group} <: Verifier
    g::G
    nr::Int32 = Int32(100)
    nv::Int32 = Int32(256)
    ne::Int32 = Int32(256)
    prghash::Hash = Hash("sha256")
    rohash::Hash = Hash("sha256")
    version::String = "3.0.4"
    sid::String = "SessionID"
    auxsid::String = "default"
end

function ProtocolSpec(protinfo::AbstractDict; auxsid = "default")

    s_H = protinfo["rohash"]
    s_PRG = protinfo["prg"]
    s_Gq = protinfo["pgroup"]

    prghash = Hash(map_hash_name(protinfo["prg"]))
    rohash = Hash(map_hash_name(protinfo["rohash"]))

    nr = parse(Int32, protinfo["statdist"])
    nv = parse(Int32, protinfo["vbitlenro"])
    ne = parse(Int32, protinfo["ebitlenro"])

    g = unmarshal(BigInt, decode(split(s_Gq, "::")[2]))

    version = String(protinfo["version"])
    sid = String(protinfo["sid"])

    return ProtocolSpec(; g, nr, nv, ne, prghash, rohash, version, sid, auxsid)
end

function ProtocolSpec(PROT_INFO::AbstractString; auxsid = "default")

    xml = String(read(PROT_INFO))
    protinfo = parse_xml(xml)

    return ProtocolSpec(protinfo; auxsid)
end

function marshal_s_Gq(g::PGroup)

    M = bitlength(order(g))

    tree = marshal(g)
    str = "ModPGroup(safe-prime modulus=2*order+1. order bit-length = $M)::" * string(tree)

    return Leaf(str)
end

function ro_prefix(spec::ProtocolSpec)

    (; version, sid, auxsid, rohash, prghash, g, nr, nv, ne) = spec

    s_PRG = map_hash_name_back(prghash)
    s_H = map_hash_name_back(rohash)

    s_Gq = marshal_s_Gq(g)

    data = (version, sid * "." * auxsid, nr, nv, ne, s_PRG, s_Gq, s_H)

    tree = Tree(data)
    binary = encode(tree)

    ρ = rohash(binary)

    return ρ
end

struct VShuffleProof{G<:Group} <: Proof
    μ::Vector{G}
    τ::Tuple{Vector{G}, G, Vector{G}, G, G, Tuple{G, G}}
    σ::Tuple{BigInt, Vector{BigInt}, BigInt, BigInt, Vector{BigInt}, BigInt}
end

==(x::VShuffleProof{G}, y::VShuffleProof{G}) where G <: Group = x.μ == y.μ && x.τ == y.τ && x.σ == y.σ

function VShuffleProof(proof::PoSProof)

    (; 𝐜, 𝐜̂, t, s) = proof

    G = typeof(𝐜[1])

    𝐮 = 𝐜
    μ = 𝐮

    (t₁, t₂, t₃, t₄, 𝐭̂) = t
    𝐁 = 𝐜̂
    𝐁′= 𝐭̂
    D′ = t₂
    A′ = t₃
    C′ = t₁
    F′ = t₄

    τ = (𝐁, A′, 𝐁′, C′, D′, F′)

    (s₁, s₂, s₃, s₄, 𝐬̂, 𝐬′) = s
    𝐤_B = 𝐬̂
    𝐤_E = 𝐬′
    k_D = s₂
    k_A = s₃
    k_C = s₁
    k_F = s₄

    σ = (k_A, 𝐤_B, k_C, k_D, 𝐤_E, k_F)

    vproof = VShuffleProof(μ, τ, σ)

    return vproof
end

function PoSProof(vproof::VShuffleProof)

    (; μ, τ, σ) = vproof

    𝐮 = μ
    𝐜 = 𝐮

    (𝐁, A′, 𝐁′, C′, D′, F′) = τ

    𝐜̂ = 𝐁
    𝐭̂ = 𝐁′
    t₂ = D′
    t₃ = A′
    t₁ = C′
    t₄ = F′

    t = (t₁, t₂, t₃, t₄, 𝐭̂)

    (k_A, 𝐤_B, k_C, k_D, 𝐤_E, k_F) = σ

    𝐬̂ = 𝐤_B
    𝐬′ = 𝐤_E
    s₂ = k_D
    s₃ = k_A
    s₁ = k_C
    s₄ = k_F

    s = (s₁, s₂, s₃, s₄, 𝐬̂, 𝐬′)

    proof = PoSProof(𝐜, 𝐜̂, t, s)

    return proof
end

function load_verificatum_proposition(basedir::AbstractString, auxsid::AbstractString)

    PUBLIC_KEY = "$basedir/publicKey"

    tree = decode(read(PUBLIC_KEY))
    pk, g = unmarshal_publickey(tree)

    NIZKP = basedir * "/dir/nizkp/$auxsid/"

    CIPHERTEXTS = "$NIZKP/Ciphertexts.bt"
    SHUFFLED_CIPHERTEXTS = "$NIZKP/ShuffledCiphertexts.bt"

    G = typeof(g)

    L_tree = decode(read(CIPHERTEXTS))
    L′_tree = decode(read(SHUFFLED_CIPHERTEXTS))

    𝔀 = convert(ElGamal{G}, L_tree) ## Is there anything I can do so that I would get a concrete type here?
    𝔀′ = convert(ElGamal{G}, L′_tree)

    return Shuffle(g, pk, 𝔀, 𝔀′)
end

function load_verificatum_proof(proofs::AbstractString, g::Group)

    PERMUTATION_COMMITMENT = "$proofs/PermutationCommitment01.bt"
    PoS_COMMITMENT = "$proofs/PoSCommitment01.bt"
    PoS_REPLY = "$proofs/PoSReply01.bt"

    G = typeof(g)

    μ_tree = decode(read(PERMUTATION_COMMITMENT))
    μ = convert(Vector{G}, μ_tree)

    τ_tree = decode(read(PoS_COMMITMENT))
    τ = convert(Tuple{Vector{G}, G, Vector{G}, G, G, Tuple{G, G}}, τ_tree)

    σ_tree = decode(read(PoS_REPLY))
    σ = convert(Tuple{BigInt, Vector{BigInt}, BigInt, BigInt, Vector{BigInt}, BigInt}, σ_tree)

    return VShuffleProof(μ, τ, σ)
end

function load_verificatum_simulator(basedir::AbstractString; auxsid = "default")

    spec = ProtocolSpec(basedir * "/protInfo.xml"; auxsid)

    proposition = load_verificatum_proposition(basedir, auxsid)

    NIZKP = basedir * "/dir/nizkp/$auxsid/"
    proof = load_verificatum_proof("$NIZKP/proofs/", proposition.g)

    simulator = Simulator(proposition, proof, spec)

    return simulator
end

### The simulator type will deal with loading the data.

struct VInit{G<:Group} #<: Verifier
    spec::ProtocolSpec{G}
    proposition::Shuffle{G} # ADD G!
    ρ::Vector{UInt8}
    𝐡::Vector{G}
end

leaf(x::String) = encode(Leaf(x))

function gen_verificatum_basis(::Type{G}, prghash::Hash, rohash::Hash, N::Integer; nr::Integer = 0, ρ = UInt8[], d = [ρ..., leaf("generators")...]) where G <: Group

    roprg = ROPRG(d, rohash, prghash)
    prg = roprg(UInt8[]) # d is a better argument than x

    return rand(prg, G, N; nr)
end

function VInit(spec::ProtocolSpec{G}, proposition::Shuffle) where G <: Group

    ρ = ro_prefix(spec) ### I can add another method there

    𝔀 = proposition.𝐞
    N = length(𝔀)

    (; g, nr, rohash, prghash)  = spec

    𝐡 = gen_verificatum_basis(G, prghash, rohash, N; nr, ρ)

    return VInit(spec, proposition, ρ, 𝐡)
end

struct VPermCommit{G<:Group} #<: Verifier
    spec::ProtocolSpec{G}
    proposition::Shuffle{G}
    ρ::Vector{UInt8}
    𝐡::Vector{G}
    s::Vector{UInt8}
    𝐞::Vector{BigInt}
end

function VPermCommit(v::VInit{G}, 𝐮::Vector{G}) where G <: Group
    (; 𝐡, ρ, spec, proposition) = v
    (; ne, prghash, rohash) = spec
    𝔀, 𝔀′ = proposition.𝐞, proposition.𝐞′
    (; g, pk) = proposition

    N = length(𝔀)

    roprg = ROPRG(ρ, rohash, prghash)

    pk_tree = (g, pk)
    tree = Tree((g, 𝐡, 𝐮, pk_tree, 𝔀, 𝔀′))
    prg = roprg(encode(tree))

    (; s) = prg

    𝐭 = rand(prg, BigInt, N; n = ne)
    𝐞 = mod.(𝐭, BigInt(2)^ne)

    return VPermCommit(spec, proposition, ρ, 𝐡, s, 𝐞)
end

struct VPoSCommit{G<:Group} #<: Verifier
    spec::ProtocolSpec{G}
    proposition::Shuffle{G} # ADD G!
    ρ::Vector{UInt8}
    𝐡::Vector{G}
    𝐞::Vector{BigInt}
    𝓿::BigInt
end

function VPoSCommit(v::VPermCommit{G}, τ::Tuple{Vector{G}, G, Vector{G}, G, G, Tuple{G, G}}) where G <: Group
    (; 𝐡, ρ, 𝐞, spec, proposition, s) = v
    (; nv, rohash) = spec

    ro_challenge = RO(rohash, nv)
    tree_challenge = Tree((Leaf(s), τ))
    𝓿 = interpret(BigInt, ro_challenge([ρ..., encode(tree_challenge)...]))

    return VPoSCommit(spec, proposition, ρ, 𝐡, 𝐞, 𝓿)
end

function VPoSCommit(v::VPermCommit{G}, 𝐜̂::Vector{G}, t::Tuple{G, G, G, Tuple{G, G}, Vector{G}}) where G <: Group
    (t₁, t₂, t₃, t₄, 𝐭̂) = t
    𝐁 = 𝐜̂

    𝐁′= 𝐭̂
    D′ = t₂
    A′ = t₃
    C′ = t₁
    F′ = t₄

    τ = (𝐁, A′, 𝐁′, C′, D′, F′)

    return VPoSCommit(v, τ)
end

PoSChallenge(verifier::VPoSCommit) = PoSChallenge(verifier.𝐡, verifier.𝐞, verifier.𝓿)

function verify(proposition::Shuffle, proof::VShuffleProof, challenge::PoSChallenge; verbose=false)

    𝐡, 𝐞, 𝓿 = challenge.𝐡, challenge.𝐮, challenge.c
    𝔀, 𝔀′ = proposition.𝐞, proposition.𝐞′
    (; g, pk) = proposition

    (; μ, τ, σ) = proof

    𝐮 = μ
    𝐁, A′, 𝐁′, C′, D′, F′ = τ
    k_A, 𝐤_B, k_C, k_D, 𝐤_E, k_F = σ

    N = length(𝔀)

    A = prod(𝐮 .^ 𝐞)

    C = prod(𝐮) / prod(𝐡)
    D = 𝐁[N] * inv(𝐡[1])^prod(𝐞)

    F = ∏(𝔀 .^ 𝐞)

    report = Report()

    report &= "A", A^𝓿 * A′ == g^k_A * prod(𝐡 .^ 𝐤_E)
    report &= "C", C^𝓿 * C′ == g^k_C
    report &= "D", D^𝓿 * D′ == g^k_D

    report &= "B", Bool[
        𝐁[1]^𝓿 * 𝐁′[1] == g^𝐤_B[1] * 𝐡[1]^𝐤_E[1],
        (𝐁[i]^𝓿 * 𝐁′[i] == g^𝐤_B[i] * 𝐁[i - 1]^𝐤_E[i] for i in 2:N)...
    ]

    enc = Enc(pk, g)
    report &= "F", F^𝓿 * F′ == enc(-k_F) * ∏(𝔀′ .^ 𝐤_E)

    if verbose || isvalid(report) == false
        println(report)
    end

    return isvalid(report)
end

function verify(proposition::Shuffle, proof::VShuffleProof, verifier::ProtocolSpec)

    v1 = VInit(verifier, proposition)

    (; μ) = proof
    v2 = VPermCommit(v1, μ)

    (; τ) = proof
    v3 = VPoSCommit(v2, τ)

    v4 = PoSChallenge(v3)

    return verify(proposition, proof, v4)
end

step(spec::ProtocolSpec, proposition::Proposition) = VInit(spec, proposition)
step(v::VInit{G}, 𝐜::Vector{G}) where G <: Group = VPermCommit(v, 𝐜)
step(v::VPermCommit, 𝐜̂, t) = VPoSCommit(v, 𝐜̂, t)

challenge(v::VInit) = (v.𝐡, v.𝐡[1])
challenge(v::VPermCommit) = v.𝐞
challenge(v::VPoSCommit) = v.𝓿