text:
abstract type AbstractMessageType <: AbstractDict{AbstractString,Any} end
const MessageType = Dict{String,Any}

struct StratumError <: Exception
    code::Int
    msg::AbstractString
    data::Any
end

function Base.showerror(io::IO, ex::StratumError)
    error_code_as_string =
    # JSON RPC 2.0 Error Codes
    # https://www.jsonrpc.org/specification#error_object
        if ex.code == -32700
            "ParseError"
        elseif ex.code == -32600
            "InvalidRequest"
        elseif ex.code == -32601
            "MethodNotFound"
        elseif ex.code == -32602
            "InvalidParams"
        elseif ex.code == -32603
            "InternalError"

            # Stratum V1 Error Codes
            # https://braiins.com/stratum-v1/docs
        elseif ex.code == 20
            "Unknown"
        elseif ex.code == 21
            "JobNotFound"
        elseif ex.code == 22
            "DuplicateShare"
        elseif ex.code == 23
            "LowDifficultyShare"
        elseif ex.code == 24
            "UnauthorizedWorker"
        elseif ex.code == 25
            "NotSubscribed"

        else
            "Unknown"
        end

    print(io, error_code_as_string)
    print(io, ": ")
    print(io, ex.msg)
    if ex.data !== nothing
        print(io, " (")
        print(io, ex.data)
        print(io, ")")
    end
end

mutable struct StratumEndpoint{IOIn<:IO,IOOut<:IO,MsgIn<:MessageType,MsgOut<:MessageType}
    pipe_in::IOIn
    pipe_out::IOOut

    name::Union{Nothing,String}

    err_handler::Union{Nothing,Function}

    msg_queue_in::Channel{MsgIn}
    msg_queue_out::Channel{MsgOut}

    outstanding_requests::Dict{String,Channel{MessageType}}

    status::Symbol

    reader::Union{Nothing,Task}
    writer::Union{Nothing,Task}
end

function StratumEndpoint(pipe_in, pipe_out; name=nothing, err_handler=nothing)
    return StratumEndpoint(
        pipe_in,
        pipe_out,
        name,
        err_handler,
        Channel{MessageType}(Inf),
        Channel{MessageType}(Inf),
        Dict{String,Channel{MessageType}}(),
        :idle,
        nothing,
        nothing,
    )
end

Base.isopen(x::StratumEndpoint) = x.status === :running && iswritable(x) && isreadable(x)

function check_isopen(x::StratumEndpoint)
    isopen(x) || error("Endpoint is not open")
    return nothing
end

Base.iswritable(x::StratumEndpoint) = isopen(x.msg_queue_out) && isopen(x.pipe_out)
Base.isreadable(x::StratumEndpoint) = isopen(x.msg_queue_in)  && isopen(x.pipe_in)

check_writable(x::StratumEndpoint) = iswritable(x) || error("Endpoint is not writable")
check_readable(x::StratumEndpoint) = isreadable(x) || error("Endpoint is not readable")

function Base.take!(x::StratumEndpoint)
    check_isopen(x)
    msg = take!(x.msg_queue_in)
    return msg
end

function Base.put!(x::StratumEndpoint, msg::AbstractDict)
    check_isopen(x)
    put!(x.msg_queue_out, msg)
    return x
end

function write_transport_layer(x::StratumEndpoint, msg::AbstractString)
    msg_utf8 = transcode(UInt8, msg)
    n = write(x.pipe_out, msg_utf8, '\n')
    flush(x.pipe_out)
    return n
end

function read_transport_layer(x::StratumEndpoint)
    msg = chomp(readline(x.pipe_in))
    return msg
end

function Base.run(x::StratumEndpoint)
    x.status == :idle   || error("Endpoint is not idle.")
    x.status == :closed && error("Endpoint is closed.")

    x.status = :running

    x.writer = @async try
        try
            while true
                msg = take!(x.msg_queue_out)
                isopen(x) || break
                write_transport_layer(x, JSON.json(msg))
                x.name === nothing ? (@debug "Sent:" msg) : (@debug "Sent ($(x.name)):" msg)
            end
        finally
            x.name === nothing ? (@debug "Closed") : (@debug "Closed ($(x.name)):")
            close(x.msg_queue_out)
            close(x.pipe_out)
        end
    catch err
        bt = catch_backtrace()
        handle_error(x, err, bt)
    end

    x.reader = @async try
        try
            while true
                msg_raw = read_transport_layer(x)
                (isopen(x) && !isempty(msg_raw)) || break

                msg = JSON.parse(msg_raw)
                x.name === nothing ? (@debug "Received:" msg) : (@debug "Received ($(x.name)):" msg)

                id = get(msg, "id", nothing)
                if haskey(x.outstanding_requests, id)
                    put!(x.outstanding_requests[id], msg)
                else
                    put!(x.msg_queue_in, msg)
                end
            end
        finally
            x.name === nothing ? (@debug "Closed") : (@debug "Closed ($(x.name)):")
            close(x.msg_queue_in)
            close(x.pipe_in)
        end
    catch err
        bt = catch_backtrace()
        handle_error(x, err, bt)
    end

    return x
end

function handle_error(x::StratumEndpoint, err, bt)
    if x.err_handler !== nothing
        x.err_handler(err, bt)
    else
        Base.display_error(stderr, err, bt)
    end
end

function send_notification(x::StratumEndpoint, method::AbstractString, params)
    check_isopen(x)

    msg = Dict("jsonrpc" => "2.0", "method" => method, "params" => params)

    put!(x, msg)

    return nothing
end

function send_request(x::StratumEndpoint, method::AbstractString, params)
    check_isopen(x)

    id = string(UUIDs.uuid4())
    msg = Dict("jsonrpc" => "2.0", "id" => id, "method" => method, "params" => params)

    x.outstanding_requests[id] = Channel{MessageType}(1)

    put!(x, msg)

    response = take!(x.outstanding_requests[id])
    delete!(x.outstanding_requests, id)

    if haskey(response, "result")
        # TODO which makes the most sense?
        # return response["result"]
        return response
    elseif haskey(response, "error")
        err_code, err_msg, err_data = response["error"]
        throw(StratumError(err_code, err_msg, err_data))
    else
        throw(StratumError(0, "ERROR AT THE TRANSPORT LEVEL", nothing))
    end
end

function Base.iterate(x::StratumEndpoint, state=nothing)
    check_isopen(x)
    try
        return take!(x), nothing
    catch err
        if err isa InvalidStateException
            return nothing
        else
            rethrow(err)
        end
    end
end

function send_success_response(x, original_request, result)
    msg = Dict("jsonrpc" => "2.0", "id" => original_request["id"], "result" => result)
    return put!(x, msg)
end

function send_error_response(x, original_request, err_code, err_msg, err_data)
    msg = Dict(
        "jsonrpc" => "2.0",
        "id" => original_request["id"],
        "error" => [err_code, err_msg, err_data],
    )
    return put!(x, msg)
end

function Base.close(x::StratumEndpoint)
    flush(x)

    x.status = :closed

    close(x.msg_queue_in)
    close(x.msg_queue_out)

    close(x.pipe_in)
    close(x.pipe_out)

    fetch(x.writer)
    fetch(x.reader)

    return nothing
end

function Base.flush(x::StratumEndpoint)
    check_isopen(x)

    while isready(x.msg_queue_out)
        yield()
    end
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 8bd459cb-20bb-483e-a849-e18caae3beef
begin
    using Pkg
	Pkg.activate(joinpath(Pkg.devdir(), "MLCourse"))
    using MLCourse, MLJ, DataFrames, MLJMultivariateStatsInterface, OpenML,
    Plots, LinearAlgebra, Statistics, Random, CSV, MLJLinearModels
    chop(x, eps = 1e-12) = abs(x) > eps ? x : zero(typeof(x))
end

# ╔═╡ 1c3aa1be-c58a-4eb9-a192-ddfffa8a2c56
using TSne

# ╔═╡ b97724e4-d7b0-4085-b88e-eb3c5bcbe441
begin
    using PlutoUI
    PlutoUI.TableOfContents()
end

# ╔═╡ 31bcc3f0-88c1-4035-b0de-b4217623dc22
md"# Reconstructing a Generative Model with PCA

Let us assume a given data set ``X`` with ``p`` predictors was generated by some mixture of ``q \leq p`` hidden (=unobserved) predictors. Under some conditions PCA can recover these unobserved predictors from the data set ``X``.

## An Artificial Toy Example

In this example we generate 2-dimensional (hidden = unobserved) data and project it to 3 dimensions. This 3-dimensional input is our measured data that we take as input for PCA. We will find that PCA recovers in this case the original data."

# ╔═╡ e24cff15-a966-405a-b083-4b86be09c6f1
begin
    Random.seed!(17)
    h = [1.5 .5] .* randn(500, 2) # 2D data with higher variance along first dimension
    G = [-.3 .4 .866
         .885 .46 .1]             # projection matrix to 3D
    x = h * G                     # projected data
end;

# ╔═╡ 9e60bd23-bd0e-4311-a059-a39a70143307
let
	plotly()
	p1 = scatter(h[:, 1], h[:, 2], aspect_ratio = 1,
                 xlabel = "H₁", ylabel = "H₂",
                 label = "h", title = "original (hidden) 2D data")
	p2 = scatter3d(x[:, 1], x[:, 2], x[:, 3], markersize = 1,
                   xlabel = "X₁", ylabel = "X₂", zlabel = "X₃",
                   title = "measured 3D data")
	plot(p1, p2, layout = (1, 2), legend = false)
end

# ╔═╡ a0bb0a63-67e2-414e-a3f0-b76875c1295d
m1 = fit!(machine(PCA(), x));

# ╔═╡ 4439f803-e70f-4231-82e9-88bb574498f4
ĥ = MLJ.transform(m1, x);

# ╔═╡ e38cdfcf-362f-4e6e-81a2-e58334bb8d5f
let
    gr()
    scatter(h[:, 1], h[:, 2],
            xlabel = "H₁", ylabel = "H₂",
            aspect_ratio = 1, label = "original hidden data")
    scatter!(ĥ.x1, ĥ.x2, label = "reconstructed hidden data", markersize = 3)
end

# ╔═╡ 16843059-269e-4557-a4b9-5d7034f99ba8
md"We can see that the reconstructed hidden data almost perfectly matches the true hidden data. We can also see that the transposed of the reconstructed projection matrix `fitted_params(m1).projection` is close to the projection matrix `G` we used."

# ╔═╡ 6e033717-5815-4ffe-b991-81b72b99df9f
(fitted_params(m1).projection', G)

# ╔═╡ b25687fc-aa11-4843-9a89-31c4c9b5bc23
md"This example is specifically constructed for PCA to succeed. The conditions on which success relies are
1. the original data set has the largest variance along ``H_1``. You can change the data set such that it has higher variance along dimension ``H_2``. What happens in this case?
2. The rows of the projection matrix `G` have norm 1. What happens if you multiply the original `G` by 2?
3. `x` was obtained by linear transformation with `G` from `h`.
"

# ╔═╡ 14bfdb2d-fd71-4c62-8823-6e292a424bf5
md"## Reconstructing Spike Trains

Here is a more realistic - but still artificial - example. Let us assume we were interested in the activity of 2 neurons somewhere in a brain, firing at different moments in time an action potential with a shape that is characteristic for the neuron (see blue and red curve in the figure below). Each one of three electrodes that are placed in the brain tissue near these two neurons sense the two neurons in a way that depends on the relative positions of the electrodes and the neurons. The goal is to reconstruct from the measurements at the electrodes the activity of the neurons.
"

# ╔═╡ f6117412-f6d8-41d7-9524-dc25fc891deb
begin
    s1(t) = t < 0 ? 0 : (exp(-t/.2) - exp(-t/.01)) * sin(10*t)
    s2(t) = t < 0 ? 0 : (exp(-t/.5) - exp(-t/.3)) * sin(4*t)
end;

# ╔═╡ 9da3a52e-23ca-4ce4-8524-0abeb4be2306
let
    gr()
    t = -1:.02:3
    plot(t, s1.(t), label = "spike of neuron 1",
         xlabel = "time [arbitrary units]", ylabel = "voltage [arbitrary units]")
    plot!(t, s2.(t), label = "spike of neuron 2",
         xlabel = "time [arbitrary units]", ylabel = "voltage [arbitrary units]")
end

# ╔═╡ 00e15b1c-a12d-4ac2-b3c6-6fece1abef31
begin
    v1 = sum(s1.((0:.02:1000) .- ts) for ts in unique(rand(0:2:1000, 200))) .+ 1e-2*randn(50001)
	v2 = sum(s2.((0:.02:1000) .- ts) for ts in unique(rand(0:2:1000, 200))) .+ 1e-2*randn(50001)
end;

# ╔═╡ 60ba84d8-4f48-4d8a-b695-11edbe91ccc8
let
    gr()
    p1 = plot(v1, xlim = (0, 1000), title = "unobserved spike train",
              label = "v1", xlabel = "time", ylabel = "voltage")
    p2 = plot(v2, xlim = (0, 1000), label = "v2", xlabel = "time", ylabel = "voltage")
    p3 = scatter(v1, v2, aspect_ratio = 1,
                 label = nothing, xlabel = "v1", ylabel = "v2")
    plot(plot(p1, p2, layout = (2, 1)), p3, layout = (1, 2))
end

# ╔═╡ b256ff25-431a-4c47-8195-d958d3d28c56
sdata = [v2 v1] * G;

# ╔═╡ abb36d40-ad0e-442c-ab87-b15837a2541c
md"In the figure above you see on the left some time steps of our artificially generated neural activity patterns. Each dot in the figure on the right shows the voltages of the two neurons at a given moment in time.
We assume now that the electrodes measure `sdata`, which is obtained by mixing the originial data in three different ways; think of the 3 columns of `G` as the mixing coefficients. In reality the dependence of the electrode signal on the neural activity might not be as simple as we assume here."

# ╔═╡ e81b5dee-dffb-4a4c-b471-0974599fd87d
let
    gr()
    p1 = plot(sdata[:, 1], xlim = (0, 1000), title = "measured signal",
              label = "electrode 1",
              xlabel = "time", ylabel = "voltage", c = :black)
    p2 = plot(sdata[:, 2], xlim = (0, 1000), label = "electrode 2",
              xlabel = "time", ylabel = "voltage", c = :black)
    p3 = plot(sdata[:, 3], xlim = (0, 1000), label = "electrode 3",
              xlabel = "time", ylabel = "voltage", c = :black)
    p4 = scatter3d(sdata[:, 1], sdata[:, 2], sdata[:, 3],
                   xlabel = "electrode 1" , ylabel = "electrode 2",
                   zlabel = "electode 3", label = nothing, c = :black)
    plot(plot(p1, p2, p3, layout = (3, 1)),
         p4, layout = (1, 2), size = (700, 500))
end

# ╔═╡ 48528bf7-8a40-4773-944c-54d0446a5cac
sma = fit!(machine(PCA(mean = 0), sdata));

# ╔═╡ 9a5e0976-1c0e-45fc-855e-b0fbec801e8c
srec = MLJ.transform(sma, sdata);

# ╔═╡ b5d2986e-3723-4463-ab05-cc176aea60d8
let
	gr()
	p1 = plot(v1, label = "original v1")
	plot!(srec.x1, label = "reconstruction v1", xlim = (0, 1000))
	p2 = plot(v2, label = "original v2")
	plot!(srec.x2, label = "reconstruction v2", xlim = (0, 1000))
    p3 = scatter(v1, v2, label = "original")
	scatter!(srec.x1, srec.x2, aspect_ratio = 1, label = "reconstruction")
	plot(plot(p1, p2, layout = (2, 1)), p3, layout = (1, 2))
end

# ╔═╡ eb421c10-a023-4faa-bcfb-374b8ebb9c38
md"Also here we see that PCA recovers the original activity patterns. Note that the almost perfect reconstruction relies again on the three conditions specified in the previous section."

# ╔═╡ c1262304-6846-4811-ae82-397863255415
md"""# PCA Finds the Directions of Largest Variance

In the figure below you see a visualization of the variance along the first, second or third principal component (loading vectors ``\phi_1, \phi_2, \phi_3``). The colored lines are parallel to the loading vector you pick in the dropdown menu. They go from each data point to the "center"-plane that goes through the origin ``(0, 0, 0)`` and stands perpendicular to the picked component (loading vector). **The length of the line that starts at data point ``i`` is given by the score ``z_{ij}`` for principal component ``j``**. The average squared length of these lines is the variance along the chosen principal component.
"""

# ╔═╡ f9a506bb-6bdf-48aa-8f81-8fb6de078533
@bind pc Select(["1 PC", "2 PC", "3 PC"])

# ╔═╡ 2b73abbb-8cef-4da0-b4cf-b640fe0ad102
md"The `PCA` function in `MLJ` has three important keyword arguments as explained in the doc string (see Live docs). If the keyword argument `pratio = 1` all principal components are computed."

# ╔═╡ 45ee714c-9b73-492e-9422-69f27be79de2
md"The report of a `PCA` machine returns the number of dimensions `indim` of the data. the number of computed principal components `outdim`, the total variance `tvar` and the `mean` of the data set, together with the variances `principalvars` along the computed principal components. The sum of the variances `principalvars` is given in `tprincipalvar` and the rest in `tresidualvar = tvar - tprincipalvar`. If you change the `pratio` to a lower value above you will see that not all principal components are computed."

# ╔═╡ ecc77003-c138-4ea9-93a7-27df46aa1e89
md"Biplots are useful to see loadings and scores at the same time. The scores of each data point are given by the coordinates of the numbers plotted in gray when reading the bottom and the left axis. For example, the score for point 22 would be approximately (0, 4.8). The loadings are given by the red labels (not the tips of the arrows) when reading the top and the right axis. For examples the loadings of the first principal component are approximately (0.95, 0.1, 0.28) (see also `fitted_params` above)."

# ╔═╡ 18ad4980-2960-4ea9-8c9a-bf17cd3546ff
function data_generator(; seed = 71,  N = 100, rng = MersenneTwister(seed))
	h = [2.5 1.5 .8] .* randn(rng, N, 3)
    a1 = rand(); a2 = rand()
	G = [1 0 0
         0 cos(a1) -sin(a1)
         0 sin(a1) cos(a1)] * [cos(a2) 0 sin(a2)
                               0       1   0
                               -sin(a2) 0 cos(a2)]
	x = h * G
end;

# ╔═╡ b2688c20-a578-445f-a5c4-335c85931862
pca_data = data_generator(seed = 241);

# ╔═╡ 927627e0-9614-4234-823d-eb2e13229784
pca = fit!(machine(PCA(pratio = 1), DataFrame(pca_data, :auto)));

# ╔═╡ c7935f29-21f3-414a-b013-a1abd0f9f502
fitted_params(pca) # shows the loadings as columns

# ╔═╡ b52d72e3-64a3-4b7e-901b-3512b97aec23
report(pca)

# ╔═╡ 89406c26-7cba-4dd2-a6a3-8d572b440e01
biplot(pca)

# ╔═╡ 743477bb-6cf0-419e-814b-e774738e8d89
m2 = fit!(machine(PCA(maxoutdim = 2), pca_data));

# ╔═╡ 4a8188a0-c809-4014-9fef-cf6b5f830b1b
begin
    function pc_vectors(m)
        proj = fitted_params(m).projection
        a1 = proj[:, 1]
        a2 = proj[:, 2]
        a3 = cross(a1, a2)
        a3 ./= norm(a3)
        (a1, a2, a3)
    end
    function pca_plane(x, y, m)
        _, _, a3 = pc_vectors(m)
        a3 ./= a3[3]
        -a3[1]*x-a3[2]*y
    end
    pca_plane(x, y) = pca_plane(x, y, m2)
end;

# ╔═╡ 884a7428-2dd4-4091-b889-3acfbfa45d5b
let
    plotly()
    scatter3d(pca_data[:, 1], pca_data[:, 2], pca_data[:, 3],
              markersize = 1, label = nothing)
    a = pc_vectors(m2)
    col = [:red, :green, :orange]
    k = Meta.parse(pc[1:1])
    var = 0.
    for i in 1:size(pca_data, 1)
        s = pca_data[i, :]
        e = s - (a[k]'*s) * a[k]
        plot3d!([s[1], e[1]], [s[2], e[2]], [s[3], e[3]],
                color = col[k], w = 4, label = nothing)
    end
    plot3d!(xlim = (-7, 7), ylim = (-7, 7), zlim = (-7, 7), aspect_ratio = 1)
end

# ╔═╡ df9a5e8d-0692-44e4-8fd0-43e8643c4b7c
md"# PCA Finds Linear Subspaces Closest to the Data

The `fitted_params` of a `PCA` machine contain as columns the the principal component vectors. The first two components span a plane (see figure below). This plane is closest to the data in the sense that the average squared distance to each data point is minimal. There is no other plane that would be closer to the data points."

# ╔═╡ d3e5ca43-f70e-4a11-a3c9-fc0755179e4c
let
    plotly()
    scatter3d(pca_data[:, 1], pca_data[:, 2], pca_data[:, 3], markersize = 1, label = "data")
    a1, a2, a3 = pc_vectors(m2)
    a1 *= 2; a2 *= 2
    wireframe!(-8:.1:8, -8:.1:8, pca_plane,
               label = nothing, w = 1)
    plot3d!([-a1[1], a1[1]], [-a1[2], a1[2]], [-a1[3], a1[3]],
            c = :red, w = 10, label = "PC1")
    plot3d!([-a2[1], a2[1]], [-a2[2], a2[2]], [-a2[3], a2[3]],
            c = :green, arrow = true, w = 10, label = "PC2")
    for i in 1:size(pca_data, 1)
        s = pca_data[i, :]
        e = s - (a3'*s) * a3
        plot3d!([s[1], e[1]], [s[2], e[2]], [s[3], e[3]],
                color = :orange, label = i == 1 ? "distance to plane" : nothing, w = 4)
    end
    plot3d!()
end

# ╔═╡ 936570b9-43da-4c32-b15a-d595f690b1a6
md"
## Reconstructions of the Data

With ``X`` the ``100\times 3``-dimensional matrix of the original data, the reconstruction is given by ``X\Phi_L\Phi_L^T = Z_L\Phi_L^T`` with ``\Phi_L`` the ``3\times L``-dimensional matrix constructed from the first ``L`` loading vectors as columns. The reconstructions are closest to the original data while being forced to vary only along the first ``L`` principal components.

Reconstruction using the first $(@bind firstkpc Select(string.(1:3))) component(s).
"

# ╔═╡ 9032ad56-123c-4585-9b67-c94e5ecb3899
let
    plotly()
    scatter3d(pca_data[:, 1], pca_data[:, 2], pca_data[:, 3],
              markersize = 1, label = "data")
    a1, a2, a3 = pc_vectors(m2)
    ϕ = if firstkpc == "1"
        a1
    elseif firstkpc == "2"
        [a1 a2]
    else
        [a1 a2 a3]
    end
    rec = pca_data * ϕ * ϕ'
    scatter3d!(rec[:, 1], rec[:, 2], rec[:, 3], legend = false,
               markersize = 1, label = "reconstruction")
end

# ╔═╡ a8188f31-204a-406f-96a2-86768918a2f9
md"## PCA for Visualization"

# ╔═╡ 265594dd-ba59-408a-925a-c2a173343df7
wine = OpenML.load(187) |> DataFrame;

# ╔═╡ 51720585-09c8-40bb-8ec2-c2f25804d731
X = float.(select(wine, Not(:class)))

# ╔═╡ d45faa6b-706b-4990-80d9-32acc001b242
mwine = fit!(machine(@pipeline(Standardizer(), PCA()), X));

# ╔═╡ 6c8c86f0-02ec-4ce7-9d16-1b4865c344df
let
	gr()
	biplot(mwine)
end

# ╔═╡ e72c09e3-91cb-4455-9880-628428730efe
md"### Scaling the Data Matters!!!"

# ╔═╡ 9f6b1079-209f-4a24-96c8-81bdcf8bfc32
mwine2 = machine(PCA(pratio = 1), X) |> fit!;

# ╔═╡ 7eac92b8-0639-46b8-966e-cbf44de7ee86
let
	gr()
	biplot(mwine2)
end

# ╔═╡ ab17ef21-ea29-4d62-b14e-b87cffc0a1f9
md"# PCA Decorrelates

A principal component analysis can be performed by computing the covariance matrix ``X^T X`` and its eigenvalues. In the cell below you can see that the covariance matrix is not diagonal and that the eigenvalues (divided by ``n-1``) correspond to the variances along the different principal components.
"

# ╔═╡ da99f084-522b-4969-a65c-31d1548d3073
let
	X = pca_data .- mean(pca_data, dims = 1) # subtract the mean of each column
	(XtX = X'X, eigen_values = eigvals(X'X))
end

# ╔═╡ d88509c3-51ca-4ff9-9c1f-c3e13ec0e133
report(pca).principalvars * (size(pca_data, 1) - 1)

# ╔═╡ 615a823f-32df-469d-b2a3-e9d2c35e36a7
md"An interesting property of the scores is that their covariance matrix is diagonal, with the eigenvalues of ``X^TX`` on the diagonal. Hence we call PCA a decorrelating transformation."

# ╔═╡ ded0e6cd-e969-47d5-81cf-e2a967bc0a34
let
    mach = fit!(machine(PCA(pratio = 1), pca_data))
    Z = MLJ.matrix(MLJ.transform(mach, pca_data))
    (ZtZ = chop.(Z'*Z),) # chop sets values close to zero to zero.
end

# ╔═╡ 9724b203-ae86-4cb6-afec-68e26c01ab83
md"## PCA for Denoising"

# ╔═╡ 715a0e0c-5105-4fee-8b93-03a966a1bd8e
data = let t = 0:.1:8
    Random.seed!(1)
    DataFrame(hcat([rand((sin, cos)).(t) .+ .2*randn(length(t))
				    for _ in 1:2000]...)', :auto)
end;

# ╔═╡ 95ba7abd-c425-46d4-8382-fedec7ed58f6
let
    gr()
	plot(Array(data[1, :]))
	for i in 2:20
		plot!(Array(data[i, :]))
	end
	plot!(legend = false, xlabel = "time", ylabel = "voltage", title = "raw data")
end

# ╔═╡ d2bfdda4-4496-4daa-bcb5-51d9a7519823
mdenoise = fit!(machine(PCA(maxoutdim = 2, mean = 0), data));

# ╔═╡ 694ad552-3088-4872-bed5-89b4dfc67219
report(mdenoise)

# ╔═╡ c54b4211-36aa-4d7c-8e7a-c67cf5bdabdc
reconstruction = MLJ.inverse_transform(mdenoise, MLJ.transform(mdenoise, data));

# ╔═╡ b9383ef0-8499-4e56-a8d6-606ec128676c
let
	plot(Array(reconstruction[1, :]))
	for i in 2:20
		plot!(Array(reconstruction[i, :]))
	end
# 	plot!(s1.(0:.02:3), w = 2)
# 	plot!(s2.(0:.02:3), w = 2)
# 	plot!(s3.(0:.02:3), w = 2)
    plot!(legend = false, xlabel = "time", ylabel = "voltage",
          title = "denoised reconstruction")
end

# ╔═╡ af03d989-b19f-44ac-a275-d8af57bbeaeb
md"## PCA for Compression"

# ╔═╡ 026a362d-3234-409b-b9cb-e6a2a50c6a0d
faces = OpenML.load(41083) |> DataFrame;

# ╔═╡ e967a9cb-ad95-4bd2-8fbd-491dc6fe0475
let
    gr()
    plot(Gray.(vcat([[hcat([[reshape(Array(faces[j*8 + i, 1:end-1]), 64, 64)' ones(64)] for i in 1:8]...); ones(8*65)'] for j in 0:5]...)), ticks = false)
end

# ╔═╡ 0cb282aa-81d4-408d-a919-0d28cdbb7bed
pca_faces = fit!(machine(PCA(), faces[:, 1:end-1]));

# ╔═╡ c70b5212-a823-4ff8-937d-10919efc766c
let
    vars = report(pca_faces).principalvars ./ report(pca_faces).tvar
    p1 = plot(vars, label = nothing, yscale = :log10,
              xlabel = "component", ylabel = "proportion of variance explained")
    p2 = plot(cumsum(vars),
              label = nothing, xlabel = "component",
              ylabel = "cumulative prop. of variance explained")
    plot(p1, p2, layout = (1, 2), size = (700, 400))
end

# ╔═╡ 0742ea40-1895-48f4-b149-d472397e83d0
let
    pcs = fitted_params(pca_faces).projection
    scale = x -> (x .- minimum(x)) ./ (maximum(x) - minimum(x))
    plots = []
    for i in [1:5; size(pcs, 2)-4:size(pcs, 2)]
        push!(plots, plot(Gray.(scale(reshape(pcs[:, i], 64, 64))'),
                          title = "PC $i"))
    end
    plot(plots..., layout = (2, 5), showaxis = false, ticks = false)
end

# ╔═╡ 7c020b97-384d-43fc-af1c-a4bf75a6f06c
md"cumulative proportion of variance explained $(@bind pratio Slider(.9:.001:.999, show_value = true))"

# ╔═╡ cddfed10-7c9c-47b0-83fc-b3322da522ed
begin
    first_few_pc_faces = fit!(MLJ.machine(PCA(pratio = pratio), faces[:, 1:end-1]))
    npc = size(fitted_params(first_few_pc_faces).projection, 2)
    md"Using the first ``L = ``$npc components. The compression ratio is ``\frac{n\times p}{(p+n)\times L}`` = $(round(400*64^2/((64^2 + 400)*npc), sigdigits = 2))."
end

# ╔═╡ 03f92562-7419-4d8e-baa2-67f5b021219e
md"Image number $(@bind imgid Select(string.(1:nrow(faces))))"

# ╔═╡ 8e466825-0c84-41c7-93c6-e3cc81a60ef5
let
    gr()
    i = Meta.parse(imgid)
    p1 = plot(Gray.(reshape(Array(faces[i, 1:end-1]), 64, 64)'),
              ticks = false, title = "orginal")
    p2 = plot(Gray.(reshape(Array(MLJ.inverse_transform(first_few_pc_faces,
                                                       MLJ.transform(first_few_pc_faces,
                                                                     faces[i:i,1:end-1]))), 64, 64)'), ticks = false, title = "reconstruction")
    plot(p1, p2, layout = (1, 2), size = (700, 500))
end

# ╔═╡ 218d4c11-45a1-4d92-8b48-d08251804638
md"## PCA Applied to the Weather Data"

# ╔═╡ 0943f000-a148-4470-ae5e-9d8d39cd8207
begin
    weather = CSV.read(joinpath(@__DIR__, "..", "data", "weather2015-2018.csv"),
                       DataFrame);
    weather = select(weather, (x -> match(r"direction", x) == nothing && match(r"time", x) == nothing).(names(weather)))
end

# ╔═╡ 2e86337e-9fb8-4dd9-ac96-c1177e7ebd16
weather_mach = fit!(machine(@pipeline(Standardizer(), PCA()), weather));

# ╔═╡ 24179171-3992-4253-88a5-15ccdbb26e13
let
    gr()
    p1 = biplot(weather_mach)
    p2 = biplot(weather_mach, pc = (3, 4))
    plot(p1, p2, layout = (1, 2), size = (700, 400))
end

# ╔═╡ 8bc5f5b1-cda1-4018-a215-5a1900c2a9d2
let
    vars = report(weather_mach).pca.principalvars ./ report(weather_mach).pca.tprincipalvar
    p1 = plot(vars, label = nothing, yscale = :log10,
              xlabel = "component", ylabel = "proportion of variance explained")
    p2 = plot(cumsum(vars),
              label = nothing, xlabel = "component",
              ylabel = "cumulative prop. of variance explained")
    plot(p1, p2, layout = (1, 2), size = (700, 400))
end

# ╔═╡ d44be4c6-a409-421a-8e68-39fa752db4c0
md"# Limitations of PCA"

# ╔═╡ d3b769d7-ec71-4496-926a-e838c3b50c2a
begin
    Random.seed!(123)
    eight_clouds = DataFrame(vcat([[.5 .5 20] .* randn(50, 3) .+ [x y 0]
                             for y in (-2, 2), x in (-6, -2, 2, 6)]...), :auto);
    eight_clouds = MLJ.transform(fit!(machine(Standardizer(), eight_clouds)))
end;

# ╔═╡ 67395b94-46bd-43e5-9e12-0bd921de8203
let
    gr()
    c = vcat([fill(i, 50) for i in 1:8]...)
    p0 = scatter3d(eight_clouds.x1, eight_clouds.x2, eight_clouds.x3,
                   c = c, title = "original 3D data")
    p1 = scatter(eight_clouds.x1, eight_clouds.x2, c = c,
                 xlabel = "x", ylabel = "y", title = "projection onto X-Y")
    pc = MLJ.transform(fit!(machine(PCA(maxoutdim = 2), eight_clouds)))
    p2 = scatter(pc.x1, pc.x2, xlabel = "PC 1", ylabel = "PC 2", c = c,
                 title = "PCA projection")
    plot(p0, p1, p2, layout = (1, 3), legend = false, size = (700, 400))
end

# ╔═╡ e4302d86-5e4f-4c56-bdca-8b4eed2af47c
tsne_proj = tsne(Array(eight_clouds), 2, 0, 2000, 50.0);

# ╔═╡ 6b319dbe-c626-4294-92b6-cab574aabfb0
let
    gr()
    scatter(tsne_proj[:, 1], tsne_proj[:, 2], legend = false,
            c = vcat([fill(i, 50) for i in 1:8]...))
end

# ╔═╡ 716da1a0-266f-41ec-8ef3-a2593ac7b74b
Markdown.parse("Also for real data PCA sometimes fails to find the clusters in
the data. We demonstrate this here with the MNIST data set.

```julia
mnist_x, mnist_y = let df = OpenML.load(554) |> DataFrame
    coerce!(df, :class => Multiclass)
    coerce!(df, Count => MLJ.Continuous)
    select(df, Not(:class)),
    df.class
end

mnist_pca = fit!(machine(PCA(maxoutdim = 2), mnist_x));

let
    mnist_proj = MLJ.transform(mnist_pca)
    scatter(mnist_proj.x1, mnist_proj.x2, c = Int.(int(mnist_y)), legend = false,
            xlabel = \"PC 1\", ylabel = \"PC 2\")
end
```
$(MLCourse.embed_figure("mnist_pca.png"))

The dots are colored according to their MNIST class label.

Looking at this it may appear puzzling that K-Means found somewhat reasonable
clusters in the MNIST dataset.

However, let us look at a TSne dimensionality reduction:

```julia
tsne_proj_mnist = tsne(Array(mnist_x[1:5000, :]), 2, 50, 1000, 20)
scatter(tsne_proj_mnist[:, 1], tsne_proj_mnist[:, 2],
        c = Int.(int(mnist_y[1:5000])),
        xlabel = \"TSne 1\", ylabel = \"TSne 2\",
        legend = false)
```
$(MLCourse.embed_figure("mnist_tsne.png"))

Here the images with different class labels are more clearly grouped into different clusters.
")

# ╔═╡ 3087540f-19e1-49cb-9b3b-c23207775ea7
md"# Principal Component Regression"

# ╔═╡ 60c37ce2-50d7-426f-9dd8-597878425f9c
pcr_data_x, pcr_data_y = let p1 = 10, p2 = 50, n = 120
	hidden = randn(n, p1)
	β = randn(p1)
	y = hidden * β .+ .1 * randn(n)
	transformation = randn(p1, p2)
	x = hidden * transformation .+ .1 * randn(n, p2)
	DataFrame(x, :auto), y
end;

# ╔═╡ 42aa3526-bd02-41fa-aae5-b5a8f3355b86
begin
	train_idx = rand((true, false), 120)
	test_idx = (!).(train_idx)
end;

# ╔═╡ 422026f0-8c87-471a-af51-ba6273f293b2
mach1 = fit!(machine(MLJLinearModels.LinearRegressor(),
		             pcr_data_x[train_idx, :],
		             pcr_data_y[train_idx]));

# ╔═╡ 3148b8bd-63fc-4052-a899-656d2c2057f8
rmse(predict(mach1, pcr_data_x[test_idx, :]), pcr_data_y[test_idx])

# ╔═╡ 434f1616-e028-4043-9aef-049e490525f5
mach2 = fit!(machine(@pipeline(PCA(), MLJLinearModels.LinearRegressor()),
		             pcr_data_x[train_idx, :],
		             pcr_data_y[train_idx]));

# ╔═╡ 9283a9d2-355b-439b-98b6-d8491904eb05
rmse(predict(mach2, pcr_data_x[test_idx, :]), pcr_data_y[test_idx])

# ╔═╡ eea3f7f8-83ab-4942-9017-504f340ad95b
md"""# Exercises

## Conceptual

**Exercise 1**

(a) The figure below shows 2-dimensional data points. Let us perform PCA on this dataset, i.e. estimate graphically the principal component loadings ``\phi_{11}, \phi_{21}, \phi_{12}, \phi_{22}``. *Hint*: start with determining graphically the direction of the first and the second principal component and normalize the direction vectors afterwards.

$(Random.seed!(818); let x = ([6 4] .* randn(400, 2)) * [-1 1
                                      -2 -.01]
    scatter(x[:, 1], x[:, 2], xlim = (-32, 32), ylim = (-16, 16),
            aspect_ratio = 1, legend = false, xlabel = "X₁", ylabel = "X₂")
  end)

(b) In the figure below you see the scores of ``n = 40`` measuments and the loadings of a principal component analysis. How many features ``p`` were measured in this data set?

$(Random.seed!(17); let data = DataFrame(Array(DataFrame(X1 = 6*randn(40), X2 = 4 * randn(40), X3 = randn(40)*2, X4 = randn(40)*2)) * randn(4, 4), :auto)
    biplot(fit!(machine(PCA(), data)))
end
)

(c) Estimate from the figure in (b) the loadings of the first two principal components and the scores of data point 35.

**Exercise 2**
The figures show biplots of principal component analysis applied to 30 three-dimensional points.

$(Random.seed!(123); let d = DataFrame(randn(30, 2) * randn(2, 3) .+ randn(30, 3) * 1e-6, :auto)
    p = fit!(machine(PCA(pratio = 1, mean = 0), d))
    plot(biplot(p), biplot(p, pc = (1, 3)), layout = (1, 2), size = (700, 350))
end)

(a) In which direction does this dataset have the highest variance? Give your answer as a vector with coefficients determined approximately from the figure.

(b) Is the two-dimensional plane that is closest to the data parallel to the plane spanned by x1 and x3? Justify your answer with a few words.

(c) How large approximately is the proportion of variance explained by the third component?

**Exercise 3** In this exercise you will explore the connection between PCA and SVD.

(a) In the lecture we said that the first principal component is the eigenvector ``\phi_1`` with the largest eigenvalue ``\lambda_1`` of the matrix ``X^T X``. From linear algebra you know that a real, symmetric matrix ``A`` can be diagonalized such that ``A = W \Lambda W^T`` holds, where ``\Lambda`` is a diagonal matrix that contains the eigenvalues and ``W`` is an orthogonal matrix that satisfies ``W^T W = I``, where ``I`` is the identity matrix. The columns of ``W`` are the eigenvectors of ``A``. Let us assume we have found the diagonalization ``X^T X = W \Lambda W^T``. Multiply this equation from the right with ``W`` and show that the columns of ``W`` contain the eigenvectors of ``X^T X``.

(b) The singular value decomposition (SVD) of ``X`` is ``X = U\Sigma V^T`` where ``U`` and ``V`` are orthogonal matrices and ``\Sigma`` is a diagonal matrix. Use the singular value decomposition to express the right-hand-side of ``X^TX = ...`` in terms of ``U``, ``V`` and ``\Sigma``, show how ``U, V, \Sigma, W`` and ``\Lambda`` relate to each other and prove the statements on slide 11 of the lecture. *Hint*: the product of two diagonal matrices is again a diagonal matrix.

## Applied

**Exercise 1** Look at the following three artificial data sets.

```julia
data1 = DataFrame(randn(50, 2) * randn(2, 3), :auto)
data2 = DataFrame(randn(50, 2) * randn(2, 3) .+ 0.5*randn(50, 3), :auto)
data3 = DataFrame(randn(50, 2) * randn(2, 10), :auto)
```

(a) Plot the datasets `data1` and `data2` with `scatter3d`. *Hint:* Plots can be displayed with different plotting backends. If you run `plotly()` in a cell of your notebook you switch to the interactive but slow plotly backendend. If you want a non-interactive but fast backend, run `gr()`.

(b) Write down your expectations on the curve of proportion of variance
explained for the tree datasets.

(c) Perform PCA and plot the proportion of variance explained.

**Exercise 2** In this exercise we look at a food consumption dataset with the relative
consumption of certain food items in European countries. The
numbers represent the percentage of the population consuming that food type.
You can load the data with
```julia
using CSV, DataFrames
data = CSV.read(download(\"https://openmv.net/file/food-consumption.csv\"), DataFrame)
```

(a) Perform PCA on this dataset and produce a two biplots: one with PC1 versus
PC2 and another with PC1 versus PC3. *Hint*: before applying PCA you may have to deal with missing data (see notebook on transformations) and convert integers to floating point numbers with `coerce!(data, Count => Continuous)`.

(b) Which are the important variables in the first three components?

(c) What does it mean for a country to have a high score in the first component?
Pick one country that has a high score in the first component and verify that
this country does indeed have this interpretation.

**Exercise 3**
Take a random image e.g.
```julia
using FileIO, Plots
img = FileIO.load(download(\"http://kenia.pordescubrir.com/wp-content/uploads/2008/08/windsurf.jpg\"))
```
Transform it to grayscale with `Gray.(img)` and to a data frame with
`DataFrame(float.(real.(Gray.(img))), :auto)`.

The goal is to compress this image with PCA.
To do so, treat each row of this image as a single data point.

(a) Perform PCA on this data.

(b) Reconstruct the image with as many principal components needed to explain 99% of variance.

(c) How many numbers are needed to store the image in this compressed form?

(d) How big is the compression factor?

"""

# ╔═╡ b574a023-21e4-4ae2-a198-7c9962368713
MLCourse.list_notebooks(@__FILE__)

# ╔═╡ 147f9c9c-72dd-4003-90b1-33fd68d762ce
MLCourse.footer()

# ╔═╡ Cell order:
# ╠═8bd459cb-20bb-483e-a849-e18caae3beef
# ╟─31bcc3f0-88c1-4035-b0de-b4217623dc22
# ╠═e24cff15-a966-405a-b083-4b86be09c6f1
# ╟─9e60bd23-bd0e-4311-a059-a39a70143307
# ╠═a0bb0a63-67e2-414e-a3f0-b76875c1295d
# ╠═4439f803-e70f-4231-82e9-88bb574498f4
# ╟─e38cdfcf-362f-4e6e-81a2-e58334bb8d5f
# ╟─16843059-269e-4557-a4b9-5d7034f99ba8
# ╠═6e033717-5815-4ffe-b991-81b72b99df9f
# ╟─b25687fc-aa11-4843-9a89-31c4c9b5bc23
# ╟─14bfdb2d-fd71-4c62-8823-6e292a424bf5
# ╠═f6117412-f6d8-41d7-9524-dc25fc891deb
# ╟─9da3a52e-23ca-4ce4-8524-0abeb4be2306
# ╟─00e15b1c-a12d-4ac2-b3c6-6fece1abef31
# ╟─60ba84d8-4f48-4d8a-b695-11edbe91ccc8
# ╠═b256ff25-431a-4c47-8195-d958d3d28c56
# ╟─abb36d40-ad0e-442c-ab87-b15837a2541c
# ╟─e81b5dee-dffb-4a4c-b471-0974599fd87d
# ╠═48528bf7-8a40-4773-944c-54d0446a5cac
# ╠═9a5e0976-1c0e-45fc-855e-b0fbec801e8c
# ╟─b5d2986e-3723-4463-ab05-cc176aea60d8
# ╟─eb421c10-a023-4faa-bcfb-374b8ebb9c38
# ╟─c1262304-6846-4811-ae82-397863255415
# ╠═b2688c20-a578-445f-a5c4-335c85931862
# ╟─f9a506bb-6bdf-48aa-8f81-8fb6de078533
# ╟─884a7428-2dd4-4091-b889-3acfbfa45d5b
# ╟─2b73abbb-8cef-4da0-b4cf-b640fe0ad102
# ╠═927627e0-9614-4234-823d-eb2e13229784
# ╠═c7935f29-21f3-414a-b013-a1abd0f9f502
# ╟─45ee714c-9b73-492e-9422-69f27be79de2
# ╠═b52d72e3-64a3-4b7e-901b-3512b97aec23
# ╟─ecc77003-c138-4ea9-93a7-27df46aa1e89
# ╠═89406c26-7cba-4dd2-a6a3-8d572b440e01
# ╟─743477bb-6cf0-419e-814b-e774738e8d89
# ╟─18ad4980-2960-4ea9-8c9a-bf17cd3546ff
# ╟─4a8188a0-c809-4014-9fef-cf6b5f830b1b
# ╟─df9a5e8d-0692-44e4-8fd0-43e8643c4b7c
# ╟─d3e5ca43-f70e-4a11-a3c9-fc0755179e4c
# ╟─936570b9-43da-4c32-b15a-d595f690b1a6
# ╟─9032ad56-123c-4585-9b67-c94e5ecb3899
# ╟─a8188f31-204a-406f-96a2-86768918a2f9
# ╠═265594dd-ba59-408a-925a-c2a173343df7
# ╠═51720585-09c8-40bb-8ec2-c2f25804d731
# ╠═d45faa6b-706b-4990-80d9-32acc001b242
# ╟─6c8c86f0-02ec-4ce7-9d16-1b4865c344df
# ╟─e72c09e3-91cb-4455-9880-628428730efe
# ╠═9f6b1079-209f-4a24-96c8-81bdcf8bfc32
# ╟─7eac92b8-0639-46b8-966e-cbf44de7ee86
# ╟─ab17ef21-ea29-4d62-b14e-b87cffc0a1f9
# ╠═da99f084-522b-4969-a65c-31d1548d3073
# ╠═d88509c3-51ca-4ff9-9c1f-c3e13ec0e133
# ╟─615a823f-32df-469d-b2a3-e9d2c35e36a7
# ╠═ded0e6cd-e969-47d5-81cf-e2a967bc0a34
# ╟─9724b203-ae86-4cb6-afec-68e26c01ab83
# ╠═715a0e0c-5105-4fee-8b93-03a966a1bd8e
# ╟─95ba7abd-c425-46d4-8382-fedec7ed58f6
# ╠═d2bfdda4-4496-4daa-bcb5-51d9a7519823
# ╠═694ad552-3088-4872-bed5-89b4dfc67219
# ╠═c54b4211-36aa-4d7c-8e7a-c67cf5bdabdc
# ╟─b9383ef0-8499-4e56-a8d6-606ec128676c
# ╟─af03d989-b19f-44ac-a275-d8af57bbeaeb
# ╠═026a362d-3234-409b-b9cb-e6a2a50c6a0d
# ╟─e967a9cb-ad95-4bd2-8fbd-491dc6fe0475
# ╠═0cb282aa-81d4-408d-a919-0d28cdbb7bed
# ╠═c70b5212-a823-4ff8-937d-10919efc766c
# ╟─0742ea40-1895-48f4-b149-d472397e83d0
# ╟─7c020b97-384d-43fc-af1c-a4bf75a6f06c
# ╟─cddfed10-7c9c-47b0-83fc-b3322da522ed
# ╟─03f92562-7419-4d8e-baa2-67f5b021219e
# ╟─8e466825-0c84-41c7-93c6-e3cc81a60ef5
# ╟─218d4c11-45a1-4d92-8b48-d08251804638
# ╠═0943f000-a148-4470-ae5e-9d8d39cd8207
# ╠═2e86337e-9fb8-4dd9-ac96-c1177e7ebd16
# ╟─24179171-3992-4253-88a5-15ccdbb26e13
# ╟─8bc5f5b1-cda1-4018-a215-5a1900c2a9d2
# ╟─d44be4c6-a409-421a-8e68-39fa752db4c0
# ╠═d3b769d7-ec71-4496-926a-e838c3b50c2a
# ╟─67395b94-46bd-43e5-9e12-0bd921de8203
# ╠═1c3aa1be-c58a-4eb9-a192-ddfffa8a2c56
# ╠═e4302d86-5e4f-4c56-bdca-8b4eed2af47c
# ╟─6b319dbe-c626-4294-92b6-cab574aabfb0
# ╟─716da1a0-266f-41ec-8ef3-a2593ac7b74b
# ╟─3087540f-19e1-49cb-9b3b-c23207775ea7
# ╠═60c37ce2-50d7-426f-9dd8-597878425f9c
# ╠═42aa3526-bd02-41fa-aae5-b5a8f3355b86
# ╠═422026f0-8c87-471a-af51-ba6273f293b2
# ╠═3148b8bd-63fc-4052-a899-656d2c2057f8
# ╠═434f1616-e028-4043-9aef-049e490525f5
# ╠═9283a9d2-355b-439b-98b6-d8491904eb05
# ╟─eea3f7f8-83ab-4942-9017-504f340ad95b
# ╟─b97724e4-d7b0-4085-b88e-eb3c5bcbe441
# ╟─b574a023-21e4-4ae2-a198-7c9962368713
# ╟─147f9c9c-72dd-4003-90b1-33fd68d762ce

ext: jl

max_stars_repo_name: jbrea/MLCourse

max_issues_repo_name: davibarreira/MLCourse

max_forks_repo_name: davibarreira/MLCourse


text:
const TEST_UUID = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
const TEST_DEP = PackageSpec(; name="Test", uuid=TEST_UUID)

@blueprint struct ProjectTest <: MustacheBlueprint
    template::TemplateFile = TemplateFile(path="test", file="runtests.jl")
    project::Bool = false
end

blueprint_view(bp::ProjectTest, ctx::Context) = Dict(
    "PKG" => ctx.kwargs.name,
)

function compile(bp::ProjectTest, ctx::Context)
    invoke(compile, Tuple{MustacheBlueprint, Context}, bp, ctx)

    f = bp.project ? make_test_project : add_test_dependency
    f(ctx.kwargs.path)
end

function make_test_project(pkg_dir::AbstractString)
    with_project(joinpath(pkg_dir, "test")) do
        Pkg.add(TEST_DEP)
    end
end

function add_test_dependency(pkg_dir::AbstractString)
    path = joinpath(pkg_dir, "Project.toml")
    toml = TOML.parsefile(path)
    get!(Dict, toml, "extras")["Test"] = TEST_UUID
    get!(Dict, toml, "targets")["test"] = ["Test"]
    write_project_toml(path, toml)

    touch(joinpath(pkg_dir, "Manifest.toml"))
    with_project(Pkg.update, pkg_dir)
end

function with_project(f, path::AbstractString)
    proj = Base.active_project()
    try
        Pkg.activate(path)
        return f()
    finally
        isnothing(proj) ? Pkg.activate() : Pkg.activate(proj)
    end
end


text:
# Fig_1_Brusselator_eigenvalues_v3.jl
# ==================================== u.p.  23.8.21

using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, OrdinaryDiffEq
using StaticArrays
using LinearAlgebra

# parameters
# ----------
a = 9
bs = [9.8, 10.2]
dvecs = [
    [4.5, 3.0, 2.0, 1.7],
    reverse([0.1, 1.0, 1.7, 2.0]),
]
fig, axs = subplots(1,2)

for (j, b) in enumerate(bs)

dvec = dvecs[j]
ivec = [1 2 6 4]   # order of colors

qmin = 0.0
qmax = 2.0
nq = 301
qvec = range(qmin,qmax; length=nq)
Re_lambda_vec = zeros(nq)
Im_lambda_vec = zeros(nq)
Re_qvec = zeros(nq)
Im_qvec = zeros(nq)

ax = axs[j]

ax.axhline(0; linewidth=1.0,color="k" )

for id = 1:length(dvec)

    dd = dvec[id]
    jj = kk = 0
    for iq = 1:nq
        cc = qvec[iq] * qvec[iq]
        b11 = b - 1 - cc
        b12 = a
        b21 = -b
        b22 = - a - dd*cc
        traceB = b11+b22
        detB = b11*b22 - b12*b21
        rootarg = traceB*traceB/4 - detB
        if rootarg < 0
            Re_lambda = traceB/2
            jj = jj + 1
            Im_qvec[jj] = qvec[iq]
            Im_lambda_vec[jj] = Re_lambda
        else
            Re_lambda = traceB/2 + sqrt(rootarg)
            if kk == 0
               kk = 1
               Re_qvec[kk] = Im_qvec[jj]
               Re_lambda_vec[kk] = Im_lambda_vec[jj]
            end
            kk = kk + 1
            Re_qvec[kk] = qvec[iq]
            Re_lambda_vec[kk] = Re_lambda
        end
    end
    ax.plot(Re_qvec[1:kk],Re_lambda_vec[1:kk], linewidth=3.0,linestyle="--",color=COLORS[ivec[id]] )
    ax.plot(Im_qvec[1:jj],Im_lambda_vec[1:jj], linewidth=3.0,linestyle="-",color=COLORS[ivec[id]] )
end

ax.set_xlim([qmin,qmax])
if j == 1
    ax.set_ylim([-2.1,2.6])
else
    ax.set_ylim([-2.1,0.4])
end
ax.legend([ "_" ; "_" ; L"d="*string(dvec[1]);
                  "_" ; L"d="*string(dvec[2]);
                  "_" ; L"d="*string(dvec[3]);
                  "_" ; L"d="*string(dvec[4]);
              #    "_" ; L"d="*string(dvec[5]);
                  ],
    fontsize = 26, handlelength = 1, handletextpad = 0.6,
    # loc = "upper left",
)

end

axs[1].set_xlabel(L"q")
axs[2].set_xlabel(L"q")
axs[1].set_ylabel(L"\max (\rm{Re} (\lambda))")

add_identifiers!(fig)
fig.tight_layout(pad=0.3)
wsave(plotsdir("11", "brusselator_eigenvalues"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
function get_units(unit_str)
    isempty(unit_str) && return 1
    # workaround for number density units
    unit_str == "#" && return 1
    # workaround stuff like kg.m/s
    unit_str = replace(unit_str, "."=>"*")
    # workaround stuff like 1/m^3
    if occursin(r"1\/([a-z]*\^[0-9]?)", unit_str)
        unit_str = replace(unit_str, r"1\/([a-z]*)\^([0-9]?)" => s"\g<1>^-\g<2>")
    end
    uparse(unit_str)
end

get_units(unit_str::NTuple) = get_units.(unit_str)


text:
using Test
using Literate
using Logging

with_logger(NullLogger()) do
    Literate.script("$(@__DIR__)/../docs/src/usage.jl", "$(@__DIR__)/build/")
end

@testset "Readme example usage" begin

    include("build/usage.jl")

end

max_issues_repo_name: Juice-jl/LogicCircuits

max_forks_repo_name: Juice-jl/LogicCircuits


text:
sevendigit = [0 1 1 1 1 1 0;
              0 1 1 1 1 1 0;
              0 0 0 0 1 1 0;
              0 0 0 1 1 0 0;
              0 0 0 1 1 0 0;
              0 0 1 1 0 0 0;
              0 0 1 1 0 0 0;
              0 1 1 0 0 0 0;
              0 1 1 0 0 0 0]

max_stars_repo_name: seanrattana/courses

max_issues_repo_name: seanrattana/courses

max_forks_repo_name: seanrattana/courses


text:
export UnitDisk, rand_unitdisk

struct UnitDisk{T} <: AbstractSites
    locs::Vector{Tuple{T,T}}
    unit::T
end

unit(ud::UnitDisk) = ud.unit

vertices(ud::UnitDisk) = 1:length(ud.locs)
Base.getindex(ud::UnitDisk, i::Int) = ud.locs[i]

isconnected(ud::UnitDisk, i, j) = distance(ud[i], ud[j]) <= ud.unit

function rand_unitdisk(n::Int, ρ::Real; ndims::Int=2)
    unit = (ρ/n)^(1/ndims)
    UnitDisk([(rand(ndims)...,) for i=1:n], unit)
end


text:
module DiffEqOperators

using Base: Number
import Base: +, -, *, /, \, size, getindex, setindex!, Matrix, convert, ==
using DiffEqBase, StaticArrays, LinearAlgebra
import LinearAlgebra: mul!, ldiv!, lmul!, rmul!, axpy!, opnorm, factorize, I
import DiffEqBase: update_coefficients!, isconstant
using SciMLBase: AbstractDiffEqLinearOperator, AbstractDiffEqCompositeOperator, DiffEqScaledOperator
import SciMLBase: getops
using SparseDiffTools
using SparseArrays, ForwardDiff, BandedMatrices, NNlib, LazyArrays, BlockBandedMatrices, LoopVectorization
using LazyBandedMatrices
using RuntimeGeneratedFunctions
using Requires
RuntimeGeneratedFunctions.init(@__MODULE__)

abstract type AbstractDiffEqAffineOperator{T} end
abstract type AbstractDerivativeOperator{T} <: AbstractDiffEqLinearOperator{T} end
abstract type AbstractMatrixFreeOperator{T} <: AbstractDiffEqLinearOperator{T} end

### Matrix-free Operators
include("matrixfree_operators.jl")
include("jacvec_operators.jl")

### Utilities
include("utils.jl")

### Exceptions
include("exceptions.jl")

### Boundary Padded Arrays
include("boundary_padded_arrays.jl")

### Boundary Operators
include("derivative_operators/bc_operators.jl")
include("derivative_operators/multi_dim_bc_operators.jl")

### Derivative Operators
include("derivative_operators/fornberg.jl")
include("derivative_operators/derivative_operator.jl")
include("derivative_operators/abstract_operator_functions.jl")
include("derivative_operators/convolutions.jl")
include("derivative_operators/ghost_derivative_operator.jl")
include("derivative_operators/derivative_operator_functions.jl")
include("derivative_operators/coefficient_functions.jl")

### Vector Calculus Operators
include("derivative_operators/vector_calculus_operators.jl")
include("derivative_operators/vector_calculus_convolutions.jl")
include("derivative_operators/vector_algebraic_operations.jl")

### Composite Operators
include("composite_operators.jl")
include("docstrings.jl")

### Concretizations
include("derivative_operators/concretization.jl")

@warn("DiffEqOperators.jl is in the process of being deprecated.\n
       - For automated PDE discretization, see MethodOfLines.jl\n
       - For MatrixFreeOperators, and other non-derivative operators, see SciMLOperators.jl\n
       - For VecJacOperators and JacVecOperators, see SparseDiffTools.jl")

# The (u,p,t) and (du,u,p,t) interface
for T in [DiffEqScaledOperator, DiffEqOperatorCombination, DiffEqOperatorComposition, GhostDerivativeOperator]
  (L::T)(u,p,t) = (update_coefficients!(L,u,p,t); L * u)
  (L::T)(du,u,p,t) = (update_coefficients!(L,u,p,t); mul!(du,L,u))
end

function __init__()
  @require Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" begin
    export VecJacOperator

    include("vecjac_operators.jl")
  end
end

export MatrixFreeOperator
export AnalyticalJacVecOperator, JacVecOperator, getops
export AbstractDerivativeOperator, DerivativeOperator,
       CenteredDifference, UpwindDifference, CompleteCenteredDifference, CompleteUpwindDifference, CompleteHalfCenteredDifference, nonlinear_diffusion, nonlinear_diffusion!,
       GradientOperator, Gradient, CurlOperator, Curl, DivergenceOperator, Divergence
export DirichletBC, Dirichlet0BC, NeumannBC, Neumann0BC, RobinBC, GeneralBC, MultiDimBC, PeriodicBC,
       MultiDimDirectionalBC, ComposedMultiDimBC
export compose, decompose, perpsize, square_norm, square_norm!, dot_product, dot_product!, cross_product,
       cross_product!

export GhostDerivativeOperator
export BoundaryConditionError
end # module


text:
mutable struct Basis{O, V, P} <: abstractBasis
    basis::O
    variables::V
    parameter::P
    f_
end

is_independent(o::Operation) = isempty(o.args)

Base.print(io::IO, x::Basis) = show(io, x)
Base.show(io::IO, x::Basis) = print(io, "$(length(x.basis)) dimensional basis in ", "$(String.([v.op.name for v in x.variables]))")

function Basis(basis::AbstractVector{Operation}, variables::AbstractVector{Operation};  parameters =  [])
    @assert all(is_independent.(variables)) "Please provide independent variables for base."

    bs = unique(basis)

    vs = sort!([b for b in [ModelingToolkit.vars(bs)...] if !b.known], by = x -> x.name)
    ps = sort!([b for b in [ModelingToolkit.vars(bs)...] if b.known], by = x -> x.name )

    f_ = ModelingToolkit.build_function(bs, vs, ps, (), simplified_expr, Val{false})[1]
    return Basis(bs, variables, parameters, f_)
end

function update!(b::Basis)
    vs = sort!([bi for bi in [ModelingToolkit.vars(b.basis)...] if !bi.known], by = x->x.name)
    ps = sort!([bi for bi in [ModelingToolkit.vars(b.basis)...] if bi.known], by = x->x.name)

    b.f_ = ModelingToolkit.build_function(b.basis, vs, ps, (), simplified_expr, Val{false})[1]
    return
end

function Base.push!(b::Basis, op::Operation)
    push!(b.basis, op)
    # Check for uniqueness
    unique!(b)
    update!(b)
    return
end

function Base.deleteat!(b::Basis, inds)
    deleteat!(b.basis, inds)
    update!(b)
    return
end

(b::Basis)(x::AbstractArray; p::AbstractArray = []) = b.f_(x, p)

Base.size(b::Basis) = size(b.basis)
ModelingToolkit.parameters(b::Basis) = b.parameter
variables(b::Basis) = b.variables

function jacobian(b::Basis)
    vs = sort!([bi for bi in [ModelingToolkit.vars(b.basis)...] if !bi.known], by = x-> x.name)
    ps = sort!([bi for bi in [ModelingToolkit.vars(b.basis)...] if bi.known], by = x-> x.name)
    j = calculate_jacobian(b.basis, variables(b))
    return ModelingToolkit.build_function(expand_derivatives.(j), vs, ps, (), simplified_expr, Val{false})[1]
end

function Base.unique!(b::Basis)
    N = length(b.basis)
    removes = Vector{Bool}()
    for i ∈ 1:N
        push!(removes, any([isequal(b.basis[i], b.basis[j]) for j in i+1:N]))
    end
    deleteat!(b, removes)
end

function Base.unique(b::Basis)
    N = length(b.basis)
    returns = Vector{Bool}()
    for i ∈ 1:N
        push!(returns, any([isequal(b.basis[i], b.basis[j]) for j in i+1:N]))
    end
    returns = [!r for r in returns]
    return Basis(b.basis[returns], variables(b), parameters = parameters(b))
end

function dynamics(b::Basis)
    return b.f_
end

function ModelingToolkit.ODESystem(b::Basis)
    # Define the time
    @parameters t
    @derivatives D'~t

    vs = similar(b.variables)
    dvs = similar(b.variables)
    for (i, vi) in enumerate(b.variables)
        vs[i] = ModelingToolkit.Operation(vi.op, [t])
        dvs[i] = D(vs[i])
    end
    eqs = dvs .~ b(vs, p = b.parameter)
    return ODESystem(eqs)
end


text:
push!(LOAD_PATH, pwd())
using NQCDynamics
using NQCDynamics.InitialConditions
using LinearAlgebra
using Unitful
using UnitfulAtomic
using Revise
using Random, Distributions
using Plots
using DelimitedFiles
#using Profile
# using DataFrames
#using BenchmarkTools

# IESH should proceed by the following steps:
#   1) initialize classical positions and momenta, initialize
#      electronic states, define surface on which to propagate
#   2) Do Verlet algorithm for the atomic positions, velocities electronic
#   3) Integrate the electronic wave function w/ 1e- Schroedinger equationa
#   4) calculate switching probability, calculate if hop should occur
#   5) If no hop, GOTO 2; if hop, see if enough energy to do it, GOTO 2

ps = 98.22694788464062
Random.seed!(13)
#n_states = 202
n_states = 42
#n_states = 502
# # number of electronic n_states
vinit = 0.005 # initial velocity
mass = 1.5 # atomic mass#
B_na = 1.0 # parameter that defines nonadiabatic coupling strength
n_DOF = 1
ntrajes = 20
tspan = (0.0, 5.0u"ps")

# 1a) Define adiabatic model, fom Model/scattering_anderson_holsteins.jl
# model = Models.ScatteringAndersonHolstein(N=n_states, a0=1, Bg=B_na,
# Cg=1, alpha=-5.0, beta=.5)
#model = Models.MiaoSubotnik(n_states=n_states, W=0.064, Gamma=0.0064)
model = Models.MiaoSubotnik(n_states=n_states, W=0.064, Gamma=0.0064)
#model = Models.TullyModelOne()
#println(model.n_states)

# 1b) Initialize atomic parameters, i.e. moving atoms, of which there is one
atoms = Atoms{Float64}([:H])

##############################################################################
#
# Try dynamics
#
###############################################################################

# create Boltzmann distribution of momenta
# kT = 9.5*10^-4 a.u. = 300 K
temperature = 5*298u"K"

# create distribution of atomic positions from Monte Carlo sampling of the PESs
# for the number of trajectories and atoms
# atoms and step size
Δ = Dict([(:H, 0.3)])
tsim = Simulation{Classical}(atoms, model; temperature=temperature)
R0 = zeros(3, length(tsim.atoms))
Δ = Dict([(:H, 0.3)])
R0 = zeros(1, length(tsim.atoms))
passes = 1000
outbolz = MetropolisHastings.run_monte_carlo_sampling(tsim, R0, Δ, passes)

# If sampled, gives an array of momenta and positions from the above initialized sampling
#bolz_pos_momenta = DynamicalDistribution(boltzmann, outbolz.R, (1, 1))

# # Initialize the simulation problem; Simulation is defined in src/simulation_constructors.jl
sim = Simulation{wave_IESH}(atoms, model)

# # Do dynamics
for i=1:ntrajes
    #i = 1
     nname=lpad(i,8,"0")
#     r, v = rand(bolz_pos_momenta)
#     v = -v
#     println(i," ", r, " ", v)

    #r = fill(-5.0, sim.DoFs, length(sim.atoms))
    r = rand(outbolz.R)
    v = fill(sqrt(rand(outbolz.energy)*2/sim.atoms.masses[1]), size(sim))
    #v = fill(5/sim.atoms.masses[1], sim.DoFs, length(sim.atoms))
    println(r, v)

    # intial state (how to initialize a vector that shows which states are occupied)
    # For now, just initialize in ground state. Makes sense, anyhow, when the molecule
    # is hurling towards the surface, inititially
    # This is defined in src/Dynamics/SurfaceHopping/surface_hopping_variables.jl
    # as: SurfaceHoppingVariables(v::AbstractArray, r::AbstractArray, n_states::Integer, state::Integer)
    k = round.(Int64, (zeros(n_states)))
    l = Int(n_states / 2)
    k[1:l] .= 1
    println(k)
    z = SurfaceHoppingVariablesIESH(v, r, n_states, k)

    #run_trajectory defined in src/Dynamics.jl
    # callbacks defined in: src/Dynamics/callbacks.jl and src/Models/Models.jl
    # Save impurity does not seem to be defined at the moment?
    #@time solution = Dynamics.run_trajectory(z, (0.0, 10000000.0), sim; output=(:save_impurity))
    @time solution = Dynamics.run_trajectory(z, (0.0, 100000.0), sim, dt=10, adaptive=false;
                                             output=(:position, :save_impurity))
    println("Finished")
    println(length(solution.t))

    aka = Int(ceil(length(solution)/10))-1
    j = 0
    #aka = Int(length(solution)
    println(aka)
    outarray = zeros(aka, Int(5 + model.n_states))
    open("trajectory_$nname.txt", "w") do fi
        write(fi, "step, r (a.u.), epot (a.u.), state, impurity population\n")
        for k = 1:length(solution)
            if (mod(k,10)==0)
                j = j + 1
                outarray[j,1] = solution.t[k]
                outarray[j,2] = solution.position[k][1]
                outarray[j,3] = solution.save_impurity[k][2]
                outarray[j,4] = solution.save_impurity[k][3]
                outarray[j,5] = solution.save_impurity[k][4]
                for imp in 1:model.n_states
                    outarray[j,5+imp] = solution.save_impurity[k][4+imp]
                end

            # outarray[i,3] = solution.save_impurity[i][2]
            # outarray[i,4] = solution.save_impurity[i][3]
            # outarray    [i,5] = solution.save_impurity[i][4]
            end
        end
        writedlm(fi, outarray,'\t')
    end

# b = plot(outarray[:,1], outarray[:,2], label="energy", marker=2)
# xlabel!("t")
# ylabel!("x")
a = plot(outarray[:,2], outarray[:,3], label="energy", marker=2)
xlabel!("x")
ylabel!("Energy (a.u.)")
b = plot(outarray[:,1], outarray[:,Int(5+42/2)], label="energy", marker=2)
xlabel!("time")
ylabel!("Impurity Population")

display(plot(b))
display(plot(a))

end
# # # savefig(a,"traj_subA_G4Em4_W10G_test5.png")
# # # savefig(b,"traj_imppop_G4Em4_W10G_test2.png")

#display(plot(b))


text:
function vertexbuffer(x, trans)
    pos = decompose(Point, x)
    return apply_transform(trans,  pos)
end

function vertexbuffer(x::Observable, p)
    return Buffer(lift(vertexbuffer, x, transform_func_obs(p)))
end

facebuffer(x) = facebuffer(GeometryBasics.faces(x))
facebuffer(x::Observable) = Buffer(lift(facebuffer, x))
function facebuffer(x::AbstractArray{GLTriangleFace})
    return x
end

function array2color(colors, cmap, crange)
    cmap = RGBAf0.(Colors.color.(to_colormap(cmap)), 1.0)
    return AbstractPlotting.interpolated_getindex.((cmap,), colors, (crange,))
end

function array2color(colors::AbstractArray{<:Colorant}, cmap, crange)
    return RGBAf0.(colors)
end

function converted_attribute(plot::AbstractPlot, key::Symbol)
    return lift(plot[key]) do value
        return convert_attribute(value, Key{key}(), Key{plotkey(plot)}())
    end
end

function create_shader(scene::Scene, plot::AbstractPlotting.Mesh)
    # Potentially per instance attributes
    mesh_signal = plot[1]
    mattributes = GeometryBasics.attributes
    get_attribute(mesh, key) = lift(x -> getproperty(x, key), mesh)
    data = mattributes(mesh_signal[])

    uniforms = Dict{Symbol,Any}()
    attributes = Dict{Symbol,Any}()

    for (key, default) in (:uv => Vec2f0(0), :normals => Vec3f0(0))
        if haskey(data, key)
            attributes[key] = Buffer(get_attribute(mesh_signal, key))
        else
            uniforms[key] = Observable(default)
        end
    end

    if haskey(data, :attributes) && data[:attributes] isa AbstractVector
        attr = get_attribute(mesh_signal, :attributes)
        attr_id = get_attribute(mesh_signal, :attribute_id)
        color = lift((c, id) -> c[Int.(id) .+ 1], attr_id)
        attributes[:color] = Buffer(color)
        uniforms[:uniform_color] = false
    else
        color_signal = converted_attribute(plot, :color)
        color = color_signal[]
        mesh_color = color_signal[]
        uniforms[:uniform_color] = Observable(false) # this is the default

        if color isa Colorant && haskey(data, :color)
            color_signal = get_attribute(mesh_signal, :color)
            color = color_signal[]
        end

        if color isa AbstractArray
            c_converted = if color isa AbstractArray{<:Colorant}
                color_signal
            elseif color isa AbstractArray{<:Number}
                lift(array2color, color_signal, plot.colormap, plot.colorrange)
            else
                error("Unsupported color type: $(typeof(color))")
            end
            if c_converted[] isa AbstractVector
                attributes[:color] = Buffer(c_converted) # per vertex colors
            else
                uniforms[:uniform_color] = Sampler(c_converted) # Texture
                !haskey(attributes, :uv) &&
                    @warn "Mesh doesn't use Texturecoordinates, but has a Texture. Colors won't map"
            end
        elseif color isa Colorant && !haskey(attributes, :color)
            uniforms[:uniform_color] = color_signal
        else
            error("Unsupported color type: $(typeof(color))")
        end
    end

    if !haskey(attributes, :color)
        uniforms[:color] = Vec4f0(0) # make sure we have a color attribute
    end

    uniforms[:shading] = plot.shading

    for key in (:ambient, :diffuse, :specular, :shininess)
        uniforms[key] = plot[key]
    end

    faces = facebuffer(mesh_signal)
    positions = vertexbuffer(mesh_signal, plot)
    instance = GeometryBasics.Mesh(GeometryBasics.meta(positions; attributes...), faces)

    get!(uniforms, :colorrange, true)
    get!(uniforms, :colormap, true)
    get!(uniforms, :model, plot.model)
    get!(uniforms, :lightposition, Vec3f0(1))

    uniforms[:normalmatrix] = map(scene.camera.view, plot.model) do v, m
        i = SOneTo(3)
        return transpose(inv(v[i, i] * m[i, i]))
    end

    return Program(WebGL(), lasset("mesh.vert"), lasset("mesh.frag"), instance; uniforms...)
end


text:
import Clustering
import Distances
import LinearAlgebra: norm

using MLJBase
using MLJClusteringInterface
using Random:seed!
using Test

const Dist = Distances

seed!(132442)
X, y = @load_crabs

####
#### KMEANS
####

@testset "Kmeans" begin
    barekm = KMeans()
    fitresult, cache, report = fit(barekm, 1, X)
    R = matrix(transform(barekm, fitresult, X))
    X_array = matrix(X)
    # distance from first point to second center
    @test R[1, 2] ≈ norm(view(X_array, 1, :) .- view(fitresult[1], :, 2))^2
    @test R[10, 3] ≈ norm(view(X_array, 10, :) .- view(fitresult[1], :, 3))^2

    p = predict(barekm, fitresult, X)
    @test argmin(R[1, :]) == p[1]
    @test argmin(R[10, :]) == p[10]

end

####
#### KMEDOIDS
####

@testset "Kmedoids" begin
    barekm = KMedoids()
    fitresult, cache, report = fit(barekm, 1, X)
    X_array = matrix(X)
    R = matrix(transform(barekm, fitresult, X))
    @test R[1, 2] ≈ Dist.evaluate(
        barekm.metric, view(X_array, 1, :), view(fitresult[1], :, 2)
    )
    @test R[10, 3] ≈ Dist.evaluate(
        barekm.metric, view(X_array, 10, :), view(fitresult[1], :, 3)
    )
    p = predict(barekm, fitresult, X)
    @test all(report.assignments .== p)
end

####
#### HierarchicalClustering
####

@testset "HierarchicalClustering" begin
    barehc = HierarchicalClustering(linkage = :complete)
    hc = transform(barehc, X)
    @test mean(int(y) .== hc(k = 2)) > .5 # that's a somewhat silly test...
end


text:
#
# function eulermat: This routine was added because it defines
#                    the rotation in the "human" way, an is thus used
#                    to set the position of the fixed molecules.
#     That means: beta is a counterclockwise rotation around x axis.
#                 gamma is a counterclockwise rotation around y axis.
#                 theta is a counterclockwise rotation around z axis.
#

function eulermat( beta :: Float64,
                   gamma :: Float64,
                   theta :: Float64,
                   deg :: String )

   if deg != "degree"
     error("ERROR: to use radians just omit the last parameter")
   end

   beta = pi * beta / 180.
   gamma = pi * gamma / 180.
   theta = pi * theta / 180.
   A = eulermat( beta, gamma, theta )
   return A

end

function eulermat( beta :: Float64,
                   gamma :: Float64,
                   theta :: Float64 )

  A = Matrix{Float64}(undef,3,3)

  c1 = cos(beta)
  s1 = sin(beta)
  c2 = cos(gamma)
  s2 = sin(gamma)
  c3 = cos(theta)
  s3 = sin(theta)

  A[1,1] = c2*c3
  A[1,2] = c1*s3 + c3*s1*s2
  A[1,3] = s1*s3 - c1*c3*s2

  A[2,1] = -c2*s3
  A[2,2] = c1*c3 - s1*s2*s3
  A[2,3] = c1*s2*s3 + c3*s1

  A[3,1] = s2
  A[3,2] = -c2*s1
  A[3,3] = c1*c2

  return A

end

max_stars_repo_name: mcubeg/JPleXL

max_issues_repo_name: mcubeg/JPleXL

max_forks_repo_name: mcubeg/JPleXL


text:
## Diagonal matrices

immutable Diagonal{T} <: AbstractMatrix{T}
    diag::Vector{T}
end
Diagonal(A::Matrix) = Diagonal(diag(A))

convert{T}(::Type{Diagonal{T}}, D::Diagonal{T}) = D
convert{T}(::Type{Diagonal{T}}, D::Diagonal) = Diagonal{T}(convert(Vector{T}, D.diag))
convert{T}(::Type{AbstractMatrix{T}}, D::Diagonal) = convert(Diagonal{T}, D)
convert{T}(::Type{Triangular}, A::Diagonal{T}) = Triangular{T, Diagonal{T}, :U, false}(A)

function similar{T}(D::Diagonal, ::Type{T}, d::(Int,Int))
    d[1] == d[2] || throw(ArgumentError("Diagonal matrix must be square"))
    return Diagonal{T}(Array(T,d[1]))
end

copy!(D1::Diagonal, D2::Diagonal) = (copy!(D1.diag, D2.diag); D1)

size(D::Diagonal) = (length(D.diag),length(D.diag))
size(D::Diagonal,d::Integer) = d<1 ? error("dimension out of range") : (d<=2 ? length(D.diag) : 1)

fill!(D::Diagonal, x) = (fill!(D.diag, x); D)

full(D::Diagonal) = diagm(D.diag)
getindex(D::Diagonal, i::Integer, j::Integer) = i == j ? D.diag[i] : zero(eltype(D.diag))

function getindex(D::Diagonal, i::Integer)
    n = length(D.diag)
    id = div(i-1, n)
    id + id * n == i-1 && return D.diag[id+1]
    zero(eltype(D.diag))
end

ishermitian(D::Diagonal) = true
issym(D::Diagonal) = true
isposdef(D::Diagonal) = all(D.diag .> 0)

tril!(D::Diagonal,i::Integer) = i == 0 ? D : zeros(D)
triu!(D::Diagonal,i::Integer) = i == 0 ? D : zeros(D)

==(Da::Diagonal, Db::Diagonal) = Da.diag == Db.diag

+(Da::Diagonal, Db::Diagonal) = Diagonal(Da.diag + Db.diag)
-(Da::Diagonal, Db::Diagonal) = Diagonal(Da.diag - Db.diag)

*{T<:Number}(x::T, D::Diagonal) = Diagonal(x * D.diag)
*{T<:Number}(D::Diagonal, x::T) = Diagonal(D.diag * x)
/{T<:Number}(D::Diagonal, x::T) = Diagonal(D.diag / x)
*(Da::Diagonal, Db::Diagonal) = Diagonal(Da.diag .* Db.diag)
*(D::Diagonal, V::Vector) = D.diag .* V
*(A::Matrix, D::Diagonal) = scale(A,D.diag)
*(D::Diagonal, A::Matrix) = scale(D.diag,A)

A_mul_B!(A::Diagonal,B::AbstractMatrix) = scale!(A.diag,B)
At_mul_B!(A::Diagonal,B::AbstractMatrix)= scale!(A.diag,B)
Ac_mul_B!(A::Diagonal,B::AbstractMatrix)= scale!(conj(A.diag),B)

/(Da::Diagonal, Db::Diagonal) = Diagonal(Da.diag ./ Db.diag )
function A_ldiv_B!{T}(D::Diagonal{T}, v::AbstractVector{T})
    length(v)==length(D.diag) || throw(DimensionMismatch(""))
    for i=1:length(D.diag)
        d = D.diag[i]
        d==zero(T) && throw(SingularException(i))
        v[i] *= inv(d)
    end
    v
end
function A_ldiv_B!{T}(D::Diagonal{T}, V::AbstractMatrix{T})
    size(V,1)==length(D.diag) || throw(DimensionMismatch(""))
    for i=1:length(D.diag)
        d = D.diag[i]
        d==zero(T) && throw(SingularException(i))
        V[i,:] *= inv(d)
    end
    V
end

conj(D::Diagonal) = Diagonal(conj(D.diag))
transpose(D::Diagonal) = D
ctranspose(D::Diagonal) = conj(D)

diag(D::Diagonal) = D.diag
trace(D::Diagonal) = sum(D.diag)
det(D::Diagonal) = prod(D.diag)
logdet{T<:Real}(D::Diagonal{T}) = sum(log(D.diag))
function logdet{T<:Complex}(D::Diagonal{T}) #Make sure branch cut is correct
    x = sum(log(D.diag))
    -pi<imag(x)<pi ? x : real(x)+(mod2pi(imag(x)+pi)-pi)*im
end
# identity matrices via eye(Diagonal{type},n)
eye{T}(::Type{Diagonal{T}}, n::Int) = Diagonal(ones(T,n))

expm(D::Diagonal) = Diagonal(exp(D.diag))
sqrtm(D::Diagonal) = Diagonal(sqrt(D.diag))

#Linear solver
function \{TD<:Number,TA<:Number}(D::Diagonal{TD}, A::AbstractArray{TA,1})
    m, n = size(A,2)==1 ? (size(A,1),1) : size(A)
    m==length(D.diag) || throw(DimensionMismatch(""))
    (m == 0 || n == 0) && return A
    C = Array(typeof(one(TD)/one(TA)),size(A))
    for j = 1:n
        for i = 1:m
            di = D.diag[i]
            di==0 && throw(SingularException(i))
            C[i,j] = A[i,j] / di
        end
    end
    return C
end
\(Da::Diagonal, Db::Diagonal) = Diagonal(Db.diag ./ Da.diag)

function inv{T}(D::Diagonal{T})
    Di = similar(D.diag)
    for i = 1:length(D.diag)
        D.diag[i]==zero(T) && throw(SingularException(i))
        Di[i]=inv(D.diag[i])
    end
    Diagonal(Di)
end

#Eigensystem
eigvals{T<:Number}(D::Diagonal{T}) = D.diag
eigvals(D::Diagonal) = [eigvals(x) for x in D.diag] #For block matrices, etc.
eigvecs(D::Diagonal) = eye(D)
eigfact(D::Diagonal) = Eigen(eigvals(D), eigvecs(D))

#Singular system
svdvals(D::Diagonal) = sort(D.diag, rev = true)
function svdfact(D::Diagonal, thin=true)
    S = abs(D.diag)
    piv = sortperm(S, rev=true)
    U = full(Diagonal(D.diag./S))
    Up= hcat([U[:,i] for i=1:length(D.diag)][piv]...)
    V = eye(D)
    Vp= hcat([V[:,i] for i=1:length(D.diag)][piv]...)
    SVD(Up, S[piv], Vp')
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
# Autogenerated wrapper script for duf_jll for i686-linux-gnu
export duf

JLLWrappers.@generate_wrapper_header("duf")
JLLWrappers.@declare_executable_product(duf)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        duf,
        "bin/duf",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Photosynthesis
using PkgUtility
using Test

include("test_Photosynthesis.jl")


text:
#@debug
function TRARC(nlp :: AbstractNLPModel,
               x₀ :: Array{Float64,1},
               TR :: TrustRegion, c :: Combi;
               s :: TStopping = TStopping(nlp),
               verbose :: Bool = true
               )

    hessian_rep,PData,solve_model,pre_process,decrease,params = extract(c)

    StopNorm = s.optimality_residual

    s = start!(s,x₀)

    α = TR.α₀
    x, xnext, d, Df = copy(x₀), copy(x₀), copy(x₀), 0.0
    xopt = x
    λ = 1.0

    n = length(x)
    ∇f = Array(Float64,n)
    ∇fnext = Array(Float64,n)

    f = obj(nlp,x)
    if isnan(f) | (f==Inf)
        OK = false
	println("f nan or ∞");
	xopt=f
	∇fopt=f
	fopt=f
        niter = Inf
        calls = [Inf, Inf, Inf, Inf]
	return xopt,∇fopt,∇fopt,fopt,niter,calls,OK,:Bad_x0
    end

    fopt = f
    grad!(nlp,x,∇f)
    norm_∇f = stop_norm(∇f)
    norm_∇f0 = norm_∇f
    ∇fopt = ∇f
    norm_∇fopt = norm_∇f
    H = hessian_rep(nlp,x)

    #optimal = (norm_g < atol) || (isinf(f) & (f<0.0))
    #tired = false

    fnext = f
    iter = 0

    optimal, unbounded, tired, elapsed_time = stop(s,iter,x,f,∇f)
    stalled = false
    finish = optimal || unbounded || tired || stalled

    verbose && display_header_iterations()
    verbose && display_success(iter, fnext, norm_g0, 0.0, α)

    succ, unsucc, verysucc, unsuccinarow = 0, 0, 0, 0

    max_unsuccinarow = TR.max_unsuccinarow

    calls = [0, 0, 0, 0]

    while ~finish

        PData = pre_process(H,∇f,params,calls,s.max_eval)
#                println("===>>> cond(H) = $(cond(H))  cond(L) = $(cond(PData.L))")
#                println(" ||H - P'LDL'P|| = $(norm(H - PData.P'*PData.L*PData.D*PData.L'*PData.P))")
        #if cond(PData.L) > 10000.0*cond(H) println("===>>> cond(H) = $(cond(H))  cond(L) = $(cond(PData.L))")
        #end

        if ~PData.OK return x,fopt,norm_∇f,Inf,Inf,[Inf,Inf,Inf,Inf],false,:PreFailed;  end

        success = false

        while ~success & (iter < s.max_iter) & (unsuccinarow < TR.max_unsuccinarow)
            try
                d,λ = solve_model(PData,α)
            catch
                println(" Problem in solve_model")

                return x,fopt,norm_∇f,norm_∇fopt,iter,calls,false,:AscentDir
            end
            #println("******* TRARC:  g⋅d = $(g⋅d), 0.5 d'*H*d = $(0.5*(H * d)⋅d)")

            Δq = -(∇f + 0.5*H*d)⋅d

            if Δq<0.0 println("*******   Ascent direction in SolveModel: Δq = $Δq")
                println("  g⋅d = $(∇f⋅d), 0.5 d'Hd = $(0.5*(H*d)⋅d)  α = $α  λ = $λ")
#@bp
                #try println(" cond(H) = $(cond(full(H)))") catch println("sparse hessian, no cond") end
                #D,Q=eig(full(H))
                #lm=findmin(D);lM=findmax(D)
                #println("λ min (H) = $lm  λ max (H) = $lM")
                #try printtln(" cond(L) = $(cond(PData.L))") catch println("Spectral") end
                #try println(" ||H - reconstructH|| = $(norm(H-reconstructH(PData)))")
                #catch
                #    println(" reconstruct H no available ")
                #end
                println("  calls  $calls   iters $verysucc; $succ; $unsucc")

                return x,fopt,norm_∇f,norm_∇fopt,iter,calls,false,:AscentDir
            end
            slope = ∇f ⋅ d
            xnext = x + d

            iter += 1

            # Trap illegal values for f, descent will backtrack until a legal value is reached
            try
                fnext = obj(nlp,xnext)
            catch
                fnext = Inf;
            end
            if isnan(fnext)
                fnext=Inf
            end

            Δf = f - fnext

            r, good_grad, ∇fnext = compute_r(nlp,f,Δf,Δq,slope,d,xnext,∇fnext)

            if r<TR.acceptance_threshold
                verbose && display_failure(iter,fnext,λ,α)
	        unsucc=unsucc+1
	        unsuccinarow = unsuccinarow +1
	        α = decrease(PData, α, TR)
                fbidon = obj(nlp,x)
	    else
	        success = true

	        unsuccinarow = 0
	        x = copy(xnext)
	        f = fnext

	        if good_grad  ∇f = copy(∇fnext)
                else grad!(nlp,x,∇f);
                end

                norm_∇f=stop_norm(∇f)
                H = hessian_rep(nlp,x)
	        if r > TR.increase_threshold
	            α = increase(PData, α, TR)
	            verbose && display_v_success(iter,f,norm_∇f,λ,α)
                    verysucc += 1
	        else
                    if r < TR.reduce_threshold
                        α = decrease(PData, α, TR)
                    end
	            verbose && display_success(iter,f,norm_∇f,λ,α)
	            succ += 1
	        end
	    end
        end
        if isnan(f)  OK=false
	    println("f nan");
	    xopt=f
	    ∇fopt=f
	    fopt=f
            niter = Inf
            calls = [Inf, Inf, Inf, Inf]
	    return xopt,∇fopt,∇fopt,fopt,niter,calls,OK
        end
        calls = [nlp.counters.neval_obj,  nlp.counters.neval_grad, nlp.counters.neval_hess, nlp.counters.neval_hprod]

        optimal, unbounded, tired, elapsed_time = stop(s,iter,x,f,∇f)
        #optimal = (norm_g < atol)| (norm_g <( rtol * norm_g0)) | (isinf(f) & (f<0.0))
        #tired = (iter >= itmax) | (sum(calls) > max_calls)
        stalled = unsuccinarow >= max_unsuccinarow

        finish = optimal || unbounded || tired || stalled
    end

    xopt = x
    fopt = f
    ∇fopt = ∇f
    norm_∇fopt = norm_∇f

    niter = verysucc+succ+unsucc;
    verbose && @printf(" %10i unsuccessful iterations \n",unsucc)
    verbose && @printf(" %10i successful iterations \n",succ)
    verbose && @printf(" %10i very successful iterations \n",verysucc)
    verbose && @printf(" %10i total successful iterations \n",verysucc+succ)
    verbose && @printf(" %10i total iterations \n",niter)

    OK = optimal

    calls = [nlp.counters.neval_obj, nlp.counters.neval_grad, nlp.counters.neval_hess, nlp.counters.neval_hprod]

    total_calls = nlp.counters.neval_obj + nlp.counters.neval_grad + n*nlp.counters.neval_hess+ nlp.counters.neval_hprod
     verbose && @printf(" %10i total f %10i g %10i H calls  %10i Hv calls  \n",nlp.counters.neval_obj, nlp.counters.neval_grad, nlp.counters.neval_hess, nlp.counters.neval_hprod)

    if OK && verbose  @printf("\n     Final cost(x) = %f  ||g|| = %9.2e  \#f/g/H evals = %i  converged\n",fopt, norm_∇fopt,total_calls)
    elseif verbose
        @printf("\n     Final cost(x) = %f  ||g|| = %9.2e  \#f/g/H evals  = %i  NOT converged\n\n",fopt, norm_∇fopt,total_calls)
    end;

    if     optimal   status = :Optimal
    elseif stalled   status = :Stalled
    elseif unbounded status = :Unbounded
    else             status = :UserLimit
    end

    return xopt,fopt,norm_∇fopt,norm(∇fopt),niter,calls,OK,status

end


text:
using Base.Test
import Base: return_types

using Compat

using Neovim
import Neovim: get_buffers, set_line, get_line
import Neovim: vim_eval, command, get_var, set_var
import Neovim: on_notify, on_request
nvim, proc = nvim_spawn()

#test buffer
buf = get_buffers(nvim)[1]
@assert isa(buf,Buffer)
set_line(buf, 0, "some text")
text = get_line(buf, 0)
@assert text == "some text"

#test high-level buffer interface
buf = current_buffer(nvim)
@assert isa(buf,Buffer)
@assert buf[:] == ["some text"]
@assert buf[1] == "some text"

buf[:] = ["alpha", "beta", "gamma"]
@assert buf[:] == ["alpha", "beta", "gamma"]
@assert buf[1] == "alpha"
@assert buf[1:3] == ["alpha", "beta", "gamma"]
@assert buf[1:1] == ["alpha"]
@assert buf[end-1:end] == ["beta", "gamma"]
@assert buf[1+end-1:end] == ["gamma"]
@assert buf[end-1+1:1+end-1] == ["gamma"]
@assert buf[end-1:2] == ["beta"]
@assert buf[end+1-1] == "gamma"
@assert buf[1+end-1] == "gamma"
@assert buf[2:1] == []
@assert buf[1:0] == []

@test_throws BoundsError buf[0]
@test_throws BoundsError buf[-1]
@test_throws BoundsError buf[end+1]
@test_throws BoundsError buf[-1:1]
@test_throws BoundsError buf[0:1]
@test_throws BoundsError buf[-1:0]
@test_throws BoundsError buf[end+1:1]
@test_throws BoundsError buf[1:end+1]
@test_throws BoundsError buf[1:1+end]
@test_throws BoundsError buf[1+end:1]
@test_throws BoundsError buf[end-1:end+1]
@test_throws BoundsError buf[end:end+1]
@test_throws BoundsError buf[end+1:end+1]

@assert (buf[2] = "beta-ish") == "beta-ish"
@assert buf[:] == ["alpha", "beta-ish", "gamma"]
@assert (buf[end] = "gamma-ish") == "gamma-ish"
@assert buf[:] == ["alpha", "beta-ish", "gamma-ish"]
@assert (buf[2:3] = ["b", "c", "d"]) == ["b", "c", "d"]
@assert buf[:] == ["alpha", "b", "c", "d"]
@assert (buf[end-1:end-2] = ["boom"]) == ["boom"]
@assert buf[:] == ["alpha", "b", "boom", "c", "d"]
@assert (buf[end-1:1+end-2] = ".") == "."
@assert buf[:] == ["alpha", "b", "boom", ".", "d"]

@test_throws BoundsError buf[0] = "twas"
@test_throws BoundsError buf[-1] = "brillyg"
@test_throws BoundsError buf[end+1] = "and"
@test_throws BoundsError buf[-1:1] = "the"
@test_throws BoundsError buf[0:1] = "slythy"
@test_throws BoundsError buf[-1:0] = "toves"
@test_throws BoundsError buf[end+1:1] = "did"
@test_throws BoundsError buf[1:end+1] = "gyre"
@test_throws BoundsError buf[1:1+end] = "and"
@test_throws BoundsError buf[1+end:1] = "gymble"
@test_throws BoundsError buf[end-1:end+1] = "in"
@test_throws BoundsError buf[end:end+1] = "the"
@test_throws BoundsError buf[end+1:end+1] = "wabe"

deleteat!(buf, 2)
@assert buf[:] == ["alpha", "boom", ".", "d"]
deleteat!(buf, 3:4)
@assert buf[:] == ["alpha", "boom"]

push!(buf, "the end")
@assert buf[:] == ["alpha", "boom", "the end"]
unshift!(buf, "stuff")
@assert buf[:] == ["stuff", "alpha", "boom", "the end"]
@assert length(buf) == 4

#test high-level cursor api
win = current_window(nvim)
cursor!(win, 1, 1)
@assert cursor(win) == (1, 1)
cursor!(win, 4, 2)
@assert cursor(win) == (4, 2)

#test eval
@assert vim_eval(nvim, "2+2") == 4

#test command and async behavior
command(nvim, "let g:test = []")
@sync for i in range(1,5)
    @async for j in range(1,5)
        command(nvim, "call add(g:test, [$i,$j])")
        if rand() > 0.5; sleep(0.001) end
    end
end

res = get_var(nvim, "test")
@assert length(res) == 25
for i in range(1,5)
    @assert all(indexin([[i,j] for j in range(1,5)], res) .> 0)
end

#test events
immutable TestHandler
    r::RemoteRef
end
on_notify(h::TestHandler, c, name, args) = put!(h.r, (name, args))
function on_request(h::TestHandler, c, serial, name, args)
    if name == "do_stuff"
        reply_result(c, serial, 10*args[1]+args[2])
    end
end

#notification
ref = RemoteRef()
nvim, proc = nvim_spawn(TestHandler(ref))
command(nvim, "call rpcnotify($(nvim.channel_id), 'mymethod', 10, 20)")
@assert take!(ref) == ("mymethod", Any[10, 20])

#request
@assert vim_eval(nvim, "100+rpcrequest($(nvim.channel_id), 'do_stuff', 2, 3)") == 123

# type stability of generated functions
@assert return_types(Neovim.get_buffers, (NvimClient,)) == [Vector{Buffer}]
@assert return_types(Neovim.command, (NvimClient,UTF8String)) == [Void]
@assert return_types(Neovim.get_current_line, (NvimClient,)) == [ByteString]
@assert return_types(Neovim.is_valid, (Tabpage,)) == [Bool]
@assert return_types(Neovim.get_height, (Window,)) == [Int]
@assert return_types(Neovim.get_mark, (Buffer,ASCIIString)) == [Tuple{Int,Int}]

#as ByteString isn't concrete anyway, this doesn't give that much really
#@assert return_types(Neovim.get_line_slice, (Buffer,Int,Int,Bool,Bool)) == [Vector{TypeVar(:_,None,ByteString)}]

# test host
hostdir = dirname(dirname(@__FILE__))
plugdir = joinpath(dirname(@__FILE__), "hosttest")
#fake initialization for :UpdateRemotePlugins
vimdir = mktempdir()
nvimrc = joinpath(vimdir, "nvimrc")
open(f->nothing,nvimrc,"w")
ENV["MYVIMRC"] = nvimrc
ENV["NEOVIM_JL_DEBUG"] = "templog"
rtp = "set rtp+=$hostdir,$plugdir"
juliap = "let g:julia_host_prog = '$(joinpath(JULIA_HOME, "julia"))'"
run(`nvim -u $nvimrc -i NONE --cmd $rtp --cmd $juliap -c UpdateRemotePlugins -c q`)
println("REGISTERED")
run(`cat templog`)
run(`cat $(joinpath(vimdir, ".nvimrc-rplugin~"))`)

try
ref = RemoteRef()
n,p = nvim_spawn(TestHandler(ref), cmd=`nvim --embed -u $nvimrc -i NONE --cmd $rtp --cmd $juliap`)

@assert vim_eval(n, "TestFun('a',3)") == "TestFun got a, 3"

command(n, "call AsyncFun($(n.channel_id), {'alfa':1, 'omega':'theend'})")
@assert take!(ref) == ("AsyncReply", Any[Dict("alfa"=>1, "omega"=>"theend")])

b = current_buffer(n)
b[1:5] = "line"
command(n, "2,3JLCommand text")
@assert b[:] == ["line","line","line","line","line","text"]
@assert get_var(n, "therange") == [2,3]

#TODO: more specific test? (somewhere in there should be "KABOOM!")
@test_throws ErrorException command(n, "Explode")

set_var(n, "zinged", 5)
command(n, "do User zing")
@assert get_var(n, "zinged") == 6
command(n, "do User yoink")
# shouldn't change
@assert get_var(n, "zinged") == 6

@assert vim_eval(n, "Global()") == 3
@assert vim_eval(n, "Local()") == 37

finally
#for debugging tests:
run(`cat templog`)
end


text:
export exactpruning!

# Exact pruning
"""
    exactpruning!(man::AbstractCutPruner, optimizer_constructor;
                  ub=Inf, lb=-Inf, epsilon=1e-5)

Remove dominated cuts in CutPruner `man`.

We use a LP solver to determine whether a cut is dominated or not.

# Arguments
* `man::AbstractCutPruner`
    Cut pruner where to remove cuts
* `optimizer_constructor`
    Optimizer used to solve LP
* `ub::Union{Float64, Vector{Float64}}`
    State x upper bound
* `lb::Union{Float64, Vector{Float64}}`
    State x lower bound
* `epsilon::Float64`
    Pruning's tolerance
"""
function exactpruning!(man::AbstractCutPruner, optimizer_constructor; ub = Inf, lb = -Inf, epsilon = 1e-5)
    K = getdominated(man.A, man.b, man.islb, man.isfun, optimizer_constructor, lb, ub, epsilon)
    removecuts!(man, K)
end

"""Return dominated cuts."""
function getdominated(A, b, islb, isfun, optimizer_constructor, lb, ub, epsilon)
    red = Int[]
    if size(A, 1) == 1
        return red
    end
    for i in 1:size(A, 1)
        if isdominated(A, b, islb, isfun, i, optimizer_constructor, lb, ub, epsilon)
            push!(red, i)
        end
    end
    red
end

"""State whether a cut is dominated with a tolerance epsilon."""
function isdominated(A, b, islb, isfun, k, optimizer_constructor, lb, ub, epsilon)
    # we use JuMP to solve the test
    # For instance, if islb & isfun, the LP becomes:
    # min - y
    #     x ∈ R^n, y ∈ R
    #     y <= (A[k, :] - A[i, :])*x + b[k] - b[i]     ∀ i != k

    # we formulate an equivalent problem as
    # min c'*z
    # s.t H*z <= h

    # get problem's dimension
    ncuts, nx = size(A)
    # allocate arrays
    h = zeros(ncuts - 1)
    H = zeros(ncuts - 1, nx + 1)
    c = zeros(nx + 1)
    c[1] = (islb) ? -1 : 1

    λk = @view A[k, :]

    ic = 0
    for ix in 1:ncuts
        if ix != k
            ic += 1
            δb = b[k] - b[ix]
            @inbounds h[ic] = ((isfun & islb) || (~isfun & ~islb)) ? δb : -δb
            @inbounds H[ic, 1] = (islb) ? 1 : -1

            for jx in 1:nx
                dl = A[ix, jx] - λk[jx]
                @inbounds H[ic, jx + 1] = (islb) ? dl : -dl
            end
        end
    end

    # update lower and upper bound if Vector
    lbx = isa(lb, Vector) ? vcat(-Inf, lb) : lb
    ubx = isa(ub, Vector) ? vcat(Inf, ub) : ub

    # solve the LP with JuMP
    model = Model(optimizer_constructor)
    z = @variable(model, [1:nx+1])
    if lb isa Vector
        set_lower_bound.(z, lb)
    end
    if ub isa Vector
        set_upper_bound.(z, ub)
    end

    @constraint(model, H * z .≤ h)
    @objective(model, Min, c ⋅ z)
    optimize!(model)

    status = termination_status(model)

    if status == MOI.INFEASIBLE
        return true
    elseif status == MOI.DUAL_INFEASIBLE
        return false
    elseif status == MOI.OPTIMAL
        res = objective_value(model)
        return (islb) ? -res < epsilon : res > -epsilon
    else
        error("Solver returned status $status: $(raw_status(model)).")
    end
end


text:
"""
  `Sorted` should be the first argument to renorm
   when the values to be re-normalized are given
   in order of nonincreasing absolute magnitude.
   The value with the largest exponent is the most
   significant value and appears as the first one;
   the value with the smallest exponent (more negative)
   is the least significant value and appears after
   any other __nonzero__ contributors the the
   extended precision significand.
"""
struct Sorted end

renorm(a::T) where T<:Real = a

renorm(a::T, b::T) where T<:Real = two_sum(a, b)

renorm(::Type{Sorted}, a::T, b::T) where T<:Real = fast_two_sum(a, b)

function renorm(a::T,b::T,c::T) where {T<:AbstractFloat}
    s, t = two_sum(b, c)
    x, u = two_sum(a, s)
    y, z = two_sum(u, t)
    x, y = two_sum(x, y)
    return x, y, z
end

function renorm_(::Type{Sorted}, a::T,b::T,c::T) where {T<:Real}
    s, t = fast_two_sum(b, c)
    x, u = fast_two_sum(a, s)
    y, z = fast_two_sum(u, t)
    x, y = fast_two_sum(x, y)
    return x, y, z
end

function renorm_3(::Type{Sorted}, ahi::T, amd::T, alo::T) where {T<:Real}
    amd, alo = fast_two_sum(amd, alo)
    ahi, amd = fast_two_sum(ahi, amd)
    amd, alo = fast_two_sum(amd, alo)
    return ahi, amd, alo
end


text:
function gettext(x)
    return fold(x, "") do s, x
        if x isa XLeaf && x[] isa AbstractString
            return s * x[]
        else
            return s
        end
    end
end


text:
# convienence routines
help() = execute("help")
help(s::AbstractString) = execute("help, "*s)
idlhelp(s::AbstractString) = execute("?"*s)
idlhelp{T<:AbstractString}(strarr::Array{T,1}) = println("IDL.idlhelp: Array input not supported")
shell_command(s::AbstractString) = println("% Shell commands not allowed in IDLRPC")
shell_command{T<:AbstractString}(strarr::Array{T,1}) = println("% Shell commands not allowed in IDLRPC")
reset() = execute(".reset_session")
full_reset() = execute(".full_reset_session")
dotrun(filename::AbstractString) = execute(".run $filename")

function execute(str::AbstractString)
    ok, strarr, msg = convert_command(str)
    ok || error(msg)
    execute_converted(strarr)
    return nothing
end

function execute_converted{T<:AbstractString}(strarr::Array{T,1})
    # does no conversion of interpolated vars, continuation chars, or newlines
    for str in strarr
        execute_converted(str) || return false
    end
    return true
end

function put_var_from_name(name::AbstractString, abort::Bool=true)
    # abort=false causes routine to not issue error
    ok = true
    msg = ""
    if !isdefined(Main, Symbol(name))
        ok = false
        msg = "IDL.put_var_from_name: undefined variable $name in Module Main"
        if abort
            error(msg)
        else
            return (ok, msg)
        end
    end
    put_var(getfield(Main, Symbol(name)), name)
    return (ok, msg)
end

function put_var{T<:AbstractString,N}(arr::Array{T,N}, name::AbstractString)
    # Sort of a HACK: do direcly since ImportNamedArray doesn't work
    execute("$name = strarr"*replace(string(size(arr)), ",)", ")"))
    for i=1:length(arr)
        j = i-1
        str = arr[i]
        execute("$name[$j] = '$str'")
    end
    return
end

const mask64 = 0x0000000000000000ffffffffffffffff
const mask32 = 0x000000000000000000000000ffffffff

function get_var(vptr::Ptr{IDL_Variable}, name::AbstractString="")
    var = unsafe_load(vptr)
    # some types not dealt with
    if (var.flags & IDL_V_FILE) != 0
        error("IDL.extract_from_vptr: $name: assoc type not setup")
    end
    ## if (var.flags & IDL_V_DYNAMIC) != 0
    ##     println("dynamic")
    ## end
    if (var.flags & IDL_V_NULL) != 0
        error("IDL.extract_from_vptr: $name: variable is null")
    end

    # array types
    if (var.flags & IDL_V_ARR) != 0
        if var.vtype == IDL_TYP_STRING
            parr = reinterpret(Ptr{IDL_Array}, convert(Int, var.buf))
            idl_arr = unsafe_load(parr)
            pdata = reinterpret(Ptr{IDL_String}, idl_arr.data)
            strarr = Array(Compat.ASCIIString, dims(idl_arr.dim, idl_arr.n_dim))
            for i=1:idl_arr.n_elts
                data = unsafe_load(pdata, i)
                strarr[i] = data.slen > 0 ? unsafe_string(data.s, Int(data.slen)) : ""
            end
            return strarr
        elseif var.vtype == IDL_TYP_STRUCT
            error("IDL.extract_from_vptr: $name: STRUCT not setup")
        elseif var.vtype == IDL_TYP_PTR
            error("IDL.extract_from_vptr: $name: PTRARR types not setup")
        elseif var.vtype == IDL_TYP_OBJREF
            error("IDL.extract_from_vptr: $name: OBJARR types not setup")
        else
            parr = reinterpret(Ptr{IDL_Array}, convert(Int, var.buf))
            idl_arr = unsafe_load(parr)
            jl_t = jl_type(var.vtype)
            pdata = reinterpret(Ptr{jl_t}, idl_arr.data)
            # not sure why this doesn't work
            # arr = unsafe_wrap(Array, pdata, dims(idl_arr.dim, idl_arr.n_dim))
            arr = Array(jl_t, dims(idl_arr.dim, idl_arr.n_dim))
            for i=1:idl_arr.n_elts
                arr[i] = unsafe_load(pdata, i)
            end
            return arr
        end
    end

    # Scalar value
    if var.vtype == IDL_TYP_UNDEF
        error("IDL.extract_from_vptr: $name: undefined variable")
    elseif var.vtype == IDL_TYP_BYTE
        return reinterpret(Int8, convert(UInt8, var.buf))
    elseif var.vtype == IDL_TYP_INT
        return reinterpret(Int16, convert(UInt16, var.buf))
    elseif var.vtype == IDL_TYP_LONG
        return reinterpret(Int32, convert(UInt32, var.buf))
    elseif var.vtype == IDL_TYP_FLOAT
        return reinterpret(Float32, convert(UInt32, var.buf))
    elseif var.vtype == IDL_TYP_DOUBLE
        return reinterpret(Float64, convert(UInt64, var.buf))
    elseif var.vtype == IDL_TYP_COMPLEX
        return complex(reinterpret(Float32, convert(UInt32, var.buf & mask32)),
                       reinterpret(Float32, convert(UInt32, var.buf >> 32)))
    elseif var.vtype == IDL_TYP_STRING
        slen = reinterpret(Int32, convert(UInt32,var.buf & mask32))
        stype = reinterpret(Int32, convert(UInt32,(var.buf & mask64) >> 32))
        println(stype)
        s = reinterpret(Ptr{Cchar}, convert(UInt64,var.buf >> 64))
        return slen > 0 ? unsafe_string(s, slen) : ""
    elseif var.vtype == IDL_TYP_STRUCT
        error("IDL.extract_from_vptr: $name: STRUCT not setup")
    elseif var.vtype == IDL_TYP_DCOMPLEX
        return complex(reinterpret(Float64, convert(UInt64, var.buf & mask64)),
                       reinterpret(Float64, convert(UInt64, var.buf >> 64)))
    elseif var.vtype == IDL_TYP_PTR
        error("IDL.extract_from_vptr: $name: PTR not setup")
    elseif var.vtype == IDL_TYP_OBJREF
        error("IDL.extract_from_vptr: $name: OBJREF not setup")
    elseif var.vtype == IDL_TYP_UINT
        return reinterpret(UInt16, convert(UInt16, var.buf))
    elseif var.vtype == IDL_TYP_ULONG
        return reinterpret(UInt32, convert(UInt32, var.buf))
    elseif var.vtype == IDL_TYP_LONG64
        return reinterpret(Int64, convert(UInt64, var.buf))
    elseif var.vtype == IDL_TYP_ULONG64
        return reinterpret(UInt64, convert(UInt64, var.buf))
    end
    # should be impossible to get here
    error("IDL.extract_from_vptr: $name: type is not setup")
end

function inside_string(pt::Int, line::AbstractString)
    for re in [r"('[^']+')", r"(\"[^\"]+\")"]
        for m in eachmatch(re, line)
            if pt >= m.offset &&
                pt < m.offset+endof(m.captures[1])
                return true
            end
        end
    end
    return false
end

function convert_continuations(line)
    # remove trailing comments and continuation lines
    # will remove continuation on final line which is invalid idl syntax
    pt = start(line)
    while (pt = first(search(line, r";|\$", pt))) > 0
        if !inside_string(pt, line)
            line = replace(line, r"(;|\$).*(\n|$)", "", 1)
        end
        if pt < endof(line)
            pt = next(line, pt)[2]
        end
    end
    return true, line, ""
end

function convert_newlines(line)
    # separates line at newline ("\n") characters into string array
    # assumes that continuation characters are already removed
    # not type stable but should not matter for repl use
    line = chomp(line)
    if in('\n', line)
        line = split(line, '\n', keep=false)
    end
    return line
end

function replace_interpolated_vars(line)
    # use %var to automatically pull var into idl and use it
    while (m = match(r"\%(\w+)", line)) != nothing
        if !inside_string(m.offset, line)
            ok, msg = put_var_from_name(ascii(m.captures[1]), false)
            if !ok
                return false, line, msg
            end
            line = line[1:m.offset-1]*line[m.offset+1:end]
        end
    end
    return true, line, ""
end

function convert_command(line)
    ok, line, msg = replace_interpolated_vars(line)
    ok || return ok, line, msg
    ok, line, msg = convert_continuations(line)
    ok || return ok, line, msg
    line = convert_newlines(line)
    return true, line, msg
end


text:
"""
The base type of a fractional integral sense, in FractionalCalculus.jl, all of the fractional integral senses belong to ```FracIntAlg```
"""
abstract type FracIntAlg end

include("RL.jl")

"""
    fracint(f, α, point, h, FracIntAlg())

The general API used for fractional integral computing.
"""
function fracint(f, α, point, h, ::FracIntAlg)
end

## Macros for convenient computing.
"""
    @fracint(f, α, point)

Return the α-order integral of **f** at specific point.

```julia-repl
julia> @fracint(x->x, 0.5, 1)
0.7522525439593486
```
"""
macro fracint(f, α, point)
    return :(fracint($f, $α, $point, 0.0001, RLInt_Approx()))
end

"""
    @semifracint(f, point)

Return the semi-integral of **f** at specific point.

```julia-repl
julia> @semifracint(x->x, 1)
0.7522525439593486
```
"""
macro semifracint(f, point)
    return :(fracint($f, 0.5, $point, 0.0001, RLInt_Approx()))
end