text:
using Random

"""
Linear channel coding methods for frequency detection; specifically to

1. create generator matrices
2. BSC decode the results
"""

get_code_lookup = Dict()

"""
Estimate the frequency from a singleton given an offset coding scheme and generator matrix G.

Arguments
---------
recvd_codeword : BitArray{1}
The sampled vector from which we are decoding

G : BitArray{2}
The generator matrix fo the code

code : Symbol
The code to use. All methods referenced must use this signature (minus "code".)

Returns
-------
k : BitArray{1}
The estimated frequency.
"""
function decode_with(recvd_codeword::BitArray{1}, G::BitArray{2}, code::Symbol)::BitArray{1}
    return get_code_lookup[code](recvd_codeword, G)
end

#LDPC Code Based on https://github.com/christianpeel/LDPCsim.jl/blob/master/src/LDPCsim.jl

function generate_LDPC_adjacency_matrix(M::Int64, N::Int64)
    # Create LDPC adjacency matrix of the Tanner graph H
    # Note that rate = M/N
    # Design rate is 1 - j/k where j is the degree of variable nodes, k degree of check nodes
    # i.e. j ones per column, k ones per row. Rounding k up means no need to headache about running out of spots
    # Pick j = 3 from convention in the literature (smallest j with good theoretical properties)
    # So with M/N = 1 - 3/k, we pick k = ceil(3 / (1 - M/N))

    # However, we must have 3N ones in total, so we actually end up using a bigger k sometimes
    # Notably, we use the larger of k or 3N/M ones per row
    # Theoretically, this actually results in a better rate (though more expensive decoding process)

    # Alg is:
    # For each column, pick three random positions to put 1's s.t. no row ends up with more than k 1's in it
    # Inputs
    #  M        : num of rows (frequently k in literature, the dimension of the linear code)
    #  N        : num of columns
    #
    # Output
    #  H        : LDPC adjacency matrix   (N - k) x N = (N - M) x N

    j = 3
    M = N - M # dimension generated here is actually N - M
    H = falses(M, N)
    num_ones = zeros(Int64, M) #number of ones per row so far
    for current_col = 1:N
        # For uniformity, we want to have 3N/M ones in total at the end
        # So as we expand, say we've filled N' columns so far. We have 3N'
        # ones, so 3N'/M ones per row. So, to preserve uniformity as we go,
        # add to rows with less than ceil(3(N'+1)/M) ones where N' is the current
        # column we're filling. current_col is N' here.

        legal = [i for (i, val) in enumerate(num_ones) if val < ((3 * current_col / M))]
        if size(legal)[1] < j
            # Need j - size(legal) more spots, pick them uniformly randomly
            d = shuffle(collect(setdiff(Set(1:M), Set(legal))))
            additional_options = d[1:j - size(legal)[1]]
            append!(legal, additional_options)
        end

        choices = shuffle(legal)[1:j]
        for idx in choices
            H[idx, current_col] = 1
            num_ones[idx] += 1
        end
    end

    # An additional optimization could be to expand the girth, but expensive and probably unnecessary
    return H
end

function generate_LDPC_generator_matrix(H::BitArray)
    # https://en.wikipedia.org/wiki/Low-density_parity-check_code#Operational_use
    # has instructions for what we're doing, but essentially:
    # 1. Gaussian eliminate s.t. H = [P | I] as a block matrix
    # 2. Return G = [I | P^\top]
    # Returns nothing if H has no corresponding generator (rank-deficient)
    H = copy(H)
    M = size(H)[1]
    N = size(H)[2]
    pivot_row = 1
    pivot_col = N - M + 1

    for _ in 1:M
        # println("In pivot positon (", pivot_row, ", ", pivot_col, ")")
        # 1. Swap a row with a 1 into this pivot position if need be
        if !H[pivot_row, pivot_col]
            success = false
            for j in (pivot_row + 1):M
                if H[j, pivot_col]
                    # switch the found row and the pivot row
                    for k in 1:N
                        H[j, k], H[pivot_row, k] = H[pivot_row, k], H[j, k]
                    end
                    success = true
                    break
                end
            end
            if !success
                return nothing #Rank failure if we reach here
            end
        end

        # 2. Kill the 1s below this in this column
        for j in (pivot_row + 1):M
            if H[j, pivot_col]
                for k in 1:N
                    H[j, k] = xor(H[j, k], H[pivot_row, k])
                end
            end
        end

        # 3. Kill the 1s above this
        for j in 1:(pivot_row - 1)
            if H[j, pivot_col]
                for k in 1:N
                    H[j, k] = xor(H[j, k], H[pivot_row, k])
                end
            end
        end

        pivot_row += 1
        pivot_col += 1
    end

    # Extract P
    k = N - M
    P = transpose(H[:, 1:k])

    # Construct and return G
    G = falses(k, N)
    G[1:k, (k + 1):N] = P
    for i in 1:k
        G[i, i] = 1
    end
    return G
end

function generate_LDPC_code(M::Int64, N::Int64)
    H = nothing
    G = nothing
    # Sometimes random check matrices will have insufficient rank (no corresponding generator)
    # So we just resample until it's good. Ideally the probability of insufficient rank goes away for larger sizes,
    # something to mention/think about later than I am writing this.
    while G == nothing
        H = generate_LDPC_adjacency_matrix(M, N)
        G = generate_LDPC_generator_matrix(H)
    end
    return H, G
end

function check_pair(G, H)
    C = G * transpose(H)
    for i in 1:size(C)[1]
        for j in 1:size(C)[2]
            if C[i, j] % 2 != 0
                return false
            end
        end
    end
    return true
end

function do_test(M::Int64, N::Int64, iters::Int64)
    for _ in 1:iters
        H, G = generate_LDPC_code(M, N)
        println(size(H), size(G))
        if !check_pair(G, H)
            println("FAIL")
        end
    end
end

function encode_LDPC(input::BitVector, G::BitArray)
    return BitVector(vec((transpose(input) * G) .% 2))
end

function decode_LDPC(y::BitVector, H::BitArray, p::Float64, iters::Int64=30)
    # y is the received codeword of length N
    # H is the check matrix as usual of dimensions (N - k) x N
    N = size(H)[2]
    k = N - size(H)[1]
    M = N - k

    variableLLR = zeros(N) # LLR for the variable nodes
    checkLLR = zeros(M) # LLR for the check nodes
    messagesToCheckNodes = zeros(N, M) # Messages passed from variable nodes to check nodes. Omega matrix
    messagesToVarNodes = zeros(M, N) # Messages passed from check nodes to variable nodes. Lambda matrix

    # 1. Initialize the variable nodes based on received information (the identity part of the codeword)
    # We initialize each variable node's LLR
    # L_i = ln(P(bit = 1|observed)/P(bit = 0|observed)) - positive (>= 0) for 1, negative for 0
    # Assuming uniform prior, L_i = ln(p/(1-p)) when y_i = 0 and ln((1-p)/p) when y_i = 1
    lnq = log(p) - log(1-p)
    for i in 1:N
        variableLLR[i] = y[i] ? -lnq : lnq
    end

    # Initial messages are just the starting beliefs
    for i in 1:N
        for j in 1:M
            if H[j, i]
                messagesToCheckNodes[i, j] = variableLLR[i]
            end
        end
    end

    # print("My initial beliefs are ")
    # display(variableLLR)
    # println()

    # print("My initial messages are ")
    # display(messagesToCheckNodes)
    # println()

    # 30 rounds of decoding based on literature and practice
    for iter in 1:iters
        # 2. Update the check node LLRs and generate new check to variable messages
        for i in 1:M
            checkLLR[i] = 0
            sign = 0 # even is negative, odd is true as per usual
            for j in 1:N
                if H[i, j]
                    # save the parity of the variable nodes' signs
                    sign += (messagesToCheckNodes[j, i] >= 0) ? 1 : 0
                    # tanh clamps to [0, 1], so each contribution is negative, and the reliability is the opposite of this
                    checkLLR[i] += log(tanh(abs(messagesToCheckNodes[j, i]) / 2))
                end
            end
            if mod(sign, 2) == 1 # make sure parity lines up
                checkLLR[i] *= -1
            end
        end

        # print("My check node LLRs are ")
        # display(checkLLR)
        # println()

        for i in 1:M
            for j in 1:N
                if H[i, j]
                    messagesToVarNodes[i, j] = 2 * atanh(exp(-abs(checkLLR[i]) - log(tanh(abs(messagesToCheckNodes[j, i] / 2)))))
                    if !xor(checkLLR[i] >= 0, messagesToCheckNodes[j, i] >= 0)
                        messagesToVarNodes[i, j] *= -1
                    end
                end
            end
        end

        # 3. Update the variable nodes and generate new variable to check messages
        for i in 1:N
            for j in 1:M
                if H[j, i]
                    variableLLR[i] += messagesToVarNodes[j, i]
                end
            end
        end

        for i in 1:N
            for j in 1:M
                messagesToCheckNodes[i, j] = variableLLR[i] - messagesToVarNodes[j, i]
            end
        end

        # 4. Quantize and halt
        estimate = falses(N)
        for i in 1:N
            if variableLLR[i] >= 0
                estimate[i] = true
            end
        end

        if BitVector(vec((transpose(estimate) * transpose(H)) .% 2)) == falses(M)
            return estimate[1:k]
        end
    end

    # Failure to decode
    return nothing
end

using Distributions

k = 100
N = 200
H, G = generate_LDPC_code(k, N)
# display(H)
input = BitVector(rand(Binomial(1, 0.5), k))
codeword = encode_LDPC(input, G)
p = 0.05
d = Binomial(1, p)
noise = BitVector(rand(d, N))
for i in 1:size(codeword)[1]
    codeword[i] = xor(codeword[i], noise[i])
end
recovered = decode_LDPC(codeword, H, p, 200)
println(input)
println(noise)
println(recovered)


text:
import Pkg
# Pkg.pkg"add CairoMakie#jk/scatter-glyphs AbstractPlotting#jk/text-layouting"
using Documenter, Literate
# avoid font caching warning in docs
using AbstractPlotting, CairoMakie, MakieLayout
CairoMakie.activate!()
scatter(rand(10), rand(10))

# generate examples
GENERATED = joinpath(@__DIR__, "src", "literate")
SOURCE_FILES = joinpath.(GENERATED, ["examples.jl"])
foreach(fn -> Literate.markdown(fn, GENERATED), SOURCE_FILES)

makedocs(
         sitename="MakieRecipes",
         pages = Any[
                     "index.md",
                     "Examples" => "literate/examples.md",
                    ]
        )

deploydocs(
    repo = "github.com/JuliaPlots/MakieRecipes.jl.git",
    push_preview = true
)

ext: jl


text:
function area_balance(
    container::OptimizationContainer,
    expression::ExpressionKey,
    area_mapping::Dict{String, Array{PSY.Bus, 1}},
    branches,
)
    time_steps = get_time_steps(container)
    nodal_net_balance = get_expression(container, expression)

    constraint = add_constraints_container!(
        container,
        AreaDispatchBalanceConstraint(),
        PSY.Area,
        keys(area_mapping),
        time_steps,
    )
    area_balance = get_variable(container, ActivePowerVariable(), PSY.Area)
    for (k, buses_in_area) in area_mapping
        for t in time_steps
            area_net =
                built_for_recurrent_solves(container) ? zero(PGAE) : JuMP.AffExpr(0.0)
            for b in buses_in_area
                JuMP.add_to_expression!(area_net, nodal_net_balance[PSY.get_number(b), t])
            end
            constraint[k, t] =
                JuMP.@constraint(container.JuMPmodel, area_balance[k, t] == area_net)
        end
    end

    expr_up = get_expression(container, EmergencyUp(), PSY.Area)
    expr_dn = get_expression(container, EmergencyDown(), PSY.Area)

    participation_assignment_up = add_constraints_container!(
        container,
        AreaParticipationAssignmentConstraint(),
        PSY.Area,
        keys(area_mapping),
        time_steps,
        meta="up",
    )
    participation_assignment_dn = add_constraints_container!(
        container,
        AreaParticipationAssignmentConstraint(),
        PSY.Area,
        keys(area_mapping),
        time_steps,
        meta="dn",
    )

    for area in keys(area_mapping), t in time_steps
        participation_assignment_up[area, t] =
            JuMP.@constraint(container.JuMPmodel, expr_up[area, t] == 0)
        participation_assignment_dn[area, t] =
            JuMP.@constraint(container.JuMPmodel, expr_dn[area, t] == 0)
    end

    return
end


text:
function coordinategenes(
    chromosome::AbstractString,
    range1::Real,
    range2::Real,
    gencode::DataFrame,
    height::Real
    )

    df = filter(x -> (x.seqnames == chromosome) && (x.end >= range1) && (x.start <= range2), gencode)
    dfg = view(df, df.feature .== "gene", :)
    dfe = view(df, df.feature .== "exon", :)
    genes = unique(dfg.gene_name)
    strand = [dfg.strand[findfirst(isequal(gene), dfg.gene_name)] for gene in genes]
    ps = Vector{Vector{Polygon}}(undef, length(genes))
    bs = Matrix{Float64}(undef, length(genes), 2)
    rows = ones(Int, length(genes))
    prop =  23000 * (range2 - range1) / 2.2e6
    for j in eachindex(genes)
        ind = findfirst(isequal(genes[j]), dfg.gene_name)
        start1 = dfg.start[ind]
        stop1 = dfg[ind, :end]
        center1 = (start1 + stop1) / 2
        label1 = center1 + (length(genes[j]) + 1) * prop
        bs[j, 1] = start1
        bs[j, 2] = stop1
        for k in (j + 1):length(genes)
            ind = findfirst(isequal(genes[k]), dfg.gene_name)
            start2 = dfg.start[ind]
            stop2 = dfg[ind, :end]
            center2 = (start2 + stop2) / 2
            label2 = center2 - (length(genes[k]) + 1) * prop
            if ((stop1 > start2) || (label1 > label2)) && (rows[j] == rows[k])
                rows[k] = rows[j] + 1
                while true
                    l = findprev(isequal(rows[k]), rows, k - 1)
                    isnothing(l) && break
                    start3 = bs[l, 1]
                    stop3 = bs[l, 2]
                    center3 = (start3 + stop3) / 2
                    label3 = center3 + (length(genes[l]) + 1) * prop
                    ((stop3 < start2) && (label3 < label2)) && break
                    rows[k] = rows[l] + 1
                end
            end
        end
    end
    for j in eachindex(genes)
        ranges = view(dfe, findall(isequal(genes[j]), dfe.gene_name), [:start, :end])
        n = size(ranges, 1)
        p = Vector{Polygon}(undef, n)
        for i = 1:n
            p[i] = Polygon([
                Point2f(ranges[i, 1], 1 - height - (rows[j] - 1) * (0.25 + height)),
                Point2f(ranges[i, 1], 1 - (rows[j] - 1) * (0.25 + height)),
                Point2f(ranges[i, 2], 1 - (rows[j] - 1) * (0.25 + height)),
                Point2f(ranges[i, 2], 1 - height - (rows[j] - 1) * (0.25 + height))
            ])
        end
        ps[j] = p
    end
    return genes, strand, ps, bs, rows
end

"""
    plotgenes!(ax::Axis, chromosome::AbstractString, range1::Real, range2::Real, gencode::DataFrame; kwargs)
    plotgenes!(ax::Axis, chromosome::AbstractString, bp::Real, gencode::DataFrame; kwargs)
    plotgenes!(ax::Axis, gene::AbstractString, gencode::DataFrame; kwargs)

Plot collapsed gene bodies for genes within a given `chromosome` and genomic range
between `range1` and `range2`.

Alternatively, plot within a given `chromosome` and a certain `window` around a
genomic coordinate `bp` or plot within a certain `window` around `gene`.

# Arguments
- `height::Real = 0.25`: the height of exons.
- `genecolor = :royalblue`: the color of genes.
- `textcolor = :black`: the color of gene labels.
- `window::Real = 1e6`: the window around `bp` or `gene`.
"""
function plotgenes!(
    ax::Axis,
    chromosome::AbstractString,
    range1::Real,
    range2::Real,
    gencode::DataFrame;
    height::Real = 0.25,
    genecolor = :royalblue,
    textcolor = :black
    )

    genes, strand, ps, bs, rows = coordinategenes(chromosome, range1, range2, gencode, height)
    if length(rows) == 0
        ax.spinewidth = 0.75
        hidexdecorations!(ax)
        hideydecorations!(ax)
        xlims!(ax, range1, range2)
        ylims!(ax, 0.875 - height, 1.375)
        rs = 18 * (0.25 + height) / 0.5
        return rs
    end
    for j in 1:size(ps, 1)
        poly!(ax, ps[j], color = genecolor, strokewidth = 0)
        lines!(ax, [bs[j, 1], bs[j, 2]],
            [1 - height / 2 - (rows[j] - 1) * (0.25 + height), 1 - height / 2 - (rows[j] - 1) * (0.25 + height)],
            color = genecolor, linewidth = 0.5)
        g = (strand[j] == "+" ? genes[j] * "→" : "←" * genes[j])
        text!(ax, "$g",
            position = ((bs[j, 1] + bs[j, 2]) / 2, 1 - (rows[j] - 1) * (0.25 + height)),
            align = (:center, :bottom), textsize = 6, color = textcolor)
    end
    ax.spinewidth = 0.75
    hidexdecorations!(ax)
    hideydecorations!(ax)
    xlims!(ax, range1, range2)
    ylims!(ax, 0.875 - height  - (maximum(rows) - 1) * (0.25 + height), 1.375)
    rs = 18 * maximum(rows) * (0.25 + height) / 0.5
    return rs
end

plotgenes!(ax::Axis, chromosome::AbstractString, bp::Real, gencode::DataFrame; window::Real = 1e6, kwargs...) =
    plotgenes!(ax, chromosome, bp - window, bp + window, gencode; kwargs...)

function plotgenes!(ax::Axis, gene::AbstractString, gencode::DataFrame; window::Real = 1e6, kwargs...)
    chr, start, stop = findgene(gene::AbstractString, gencode::DataFrame)
    plotgenes!(ax, chr, start - window, stop + window, gencode; kwargs...)
end

"""
    plotgenes!(ax::Axis, chromosome::AbstractString, range1::Real, range2::Real, highlight::Tuple{AbstractVector, AbstractVector}, gencode::DataFrame; height::Real)
    plotgenes!(ax::Axis, chromosome::AbstractString, bp::Real, highlight::Tuple{AbstractVector, AbstractVector}, gencode::DataFrame; window::Real, height::Real)
    plotgenes!(ax::Axis, gene::AbstractString, highlight::Tuple{AbstractVector, AbstractVector}, gencode::DataFrame; window::Real, height::Real)

Plot gene bodies with a vector of genes highlighted by a vector of colors via `highlight`.
"""
function plotgenes!(
    ax::Axis,
    chromosome::AbstractString,
    range1::Real,
    range2::Real,
    highlight::Tuple{AbstractVector, AbstractVector},
    gencode::DataFrame;
    height::Real = 0.25
)

    genes, strand, ps, bs, rows = coordinategenes(chromosome, range1, range2, gencode, height)
    if length(rows) == 0
        ax.spinewidth = 0.75
        hidexdecorations!(ax)
        hideydecorations!(ax)
        xlims!(ax, range1, range2)
        ylims!(ax, 0.875 - height, 1.375)
        rs = 18 * (0.25 + height) / 0.5
        return rs
    end
    for j in 1:size(ps, 1)
        ind = findfirst(isequal(genes[j]), highlight[1])
        if !isnothing(ind)
            poly!(ax, ps[j], color = highlight[2][ind], strokewidth = 0)
            lines!(ax, [bs[j, 1], bs[j, 2]],
                [1 - height / 2 - (rows[j] - 1) * (0.25 + height), 1 - height / 2 - (rows[j] - 1) * (0.25 + height)],
                color = highlight[2][ind], linewidth = 0.5)
            g = (strand[j] == "+" ? genes[j] * "→" : "←" * genes[j])
            text!(ax, "$g",
                position = ((bs[j, 1] + bs[j, 2]) / 2, 1 - (rows[j] - 1) * (0.25 + height)),
                align = (:center, :bottom), textsize = 6, color = highlight[2][ind])
        else
            poly!(ax, ps[j], color = "gray60", strokewidth = 0)
            lines!(ax, [bs[j, 1], bs[j, 2]],
                [1 - height / 2 - (rows[j] - 1) * (0.25 + height), 1 - height / 2 - (rows[j] - 1) * (0.25 + height)],
                color = "gray60", linewidth = 0.5)
            g = (strand[j] == "+" ? genes[j] * "→" : "←" * genes[j])
            text!(ax, "$g",
                position = ((bs[j, 1] + bs[j, 2]) / 2, 1 - (rows[j] - 1) * (0.25 + height)),
                align = (:center, :bottom), textsize = 6, color = "gray60")
        end
    end
    ax.spinewidth = 0.75
    hidexdecorations!(ax)
    hideydecorations!(ax)
    xlims!(ax, range1, range2)
    ylims!(ax, 0.875 - height  - (maximum(rows) - 1) * (0.25 + height), 1.375)
    rs = 18 * maximum(rows) * (0.25 + height) / 0.5
    return rs
end

plotgenes!(ax::Axis, chromosome::AbstractString, bp::Real, highlight::Tuple{AbstractVector, AbstractVector}, gencode::DataFrame; window::Real = 1e6, kwargs...) =
    plotgenes!(ax, chromosome, bp - window, bp + window, highlight, gencode; kwargs...)

function plotgenes!(ax::Axis, gene::AbstractString, highlight::Tuple{AbstractVector, AbstractVector}, gencode::DataFrame; window::Real = 1e6, kwargs...)
    chr, start, stop = findgene(gene::AbstractString, gencode::DataFrame)
    plotgenes!(ax, chr, start - window, stop + window, highlight, gencode; kwargs...)
end

plotgenes!(ax::Axis, chromosome::AbstractString, range1::Real, range2::Real, highlight::Tuple{AbstractString, Any}, gencode::DataFrame; kwargs...) =
    plotgenes!(ax, chromosome, range1, range2, ([highlight[1]], [highlight[2]]), gencode; kwargs...)

plotgenes!(ax::Axis, chromosome::AbstractString, bp::Real, highlight::Tuple{AbstractString, Any}, gencode::DataFrame; kwargs...) =
    plotgenes!(ax, chromosome, bp, ([highlight[1]], [highlight[2]]), gencode; kwargs...)

plotgenes!(ax::Axis, gene::AbstractString, highlight::Tuple{AbstractString, Any}, gencode::DataFrame; kwargs...) =
    plotgenes!(ax, gene, ([highlight[1]], [highlight[2]]), gencode; kwargs...)

"""
    labelgenome(g::GridPosition, chromosome::AbstractString, range1::Real, range2::Real)

Label `g` with a given `chromosome` and genomic range between `range1` and `range2`.
"""
function labelgenome(g::GridPosition, chromosome::AbstractString, range1::Real, range2::Real)
    Label(g, "~$(round(range1 / 1e6; digits = 1)) Mb", textsize = 6, halign = :left)
    Label(g, "Chr $(chromosome)", textsize = 6, halign = :center)
    Label(g, "~$(round(range2 / 1e6; digits = 1)) Mb", textsize = 6, halign = :right)
end


text:
using ITensors,
      Test,
      Random,
      LinearAlgebra

@testset "SVD Algorithms" begin

  @testset "Test Orthog" begin
    M1 = [1.0 2.0 5.0 4.0;
         1.0 1.0 1.0 1.0;
         0.0 0.5 0.5 1.0;
         0.0 1.0 1.0 2.0]
    orthog!(M1)
    @test norm(M1'*M1-Diagonal(ones(size(M1,1)))) < 1E-13

    M2 = rand(10,10)
    orthog!(M2)
    @test norm(M2'*M2-Diagonal(ones(size(M2,1)))) < 1E-13
  end

  @testset "Matrix With Zero Sing Val" begin
    M = [1.0 2.0 5.0 4.0;
         1.0 1.0 1.0 1.0;
         0.0 0.5 0.5 1.0;
         0.0 1.0 1.0 2.0]
    U,S,V = recursiveSVD(M)
    @test norm(U*Diagonal(S)*V'-M) < 1E-13
  end

  @testset "Real Matrix" begin
    M = rand(10,20)
    U,S,V = recursiveSVD(M)
    @test norm(U*Diagonal(S)*V'-M) < 1E-13

    M = rand(20,10)
    U,S,V = recursiveSVD(M)
    @test norm(U*Diagonal(S)*V'-M) < 1E-13
  end

  @testset "Cplx Matrix" begin
    M = rand(ComplexF64,10,15)
    U,S,V = recursiveSVD(M)
    @test norm(U*Diagonal(S)*V'-M) < 1E-13

    M = rand(ComplexF64,15,10)
    U,S,V = recursiveSVD(M)
    @test norm(U*Diagonal(S)*V'-M) < 1E-13
  end

end


text:
module StaticLinearMaps
using StaticArrays

import Base: +, -, *, \, /, ==, transpose
import StaticArrays: Size
export SAdjoint, SScaling, SComposition, domain, codomain, MProjection

if VERSION >= v"0.7.0-DEV.1415"
    const adjoint = Base.adjoint
else
    const adjoint = Base.ctranspose
end

abstract type SLinearMap{M,N} end

Size(::SLinearMap{M,N}) where {M,N} = (M,N)
Size(::Type{<:SLinearMap{M,N}}) where {M,N} = (M,N)

struct SScaling{N,T} <: SLinearMap{N,N}
    λ::T
end

*(a::SScaling, x::StaticVector) = a.λ*x
#*(x::StaticVector, a::SScaling) = x*a.λ
*(x::RowVector{T,S}, a::SScaling) where {T,S<:StaticVector} = x*a.λ

struct SAdjoint{N,M,T<:SLinearMap{M,N}} <: SLinearMap{N,M}
    f::T
end

adjoint(t::T) where {T<:SLinearMap{M,N}} where {M,N} = SAdjoint{N,M,T}(t)

#*(a::SAdjoint, x) = (x'*a.f)' # this fallback is too likely to produce stackoverflow errors
*(x, a::SAdjoint) = (a.f*x')'
*(x, a::T) where {T<:SLinearMap} = (a'*x')'

struct SComposition{M,N,P,S<:SLinearMap{M,N},T<:SLinearMap{N,P}} <: SLinearMap{M,P}
    f::S
    g::T
end

*(f::S, g::T) where {S<:SLinearMap{M,N},T<:SLinearMap{N,P}} where {M,N,P} = SComposition{M,N,P,S,T}(f,g)

*(c::SComposition, x) = c.f*(c.g*x)
*(y, c::SComposition) = (y*c.f)*c.g

struct SProjection{M,N,S<:SVector{M,Int}} <: SLinearMap{M,N}
    p::S
end

@generated function *(P::SProjection{M,N}, x::StaticVector{N,T}) where {M,N,T}
    expr = [:(x[P.p[$i]]) for i in 1:M]
    :(SVector{M,T}($(expr...)))
end

struct SInjection{M,N,S<:SVector{M,Int}} <: SLinearMap{M,N}
    p::S
end

@generated function *(P::SInjection{M,N}, x::StaticVector{N,T}) where {M,N,T}
    expr = [:(P.p[$i] > 0 ? x[P.p[$i]]' : zero(T)') for i in 1:M]
    :(SVector{M,T}($(expr...)))
end

function *(PT::SAdjoint{N,M,Q}, x::StaticVector{M,T}) where {N,M,T,Q<:SProjection}
    P = PT.f
    out = zero(MVector{N,T})
    for i in 1:length(P.p)
        out[P.p[i]] = x[i]
    end
    out
end

struct MProjection{M,N,S} <: SLinearMap{M,N}
    d::S
end

function *(P::MProjection{M,N}, x::StaticVector{N,T}) where {M,N,T}
    out = zero(MVector{M, T})
    for i in 1:length(P.d)
        out[P.d[i][2]] = x[P.d[i][1]]
    end
    out
end
function *(PT::SAdjoint{N,M,Q}, x::StaticVector{M,T}) where {N,M,T,Q<:MProjection}
    P = PT.f
    out = zero(MVector{N,T})
    for i in 1:length(P.d)
        out[P.d[i][1]] = x[P.d[i][2]]
    end
    out
end

end # module


text:
"""
    unsafe_full_panel_view(x, i, j)

NOTE: This function might be removed in the future. Use `full_panel_view` where possible.

Get an unsafe `Panel` view of panel (i,j) in matrix `x`. An entire panel is
returned even if the matrix does not cover this entire panel.

For performance reasons, a pointer is used instead of a proper object reference.
This means that the original PanelMatrix must be protected from garbage collection
while an unsafe `Panel` is in scope, or memory corruption may occur.

```
   GC.@preserve A let a = unsafe_full_panel_view(A, i, j)
       # do stuff with a, but don't let it escape the scope of the block.
   end
```
"""
@inline function unsafe_full_panel_view(x::PanelMatrix{T}, i::Integer, j::Integer) where T
    length = prod(x.panel_size)
    zeroth = ((i - 1) + (j - 1) * x.panel_stride) * length
    @boundscheck all(1 .<= (i, j) .* x.panel_size .<= x.size .+ x.pad_first .+ x.pad_last) || throw(panel_index_error)
    Panel{T}(pointer(x.data, zeroth+1), x.panel_size, x.panel_class, static.((0, 0)), static.((0, 0)))
end

# Copied from https://discourse.julialang.org/t/benchmarking-profiling-cache-use/7260
# and modified to fix https://github.com/JuliaLang/julia/issues/24164
@inline function prefetch(address::Ptr)
    Base.llvmcall(("declare void @llvm.prefetch(i8* , i32 , i32 , i32 )",
        """
        %2 = inttoptr i64 %0 to i8*
        call void @llvm.prefetch(i8 * %2, i32 0, i32 0, i32 1)
        ret void
        """), Nothing, Tuple{Int64}, Int64(address))
end


text:
using Julie
using Documenter: deploydocs, makedocs

makedocs(sitename = "Julie.jl", modules = [Julie], doctest = false)
deploydocs(repo = "github.com/ErikQQY/Julie.jl.git")

max_stars_repo_name: ErikQQY/Julie

max_issues_repo_name: ErikQQY/Julie

max_forks_repo_name: ErikQQY/Julie


text:
# Julia1_1.1.jl Why Learn Julia?
# running OS utilities using run() and read()
pwd()
run(`echo hello`)  # note (`) is a backtick`
files=String(read(`ls`))
files=split(files,"\n")
println(files)

# Compact
using LinearAlgebra
x = [1 2; 3 4]
y = x^2 .+ 3x .+ 1.1       # .+ corresponding component ops

eye=[1 0; 0 1.0] # unit matrix
z = eye/y                  # matrix inversion

# Notice compactness when working with matrices
# In this example several nested loops are needed in C
# the simple formula representation would also be lost.
# also notice implicit multiplication for constant

x=2.1+3.2im # complex floating-point number
# versatile polymorphic type system

# built in Int8-Int128, UInt8-UInt128, Float16-Float64
# other compound types made from these
# designed for linear algebra
# automatic conversions and promotions

# There is no meaningful concept of compile-time type,
# all types are actual at runtime.
# Only values, not variables have types.

# multicore parallel example, 1=head, 0=tail
# How many cores do we have>
Sys.CPU_THREADS

using Distributed
  heads = @distributed (+) for i=1:200000000
  Int(floor(rand()+0.5))
end
print(heads,"\n")

# Calling a c lib function
#using Distributed already called
t1=ccall((:clock, "libc"), UInt32,())
heads = @distributed (+) for i=1:200000000
  Int(floor(rand()+0.5))
end
print(heads,"\n")
t2=ccall((:clock, "libc"), UInt32,())
print((t2-t1)/1000000.," seconds\n")
# Don't need the C clock function since Julia has a
# @time macro and other timing operations.

# unicode characters in strings
s="abπ∑ef\n"
print(s)

# ok to use unicode characters in formulas
# generally Roman and Greek characters
ß=2π/3
µ=3+3.3im
Y=[1 2im; 3.3 6+4im]
Q = π*µ .+ ß*Y^2

# the C swap macro in Julia might be:
macro swap(a, b) # this will work for any type also mixed

    quote

        local temp=$(esc(a))

        $(esc(a))=$(esc(b))
        $(esc(b))=temp
    end
end
# Don't need swap macro
a,b=b,a

# metaprogramming machinery exposed to the user in Julia
# eval(), Expr(), parse(), Symbol(), typeof(), Regex()

# always check what Julia has to offer before rolling
# your own functions, macros and extensions.


text:
struct TrialFESpace{S} <: SingleFieldFESpace
  dirichlet_values::AbstractVector
  space::S
  function TrialFESpace(dirichlet_values::AbstractVector,space::SingleFieldFESpace)
    new{typeof(space)}(dirichlet_values,space)
  end
end

"""
"""
function TrialFESpace(space::SingleFieldFESpace)
  space
end

#function TrialFESpace(space::TrialFESpace)
#  space
#end

"""
"""
function TrialFESpace(space::SingleFieldFESpace,objects)
  dirichlet_values = compute_dirichlet_values_for_tags(space,objects)
  TrialFESpace(dirichlet_values,space)
end

"""
"""
function TrialFESpace!(dir_values::AbstractVector,space::SingleFieldFESpace,objects)
  dir_values_scratch = zero_dirichlet_values(space)
  dir_values = compute_dirichlet_values_for_tags!(dir_values,dir_values_scratch,space,objects)
  TrialFESpace(dir_values,space)
end

"""
"""
function TrialFESpace!(space::TrialFESpace,objects)
  dir_values = get_dirichlet_dof_values(space)
  dir_values_scratch = zero_dirichlet_values(space)
  dir_values = compute_dirichlet_values_for_tags!(dir_values,dir_values_scratch,space,objects)
  space
end

# Allow do-block syntax

function TrialFESpace(f::Function,space::SingleFieldFESpace)
  TrialFESpace(space,f)
end

function TrialFESpace!(f::Function,dir_values::AbstractVector,space::SingleFieldFESpace)
  TrialFESpace!(dir_values,space,f)
end

function TrialFESpace!(f::Function,space::TrialFESpace)
  TrialFESpace!(space,f)
end

# Remove Dirichlet from the given space

function HomogeneousTrialFESpace(U::SingleFieldFESpace)
  dirichlet_values = zero_dirichlet_values(U)
  TrialFESpace(dirichlet_values,U)
end

function HomogeneousTrialFESpace!(dirichlet_values::AbstractVector,U::SingleFieldFESpace)
  fill!(dirichlet_values,zero(eltype(dirichlet_values)))
  TrialFESpace(dirichlet_values,U)
end

# Genuine functions

get_dirichlet_dof_values(f::TrialFESpace) = f.dirichlet_values

# Delegated functions

get_free_dof_ids(f::TrialFESpace) = get_free_dof_ids(f.space)

zero_free_values(f::TrialFESpace) = zero_free_values(f.space)

get_triangulation(f::TrialFESpace) = get_triangulation(f.space)

get_dof_value_type(f::TrialFESpace) = get_dof_value_type(f.space)

get_vector_type(f::TrialFESpace) = get_vector_type(f.space)

get_cell_dof_ids(f::TrialFESpace) = get_cell_dof_ids(f.space)

get_fe_basis(f::TrialFESpace) = get_fe_basis(f.space)

get_trial_fe_basis(f::TrialFESpace) = get_trial_fe_basis(f.space)

get_fe_dof_basis(f::TrialFESpace) = get_fe_dof_basis(f.space)

ConstraintStyle(::Type{<:TrialFESpace{B}}) where B = ConstraintStyle(B)

get_cell_isconstrained(f::TrialFESpace) = get_cell_isconstrained(f.space)

get_cell_constraints(f::TrialFESpace) = get_cell_constraints(f.space)

get_dirichlet_dof_ids(f::TrialFESpace) = get_dirichlet_dof_ids(f.space)

get_cell_is_dirichlet(f::TrialFESpace) = get_cell_is_dirichlet(f.space)

zero_dirichlet_values(f::TrialFESpace) = zero_dirichlet_values(f.space)

num_dirichlet_tags(f::TrialFESpace) = num_dirichlet_tags(f.space)

get_dirichlet_dof_tag(f::TrialFESpace) = get_dirichlet_dof_tag(f.space)

scatter_free_and_dirichlet_values(f::TrialFESpace,fv,dv) = scatter_free_and_dirichlet_values(f.space,fv,dv)

gather_free_and_dirichlet_values(f::TrialFESpace,cv) = gather_free_and_dirichlet_values(f.space,cv)

gather_free_and_dirichlet_values!(fv,dv,f::TrialFESpace,cv) = gather_free_and_dirichlet_values!(fv,dv,f.space,cv)

gather_dirichlet_values(f::TrialFESpace,cv) = gather_dirichlet_values(f.space,cv)

gather_dirichlet_values!(dv,f::TrialFESpace,cv) = gather_dirichlet_values!(dv,f.space,cv)

gather_free_values(f::TrialFESpace,cv) = gather_free_values(f.space,cv)

gather_free_values!(fv,f::TrialFESpace,cv) = gather_free_values!(fv,f.space,cv)


text:
import Bukdu: Endpoint

Bukdu.config(:app, Endpoint,
        http= (:port => 80),
        https= (:port => 443),
        url= (:host => "example.com", :path => "/api"),
        static_url= (:host => "static.example.com"))

using Base.Test
@test Endpoint[:url] == Dict(:host => "example.com", :path => "/api")


text:
using NlmeData: dataset, datasets
using Test

@testset "NlmeData" begin
    ds = datasets()
    @test length(ds) == 48
    @test eltype(ds) == String
    oxboys = dataset(:Oxboys)
    @test length(oxboys) == 4
    @test keys(oxboys) == [:subject, :age, :height, :occasion]
    @test eltype(oxboys.age) == Float64
    @test_throws ArgumentError dataset("foo")
end


text:
function poli_sturm(w, d, e)

    n = size(d,1);
    pm = zeros(n+1,1);

    pm[1] = 1;
    pm[2] = w-d[1];

    for i = 3:n+1
        pm[i] = (w-d[i])*pm[i-1] - (e[i-1]^2) * pm[i-2];
    end
    return pm;
end


text:
using CoordinateTransformations
using Accessors

# P = AbstractParticle does not suffice
for P in [EGSParticle, IAEAParticle, CompressedIAEAParticle]
    @eval function (t::Translation)(p::$P)
        pos = t(position(p))
        set_position(p, pos)
    end
    @eval function (t::LinearMap)(p::$P)
        pos = t(position(p))
        dir = t(direction(p))
        set_position_direction(p, pos, dir)
    end
    @eval function (t::AffineMap)(p::$P)
        pos = t(position(p))
        dir = t.linear*direction(p)
        set_position_direction(p, pos, dir)
    end
end

"""
    propagate_z(p, z)

Propagate particle `p` such that `p.z=z`
"""
function propagate_z(p::P, z) where {P}
    if P <: AbstractParticle
        @argcheck has_z(P)
    end
    pos = position(p)
    dir = direction(p)
    z_to = z
    z_from = pos[3]
    t = (z_to - z_from) / dir[3]
    @set position(p) = pos + t * dir
end

export SetZ
"""

    SetZ(z)(p::EGSParticle)::EGSParticleZ

Set the `z` coordinate of `p`.
"""
struct SetZ
    z::Float32
end

(f::SetZ)(p::EGSParticle) = @set p.z = f.z


text:
const threshold = 5

function normalize_text(path::String)
    println("Tokenizing...")
    text = readstring(`java -cp stanford-ner.jar edu.stanford.nlp.process.PTBTokenizer $path -preserveLines`)

    println("Normalizing...")
    words = String[]
    countdict = Dict{String,Int}()
    for w in split(text," ")
        w = normalize_string(w, :NFKC)
        w = lowercase(w)
        w = replace(w, r"[0-9]", '0') # ex. 123.45 -> 000.00
        push!(words, w)
        if haskey(countdict, w)
            countdict[w] += 1
        else
            countdict[w] = 1
        end
    end
    println("Original vocabulary size: $(length(countdict))")

    unkdict = Dict{String,String}()
    for (k,v) in countdict
        v < threshold && (unkdict[k] = k)
    end
    c = length(countdict) - length(unkdict)
    println("Normalized (>= $threshold) vocabulary size: $c")

    words = map(words) do w
        haskey(unkdict,w) ? "UNKNOWN" : w
    end
    str = join(words, " ")
    open("$path.clean","w") do io
        println(io, str)
    end
end

normalize_text(ARGS[1])


text:
__precompile__(true)

"""
Module `Arrays.jl` -- Auxiliary machinery for vectors and matrices.
"""
module Arrays

using LinearAlgebra, SparseArrays
using ..Assertions: @assert, activate_assertions
# activate assertions by default
activate_assertions(Arrays)
using ..LazySets: _geq

include("matrix_operations.jl")
include("vector_operations.jl")
include("matrix_vector_operations.jl")
include("SingleEntryVector.jl")

end # module


text:
export getSensTMatVecbatch

function getSensTMatVecbatch(v::Vector,srcind::Vector,recind::Vector,pFor::EikonalInvParam)

   # extract pointers
    Q     = pFor.Sources
    P     = pFor.Receivers
    nsrc = size(Q,2); nrec = size(P,2)

    n = length(srcind);
    nt = size(Q,1);

    # allocate space for result
    JTv  = zeros(nt)
	JTvi = zeros(nt)
	pEiks = pFor.eikonalParams;

	#if length(v)!=nsrc*nrec
	#	warn("wrong size of vector v in EikonalInv::getSensTMatVec.jl");
	#end

    for i=1:n
		tdum = zeros(nrec);
		tdum[recind[i]] = 1.0;
		#print(v[(srcind[i]-1)*nrec + recind[i]]*tdum);
		t = P*(v[(srcind[i]-1)*nrec + recind[i]]*tdum); #P*v[(recind[i]-1)*nrec + 1 : recind[i]*nrec];
		getSensTMatVecEik(t,JTvi,pEiks[srcind[i]]);
        JTv += JTvi;
    end
    #print(size(JTv));
    return JTv;
end


text:
# Functions to compute linear regression statistics.
using GLM, Distributions

function lin_reg( df::DataFrame, Y::Symbol, X::Symbol )
  println("X: ",X,"  Y: ",Y)
  @assert Y in names(df)
  @assert X in names(df)
  d = TDist(size(df)[1]-2)
  global ddf
  ddf = df
  #lreg = lm(@formula( evolvable_count ~ complexity), ddf )
  cmd="lm(@formula($(String(Y)) ~ $(String(X))), ddf )"
  println("cmd: ",cmd)
  lreg = eval(Meta.parse(cmd))
  slope = coef(lreg)[2]
  tempdf = DataFrame( X=>[0.0])
  intercept = predict(lreg,tempdf)[1]
  stderr = stderror(lreg)[2]
  t = slope/stderr
  pvalue = slope > 0.0 ? ccdf(d,t) : cdf(d,t)
  qminus = quantile(d,0.025)
  qplus = quantile(d,0.975)
  lowerci = qplus*stderr + slope
  upperci = qminus*stderr + slope
  numinputs = :evo_count in names(df) ? df.numinputs[1] : 3
  ( Y, X, numinputs, slope, intercept, stderr, t, pvalue, lowerci, upperci )
end

function run_lin_reg( df_field_pair_list::Vector{Tuple{DataFrame,Symbol,Symbol}};
    csvfile::String="" )
  lrdf = DataFrame()
  lrdf.Y = Symbol[]
  lrdf.X = Symbol[]
  lrdf.numinputs = Int64[]
  lrdf.slope = Float64[]
  lrdf.intercept = Float64[]
  lrdf.stderr = Float64[]
  lrdf.t = Float64[]
  lrdf.pvalue = Float64[]
  lrdf.lower_ci = Float64[]
  lrdf.upper_ci = Float64[]
  for p in df_field_pair_list
    lrdf_row = lin_reg( p[1], p[2], p[3] )
    push!(lrdf, lrdf_row)
  end
  if length(csvfile) > 0
    hostname = chomp(open("/etc/hostname") do f read(f,String) end)
    open( csvfile, "w" ) do f
      println(f,"# date and time: ",Dates.now())
      println(f,"# host: ",hostname," with ",nprocs()-1,"  processes: " )
      println(f,"# funcs: ", Main.CGP.default_funcs(lrdf.numinputs[end]))
      CSV.write( f, lrdf, append=true, writeheader=true )
    end
  end
  lrdf
end


text:
using DelftBlue
using Test

@testset "DelftBlue.jl" begin
    # Write your tests here.
end


text:
using Gnuplot, ColorSchemes
x = -1:0.05:1
y = -1.5:0.05:2
egg(x,y) = x^2 + y^2/(1.4 + y/5)^2
segg = [egg(x,y) for x in x, y in y]
@gsp x y segg "w l lc palette" palette(:cool; rev= true) "set view map"
@gsp :- "set contour base;set key off" "set auto fix"
@gsp :- "set cntrparam levels incremental  0,0.01,1" "unset surface"
@gsp :- xrange = (-1.2,1.2) yrange = (-1.5,2) cbrange =(0,1)
@gsp :- xlab = "x" ylab = "y" "set size ratio -1"
#@gsp :- "unset xtics; unset ytics; unset border" # hide
save(term = "pngcairo size 600,800", output = "egg.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
function get_trajopt_data_semidirect(mechanism::Mechanism,env::Environment,Δt::Real,N::Int;
                                    relax_comp=false)
    vs = VariableSelector()

    for n = 1:N
        add_var!(vs, Symbol("q", n), num_positions(mechanism))
        add_var!(vs, Symbol("v", n), num_velocities(mechanism))
        if n < N
            add_var!(vs, Symbol("u", n), num_velocities(mechanism))
            add_var!(vs, Symbol("h", n), 1)
            for i = 1:length(env.contacts)
                if relax_comp
                    add_var!(vs, Symbol("slack", i, "_", n), 1)
                end
                add_var!(vs, Symbol("c_n", i, "_", n), 1)
            end
        end
    end

    cs = ConstraintSelector()

    for n = 1:N-1
        add_eq!(cs, Symbol("kin", n), num_positions(mechanism))
        add_eq!(cs, Symbol("dyn", n), num_velocities(mechanism))
        add_ineq!(cs, Symbol("h_pos", n), 1)
        for i = 1:length(env.contacts)
            if relax_comp
                add_ineq!(cs, Symbol("slack_pos", i, "_", n), 1)
                add_ineq!(cs, Symbol("ϕ_c_n_comp", i, "_", n), 1)
            else
                add_eq!(cs, Symbol("ϕ_c_n_comp", i, "_", n), 1)
            end
            add_ineq!(cs, Symbol("c_n_pos", i, "_", n), 1)
            add_ineq!(cs, Symbol("ϕ_pos", i, "_", n), 1)
        end
    end

    state_cache = [StateCache(mechanism) for n = 1:2]
    envj_cache = [EnvironmentJacobianCache(env) for n = 1:2]

    generate_solver_fn = :generate_solver_fn_trajopt_semidirect
    extract_sol = :extract_sol_trajopt_semidirect

    fric_vs = []
    fric_cs = []
    fric_options = []

    f_options = Dict{String, Any}()
    f_options["num_fosteps"] = 0
    f_options["num_sosteps"] = 20
    f_options["c"] = 1.
    f_options["c_fos"] = 1.
    f_options["c_sos"] = 1.

    for n = 1:N-1
        f_vs = VariableSelector()
        for i = 1:length(env.contacts)
            β_dim = size(env.contacts[i].obstacle.basis,2)
            add_var!(f_vs, Symbol("β", i), β_dim)
        end

        f_cs = ConstraintSelector()
        for i = 1:length(env.contacts)
            β_dim = size(env.contacts[i].obstacle.basis,2)
            add_ineq!(f_cs, Symbol("β_pos", i), β_dim)
            add_ineq!(f_cs, Symbol("fric_cone", i), 1)
        end

        push!(fric_vs, f_vs)
        push!(fric_cs, f_cs)
        push!(fric_options, f_options)
    end

    fric_models = []
    for n = 1:N-1
        fric_P = sparse(ones(4,4))
        fric_q = ones(4)
        fric_A = sparse(ones(5,4))
        fric_l = ones(5)
        fric_u = ones(5)
        fric_model = OSQP.Model()
        OSQP.setup!(fric_model, P=fric_P, q=fric_q, A=fric_A, l=fric_l, u=fric_u, verbose=0)

        push!(fric_models, fric_model)
    end

    sim_data = SimData(mechanism,env,
                       state_cache,envj_cache,
                       Δt,vs,cs,generate_solver_fn,extract_sol,
                       [],[],[],
                       [],[],[],
                       fric_vs,fric_cs,fric_options,
                       N,[],[],fric_models,[])

    sim_data
end

function extract_sol_trajopt_semidirect(sim_data::SimData, xopt::AbstractArray{T}) where T
    N = sim_data.N
    vs = sim_data.vs
    env = sim_data.env
    relax_comp = haskey(vs.vars, :slack1)

    qtraj = Array{Array{Float64,1},1}(undef, 0)
    vtraj = Array{Array{Float64,1},1}(undef, 0)
    utraj = Array{Array{Float64,1},1}(undef, 0)
    htraj = Array{Float64,1}(undef, 0)
    contact_traj = Array{Array{Float64,1},1}(undef, 0)
    slack_traj = Array{Array{Float64,1},1}(undef, 0)
    x = MechanismState(sim_data.mechanism)
    for n = 1:N
        set_configuration!(x, vs(xopt, Symbol("q", n)))
        normalize_configuration!(x)
        q = configuration(x)
        push!(qtraj, q)
        push!(vtraj, vs(xopt, Symbol("v", n)))
        if n < N
            push!(utraj, vs(xopt, Symbol("u", n)))
            push!(htraj, vs(xopt, Symbol("h", n))[1])
            if relax_comp
                push!(slack_traj, vs(xopt, Symbol("slack", n)))
            end
            contact_sol = []
            for i = 1:length(env.contacts)
                contact_sol = vcat(contact_sol,
                                   vs(xopt, Symbol("c_n", i, "_", n)))
            end
            push!(contact_traj, contact_sol)
        end
    end

    # some other useful vectors
    ttraj = vcat(0., cumsum(htraj)...)
    qv_mat = vcat(hcat(qtraj...),hcat(vtraj...))

    qtraj, vtraj, utraj, htraj, contact_traj, slack_traj, ttraj, qv_mat, xopt
end

function generate_solver_fn_trajopt_semidirect(sim_data::SimData)
    N = sim_data.N
    vs = sim_data.vs
    cs = sim_data.cs

    num_contacts = length(sim_data.env.contacts)
    num_vel = num_velocities(sim_data.mechanism)
    world = root_body(sim_data.mechanism)
    world_frame = default_frame(world)
    relax_comp = haskey(vs.vars, :slack1)

    function eval_obj(x::AbstractArray{T}) where T
        f = 0.

        if relax_comp
            for n = 1:N-1
                for i = 1:num_contacts
                    slack = vs(x, Symbol("slack", i, "_", n))
                    f += sum(slack)
                end
            end
        end

        # extra user-defined objective
        for i = 1:length(sim_data.obj_fns)
            obj_name, obj_fn = sim_data.obj_fns[i]
            f += obj_fn(x)
        end

        f
    end

    function eval_cons(x::AbstractArray{T}) where T
        g = Vector{T}(undef, cs.num_eqs + cs.num_ineqs) # TODO preallocate

#         @threads for n = 1:(N-1)
        for n = 1:(N-1)
            q0 = vs(x, Symbol("q", n))
            v0 = vs(x, Symbol("v", n))
            u0 = vs(x, Symbol("u", n))
            h = vs(x, Symbol("h", n))[1]
            qnext = vs(x, Symbol("q", n+1))
            vnext = vs(x, Symbol("v", n+1))
            if relax_comp
                slack = vs(x, Symbol("slack", n))
            end

            x0 = sim_data.state_cache[1][T]
            xn = sim_data.state_cache[2][T]
            envj = sim_data.envj_cache[2][T]

            set_configuration!(x0, q0)
            set_velocity!(x0, v0)
            setdirty!(x0)
            set_configuration!(xn, qnext)
            set_velocity!(xn, vnext)
            setdirty!(xn)

            H = mass_matrix(x0)
            Hi = inv(H)
            config_derivative = configuration_derivative(xn) # TODO preallocate
            dyn_bias = dynamics_bias(xn) # TODO preallocate

            contact_jacobian!(envj, xn)
            contact_bias = Vector{T}(undef, num_vel)
            if (num_contacts > 0)
                # compute friction forces
                # contact_friction_τ_direct!(contact_bias, sim_data, h, Hi, envj, dyn_bias, u0, v0, x, n)
                contact_friction_τ_direct_osqp!(contact_bias, sim_data, h, Hi, envj, dyn_bias, u0, v0, x, n)
            end

            g[cs(Symbol("kin", n))] .= qnext .- q0 .- h .* config_derivative
            g[cs(Symbol("dyn", n))] .= H * (vnext - v0) .- h .* (u0 .- dyn_bias .- contact_bias)
            g[cs(Symbol("h_pos", n))] .= [-h]
            for i = 1:num_contacts
                c_n = vs(x, Symbol("c_n", i, "_", n))
                g[cs(Symbol("c_n_pos", i, "_", n))] .= -c_n
                g[cs(Symbol("ϕ_pos", i, "_", n))] .= -envj.contact_jacobians[i].ϕ
                if !relax_comp
                    g[cs(Symbol("ϕ_c_n_comp", i, "_", n))] .= envj.contact_jacobians[i].ϕ .* c_n
                else
                    slack = vs(x, Symbol("slack", i, "_", n))
                    g[cs(Symbol("slack_pos", i, "_", n))] .= -slack
                    g[cs(Symbol("ϕ_c_n_comp", i, "_", n))] .= envj.contact_jacobians[i].ϕ .* c_n .- slack
                end
            end
        end

        # extra user-defined constraints
        for i = 1:length(sim_data.con_fns)
            con_name, con_fn = sim_data.con_fns[i]
            g[cs(con_name)] .= con_fn(x)
        end

        g
    end

    generate_autodiff_solver_fn(eval_obj,eval_cons,cs.eqs,cs.ineqs,vs.num_vars)
end


text:
function gradient(H::LocalHamiltonian, Ψρs::InfiniteCMPSData, HL = nothing, HR = nothing;
                    kwargs...)
    Ψ, ρL, ρR = Ψρs
    if isnothing(HL)
        HL, = leftenv(H, Ψρs; kwargs...)
    end
    if isnothing(HR)
        HR, = rightenv(H, Ψρs; kwargs...)
    end

    Q = Ψ.Q
    Rs = Ψ.Rs

    # gradQ = ∑(coeff * localgradientQ)  +  HL*ρR + ρL*HR
    gradQ = zero(Q)
    for (coeff, op) in zip(coefficients(H.h), operators(H.h))
        if coeff isa Number
            axpy!(coeff, localgradientQ(op, Q, Rs, ρL, ρR), gradQ)
        else
            mul!(gradQ, coeff, localgradientQ(op, Q, Rs, ρL, ρR), 1, 1)
        end
    end
    mul!(gradQ, HL, ρR, 1, 1)
    mul!(gradQ, ρL, HR, 1, 1)

    # gradR = ∑(coeff * localgradientR)  +  HL*R*ρR + ρL*R*HR
    gradRs = zero.(Rs)
    for (coeff, op) in zip(coefficients(H.h), operators(H.h))
        if coeff isa Number
            axpy!.(coeff, localgradientRs(op, Q, Rs, ρL, ρR), gradRs)
            if op isa ContainsDifferentiatedCreation && !(Q isa Constant)
                grad∂Rs = localgradient∂Rs(op, Q, Rs, ρL, ρR)
                axpy!.(-coeff, ∂.(grad∂Rs), gradRs)
            end
        else
            mul!.(gradRs, (coeff,), localgradientRs(op, Q, Rs, ρL, ρR), 1, 1)
            if op isa ContainsDifferentiatedCreation && !(Q isa Constant)
                grad∂Rs = localgradient∂Rs(op, Q, Rs, ρL, ρR)
                mul!.(gradRs, (-coeff,), ∂.(grad∂Rs), 1, 1)
            end
        end
    end
    mul!.(gradRs, (HL,), Rs .* (ρR,), 1, 1)
    mul!.(gradRs, (ρL,), Rs .* (HR,), 1, 1)

    return gradQ, gradRs
end


text:
using Sundials

## f routine. Compute function f(t,y).

function f(t, y, ydot)
    ydot[1] = -0.04*y[1] + 1.0e4*y[2]*y[3]
    ydot[3] = 3.0e7*y[2]*y[2]
    ydot[2] = -ydot[1] - ydot[3]
    return Sundials.CV_SUCCESS
end

t = [0.0; 4 * exp10.(range(-1., stop=7., length=9))]
y0 = [1.0, 0.0, 0.0]
res = Sundials.cvode(f, y0, t)


text:
distrname(X::ElementaryRandVar) = string(typeof(X).name)
distrname(X::SymbolicRandVar) = "CompositeRandVar"

"For SymbolcRandVar just retrn ndims"
get_params(X::SymbolicRandVar) = "(ndims=$(ndims(X)))"

"Parameters in form p1=value, e.g. (dim=1, lb=0.0, ub=1.0)"
function get_params(X::ElementaryRandVar)
  params = ["$field=$(getfield(X,field))" for field in fieldnames(X)]
  string("(",join(params, ", "),")")
end

"Show a RandVar with its parameters e.g.,(dim=1, lb=0.0, ub=1.0)"
function show{T}(io::IO, X::SymbolicRandVar{T})
  print(io, "$(distrname(X)){$T}$(get_params(X))")
end

show{T}(io::IO, X::ConstantRandVar{T}) = print(io, X.val)


text:
# Let's make the stroud stuff @generated.

function stroud_eval1d(deg, coeffs, result, qpts)
    nqp1d = size(qpts, 1)
    for i=1:nqp1d
        result[i] = 0.0
        xi = qpts[i, 1]
        s = 1-xi
        r = xi / s
        w = s^deg
        for alpha=0:deg
            result[i] += w * coeffs[alpha+1]
            w *= r * (deg-alpha) / (alpha+1.0)
        end
    end
    return
end

function stroud_eval2d(deg, coeffs, result, qpts)
    nqp1d = size(qpts, 1)
    result[:] = 0.0
    tmp = zeros(nqp1d, deg+1)
    for i2=1:nqp1d
        xi = qpts[i2, 2]
        s = 1 - xi
        r = xi / s
        for alpha1=0:deg
            w = s^(deg-alpha1)
            aind = polydim(2,deg) - polydim(2, deg-alpha1) + 1
            for alpha2=0:deg-alpha1
                tmp[i2, alpha1+1] += w * coeffs[aind+alpha2]
                w = w * r * (deg-alpha1-alpha2) / (1.0 + alpha2)
            end
        end
    end
    for i1=1:nqp1d
        xi = qpts[i1, 1]
        s = 1-xi
        r = xi / s
        w = s^deg
        for alpha1=0:deg
            for i2=1:nqp1d
                result[(i1-1)*nqp1d + i2] += w*tmp[i2, alpha1+1]
            end
            w = w * r * (deg-alpha1) / (1. + alpha1)
        end
    end

    return
end

val2val{d}(::Type{Val{d}}) = d

@generated function stroud_eval{sdim}(::Type{Val{sdim}}, deg, coeffs, result, qpts)
    return stroud_eval_impl(Val{sdim}, deg, coeffs, result, qpts)
end

function make_stroud_eval_constructor(sd)
    # return sdim = 1 code as a special case
    if sd == 1
        function makeit(deg, qpts)
            function doit(coeffs, result)
                nqp1d = size(qpts, 1)
                for i=1:nqp1d
                    result[i] = 0.0
                    xi = qpts[i, 1]
                    s = 1-xi
                    r = xi / s
                    w = s^deg
                    for alpha=0:deg
                        result[i] += w * coeffs[alpha+1]
                        w *= r * (deg-alpha) / (alpha+1.0)
                    end
                end
                return
            end
            return doit
        end
        return makeit
    else # sd >= 2
        # first we have to generate the function body.
        ai = [symbol("alpha",i) for i=1:sd]
        ii = [symbol("i",i) for i=1:sd]

        degminalphas = Dict()
        degminalphas[0] = :deg
        for i=1:sd
            degminalphas[i] = :($(degminalphas[i-1]) - $(ai[i]))
        end

        # work on the first stage, which is different since it
        # reads from the input

        tmp_inds = vcat([:($(ii[sd]))],
                        [:($(ai[i])+1) for i=(sd-1):-1:1],
                        [:1])
        tmp_ref = Expr(:ref,
                       :tmp)

        for ti in tmp_inds
            push!(tmp_ref.args, ti)
        end

        first_stage = Expr(:for,
                           :($(ai[sd])=0:$(degminalphas[sd-1])),
                           Expr(:block,
                                :($(tmp_ref) += w*coeffs[aind+$(ai[sd])]),
                                :(w = w * r * $(degminalphas[sd]) / (1.0 + $(ai[sd])))))

        aind = :(binomial($(sd)+deg, $(sd)) - binomial($(sd)+deg-alpha1, $(sd)))
        for i=(sd-1):-1:2
            aind = :($(aind) + binomial($(i)+ $(degminalphas[sd-i]),$(i))
                     - binomial($(i)+ $(degminalphas[sd-i+1]), $(i)))
        end
        aind = :($(aind)+1)
        first_stage = Expr(:block,
                           :(aind = $(aind)),
                           :(w=s^($(degminalphas[sd-1]))),
                           first_stage)

        for i=(sd-1):-1:1
            first_stage = Expr(:for,
                               :($(ai[i])=0:$(degminalphas[i-1])),
                               first_stage)
        end

        first_stage = Expr(:for,
                           :($(ii[sd])=1:nqp1d),
                           Expr(:block,
                                :(xi=qpts[$(ii[sd]), $(sd)]),
                                :(s=1-xi),
                                :(r=xi/s),
                                first_stage))

        # middle stages:
        # on stage k, we have k-1 innermost i loops over last i variables
        # tucked inside of (sd-k+1) alpha loops
        # tucked inside of a loop over ik

        middle_stages = []

        for k=2:(sd-1)
            # start wth write location for tmp
            iinds = [:($(ii[sd-i+1])) for i=1:k]
            ainds = [:($(ai[i])+1) for i=(sd-k):-1:1]
            tmp_inds = vcat(iinds, ainds, [:($k)])
            tmp_write_ref = Expr(:ref,
                                 :tmp)

            for ti in tmp_inds
                push!(tmp_write_ref.args, ti)
            end

            # read location for tmp
            iinds = [:($(ii[sd-i+1])) for i=1:(k-1)]
            ainds = [:($(ai[i])+1) for i=(sd-k+1):-1:1]
            tmp_inds = vcat(iinds, ainds, [:($(k-1))])
            tmp_read_ref = Expr(:ref,
                                :tmp)

            for ti in tmp_inds
                push!(tmp_read_ref.args, ti)
            end

            ms = :($(tmp_write_ref) += w * $(tmp_read_ref))

            # wrap ms into the i loops
            for i=1:(k-1)
                ms = Expr(:for,
                          :($(ii[sd-i+1])=1:nqp1d),
                          ms)
            end

            # the innermost alpha is sd - k + 1
            j = sd - k
            ms = Expr(:for,
                      :($(ai[j+1])=0:$(degminalphas[j])),
                      Expr(:block,
                           ms,
                           :(w=w*r*($(degminalphas[j+1])) / (1.0+$(ai[j+1])))))

            ms = Expr(:block,
                      :(w=s^($(degminalphas[j]))),
                      ms)

            for i=k-1:-1:1
                ms = Expr(:for,
                          :($(ai[i])=0: $(degminalphas[i-1])),
                          ms)
            end

            ms = Expr(:for,
                      :($(ii[sd-k+1])=1:nqp1d),
                      Expr(:block,
                           :(xi=qpts[$(ii[sd-k+1]), $(sd-k+1)]),
                           :(s=1-xi),
                           :(r=xi/s),
                           ms)
                      )

            push!(middle_stages, ms)

        end

        # last stage is different because we're writing into the result
        # we have an outer loop over i1
        write_idx = :($(ii[sd]))
        for i=(sd-1):-1:1
            write_idx = :(($(ii[i])-1)*nqp1d^($(sd-i)) + $(write_idx))
        end

        tmp_idx = Expr(:ref,
                       :tmp)
        for i=sd:(-1):2
            push!(tmp_idx.args, :($(ii[i])))
        end
        push!(tmp_idx.args, :(alpha1+1))
        push!(tmp_idx.args, :($(sd-1)))

        last_stage = Expr(:+=,
                          Expr(:ref,
                               :result,
                               write_idx),
                          :(w * $(tmp_idx)))

        for i=sd:-1:2
            last_stage = Expr(:for,
                              :($(ii[i])=1:nqp1d),
                              last_stage)
        end

        last_stage = Expr(:for,
                          :(i1=1:nqp1d),
                          Expr(:block,
                               :(xi=qpts[i1,1]),
                               :(s=1-xi),
                               :(r=xi/s),
                               :(w=s^deg),
                               Expr(:for,
                                    :(alpha1=0:deg),
                                    Expr(:block,
                                         last_stage,
                                         :(w=w*r*(deg-alpha1)/(1.+alpha1))))))

        make_tmp = Expr(:call, :zeros)
        push!(make_tmp.args, :nqp1d)
        for i=1:sd-1
           push!(make_tmp.args, :(deg+1))
        end
        push!(make_tmp.args, sd-1)

        doit_func = Expr(:function,
                         :(doit(coeffs, result)))

        doit_func_body = Expr(:block)

        push!(doit_func_body.args,
              :(for i=1:length(result) result[i] = 0.0 end))
        push!(doit_func_body.args,
              first_stage)
        for ms in middle_stages
            push!(doit_func_body.args, ms)
        end
        push!(doit_func_body.args, last_stage)

        push!(doit_func.args, doit_func_body)

        func_impl = quote
            function makeit(deg, qpts)
                nqp1d = size(qpts, 1)
                tmp = $(make_tmp)
                $(doit_func)
                return doit
            end
        end

        println(func_impl)

        return eval(func_impl)
end

end


text:
using ClickHouse: Column, chwrite, chread,
         read_col, VarUInt, parse_typestring, result_type
using Dates
using CategoricalArrays
using UUIDs
import Sockets
using Sockets: IPv4, IPv6
using DecFP

@testset "Parse type" begin
    r = parse_typestring("Int32")
    @test r.name == :Int32
    @test_throws ErrorException parse_typestring("KKKK")

    r = parse_typestring("   String  ")
    @test r.name == :String
    @test result_type(r) == Vector{String}

    r = parse_typestring("   Enum8('a' = 10, 'b'=1, 'addd' = 45)  ")

    @test r.name == :Enum8
    @test length(r.args) == 3
    @test r.args[1] == "'a' = 10"
    @test r.args[2] == "'b'=1"
    @test r.args[3] == "'addd' = 45"
    @test result_type(r) == CategoricalVector{String}

    r = parse_typestring(" FixedString(4)")
    @test r.name == :FixedString
    @test r.args[1] == "4"
    r = parse_typestring(" FixedString(44)")
    @test r.name == :FixedString
    @test r.args[1] == "44"
    @test result_type(r) == Vector{String}

    r = parse_typestring("Tuple(Int64, String)")
    @test r.name == :Tuple
    @test r.args[1].name == :Int64
    @test r.args[2].name == :String
    @test result_type(r) == Vector{Tuple{Int64, String}}

    r = parse_typestring("Tuple(Enum16('a' = 10), Tuple(Int32, Float32))")
    @test r.name == :Tuple
    @test r.args[1].name == :Enum16
    @test r.args[1].args[1] == "'a' = 10"
    @test r.args[2].name == :Tuple
    @test r.args[2].args[1].name == :Int32
    @test r.args[2].args[2].name == :Float32
    @test result_type(r) == Vector{
        Tuple{
            CategoricalValue{String},
            Tuple{Int32, Float32}
            }
        }

    r = parse_typestring("LowCardinality(String)")
    @test result_type(r) == CategoricalVector{String}

    r = parse_typestring("Array(Array(Nullable(Int32)))")
    @test result_type(r) == Vector{
        Vector{
            Vector{Union{Missing, Int32}}
        }
    }

end

@testset "Int columns" begin

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = Sockets.IPv4.(rand(UInt32, nrows))
    column = Column("test", "IPv4", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test res == column

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = Sockets.IPv6.(rand(UInt128, nrows))
    column = Column("test", "IPv6", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test res == column

end

@testset "IP columns" begin

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = rand(Int64, nrows)
    column = Column("test", "Int64", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test res == column

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = rand(Int32, nrows)
    column = Column("test", "Int64", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test res == column

end

@testset "String columns" begin

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = string.(rand(Int64, nrows))
    column = Column("test", "String", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test res == column

end

@testset "Fixed String columns" begin

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = string.(rand(["aaaa", "bbbb", "cccc", "dddd"], nrows))
    column = Column("test", "FixedString(4)", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test res == column

    sock = ClickHouseSock(PipeBuffer())
    data = string.(rand(["aaaaa", "bbbb", "cccc", "dddd"], nrows))
    column = Column("test", "FixedString(4)", data)
    @test_throws ErrorException chwrite(sock, column)

    sock = ClickHouseSock(PipeBuffer())
    data = string.(rand(["aaaa", "bbbb", "cccc", "dddd"], nrows))
    column = Column("test", "FixedString(5)", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test all(s->length(s)==5, res.data)
    @test res.data == data.* " "
end

@testset "Date columns" begin

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = Date.(rand(2010:2020, nrows), rand(1:12, nrows), rand(1:20, nrows))
    column = Column("test", "Date", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test res == column

end

@testset "DateTime columns" begin

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = DateTime.(rand(2010:2020, nrows), rand(1:12, nrows), rand(1:20, nrows),
    rand(0:23, nrows), rand(0:59, nrows), rand(0:59, nrows))
    column = Column("test", "DateTime", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test res == column

end

@testset "DateTime64 columns" begin

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = DateTime.(rand(2010:2020, nrows), rand(1:12, nrows), rand(1:20, nrows),
    rand(0:23, nrows), rand(0:59, nrows), rand(0:59, nrows))
    column = Column("test", "DateTime64(0)", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test res == column

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = DateTime.(rand(2010:2020, nrows), rand(1:12, nrows), rand(1:20, nrows),
    rand(0:23, nrows), rand(0:59, nrows), rand(0:59, nrows))
    column = Column("test", "DateTime64(2)", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test res == column

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = DateTime.(rand(2010:2020, nrows), rand(1:12, nrows), rand(1:20, nrows),
    rand(0:23, nrows), rand(0:59, nrows), rand(0:59, nrows),
    rand(1:999))
    column = Column("test", "DateTime64(3)", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test res == column

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = DateTime.(rand(2010:2020, nrows), rand(1:12, nrows), rand(1:20, nrows),
    rand(0:23, nrows), rand(0:59, nrows), rand(0:59, nrows),
    rand(1:999))
    column = Column("test", "DateTime64(6)", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test res == column

end

@testset "Enum columns" begin

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = rand(["a","b","c"], nrows)
    column = Column("test", "Enum8('a'=1,'b'=3,'c'=10)", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test res == column

end

@testset "Enum columns categorial in" begin

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = CategoricalVector(rand(["a","b","c"], nrows))
    column = Column("test", "Enum8('a'=1,'b'=3,'c'=10)", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))

    @test res == column

end

@testset "UUID columns" begin

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = Vector{UUID}(undef, nrows)
    data .= uuid4()
    column = Column("test", "UUID", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test res == column

end

@testset "Tuple columns" begin

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = tuple.(rand(Int64, nrows), rand(Int8, nrows))
    column = Column("test", "Tuple(Int64, Int8)", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test res == column

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = tuple.(rand(Int64, nrows), string.(rand(Int8, nrows)))
    column = Column("test", "Tuple(Int64, String)", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test res == column

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = tuple.(
        rand(["aa", "bb", "ccc"], nrows),
        tuple.(rand(Int16, nrows))
        )
    column = Column("test", "Tuple(Enum16('aa' = 1, 'bb' = 2, 'ccc' = 10), Tuple(Int16))", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test res == column

end

@testset "Nullable columns" begin

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = rand(Int64, nrows)
    data = convert(Vector{Union{Int64, Missing}}, data)
    data[rand(1:nrows, 20)] .= missing
    column = Column("test", "Nullable(Int64)", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test all(zip(data, res.data)) do t
        (a, b) = t
        return (ismissing(a) && ismissing(b)) ||
        (a == b)
    end

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = rand(Float64, nrows)
    data = convert(Vector{Union{Float64, Missing}}, data)
    data[rand(1:nrows, 20)] .= missing
    column = Column("test", "Nullable(Float64)", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test all(zip(data, res.data)) do t
        (a, b) = t
        return (ismissing(a) && ismissing(b)) ||
        (a ≈ b)
    end

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = string.(rand(Int64, nrows))
    data = convert(Vector{Union{String, Missing}}, data)
    data[rand(1:nrows, 20)] .= missing
    column = Column("test", "Nullable(String)", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test all(zip(data, res.data)) do t
        (a, b) = t
        return (ismissing(a) && ismissing(b)) ||
        (a == b)
    end

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = CategoricalVector(rand(["a","b","c", missing], nrows))

    data[rand(1:nrows, 20)] .= missing

    column = Column("test", "Nullable(Enum8('a'=1,'b'=3,'c'=10))", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))

    @test all(zip(data, res.data)) do t
        (a, b) = t
        return (ismissing(a) && ismissing(b)) ||
        (a == b)
    end

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = Date.(rand(2010:2020, nrows), rand(1:12, nrows), rand(1:20, nrows))
    data = convert(Vector{Union{Date, Missing}}, data)
    data[rand(1:nrows, 20)] .= missing
    column = Column("test", "Nullable(Date)", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test all(zip(data, res.data)) do t
        (a, b) = t
        return (ismissing(a) && ismissing(b)) ||
        (a == b)
    end

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = DateTime.(rand(2010:2020, nrows), rand(1:12, nrows), rand(1:20, nrows),
    rand(0:23, nrows), rand(0:59, nrows), rand(0:59, nrows))
    data = convert(Vector{Union{DateTime, Missing}}, data)
    data[rand(1:nrows, 20)] .= missing
    column = Column("test", "Nullable(DateTime)", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test all(zip(data, res.data)) do t
        (a, b) = t
        return (ismissing(a) && ismissing(b)) ||
        (a == b)
    end
end

@testset "LowCardinality columns" begin

    sock = ClickHouseSock(PipeBuffer())
    nrows = 10
    data = rand(1:10, nrows)
    column = Column("test", "LowCardinality(Int64)", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))

    @test res.data == data

    sock = ClickHouseSock(PipeBuffer())
    data = rand(1:10, nrows)
    data = convert(Vector{Union{Int64, Missing}}, data)
    data[rand(1:nrows, 5)] .= missing
    column = Column("test", "LowCardinality(Nullable(Int64))", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))

    @test all(zip(data, res.data)) do t
        (a, b) = t
        return (ismissing(a) && ismissing(b)) ||
        (a == b)
    end

    sock = ClickHouseSock(PipeBuffer())
    data = rand(["a", "b", "c"], nrows)
    column = Column("test", "LowCardinality(String)", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))

    @test res.data == data

    sock = ClickHouseSock(PipeBuffer())
    data = CategoricalVector(rand(["a","b","c"], nrows))

    column = Column("test", "LowCardinality(Enum8('a'=1,'b'=3,'c'=10))", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test res.data == data

end

@testset "Array collumns" begin

    nrows = 1000
    sock = ClickHouseSock(PipeBuffer())
    data = rand([[[1,2],[3]], [[3],[4,5]], [[6],[7,8]]], nrows)
    column = Column("test", "Array(Array(Int64))", data)
    chwrite(sock, column)

    res = read_col(sock, VarUInt(nrows))
    @test res.data == data

    sock = ClickHouseSock(PipeBuffer())
    data = rand([
        ["ab", "bc", "cd"],
        ["ab", "ed", "ab"]
    ], nrows)
    column = Column("test", "Array(String)", data)
    chwrite(sock, column)

    res = read_col(sock, VarUInt(nrows))
    @test res.data == data

    sock = ClickHouseSock(PipeBuffer())
    data = rand([
        ["ab", "bc", "cd"],
        ["ab", "ed", "ab"]
    ], nrows)
    column = Column("test", "Array(LowCardinality(String))", data)
    chwrite(sock, column)

    res = read_col(sock, VarUInt(nrows))
    @test res.data == data
    sock = ClickHouseSock(PipeBuffer())
    data = rand([
        [["ab"], [missing, "cd"]],
        [["ab", "ac"], [missing, "ab"]]
    ], nrows)
    sock = ClickHouseSock(PipeBuffer())
    column = Column("test", "Array(Array(Nullable(String)))", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test string(data) == string(res.data)
    @test recursive_miss_cmp(data, res.data)

    sock = ClickHouseSock(PipeBuffer())
    data = rand([
        [["ab"], [missing, "cd"]],
        [["ab", "ac"], [missing, "ab"]]
    ], nrows)
    sock = ClickHouseSock(PipeBuffer())
    column = Column("test", "Array(Array(LowCardinality(Nullable(String))))", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test recursive_miss_cmp(data, res.data)
end

@testset "Nothing column" begin
    sock = ClickHouseSock(PipeBuffer())
    data = [missing, missing, missing, missing]
    column = Column("test", "Nullable(Nothing)", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(4))
    @test all(ismissing.(data))

    sock = ClickHouseSock(PipeBuffer())
    data = [[], [], [], []]
    column = Column("test", "Array(Nothing)", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(4))
    @test all(data .== Ref(Missing[]))
end

@testset "Int columns" begin

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = rand(Int64, nrows)
    column = Column("test", "SimpleAggregateFunction(sum, Int64)", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test res == column

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = rand(Int32, nrows)
    column = Column("test", "Int64", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test res == column

end

@testset "Decimal columns" begin

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = Dec32.(rand(1000:9999, nrows), -3)
    column = Column("test", "Decimal32(3)", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test res == column

    sock = ClickHouseSock(PipeBuffer())
    column = Column("test", "Decimal(4,3)", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test res == column

    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = Dec64.(rand(1000000000:9999999999, nrows), -4)
    column = Column("test", "Decimal64(4)", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test res == column

    sock = ClickHouseSock(PipeBuffer())
    column = Column("test", "Decimal(10,4)", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test res == column
    sock = ClickHouseSock(PipeBuffer())
    nrows = 100
    data = Dec128.(rand(Int128(10)^20:(Int128(10)^21 - 1), nrows), -14)
    column = Column("test", "Decimal128(14)", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test res == column

    sock = ClickHouseSock(PipeBuffer())
    column = Column("test", "Decimal(20,14)", data)
    chwrite(sock, column)
    res = read_col(sock, VarUInt(nrows))
    @test res == column

end


text:
# This version shows how verifiability for proof of shuffle is added

using ShuffleProofs
import ShuffleProofs: PrimeGenerator, Enc, Dec, shuffle, ElGamal, ProtocolSpec, verify, Shuffle

sk = 123 # Only authorithy who do decryption would know

g = PrimeGenerator(3, 23)
options = [g, g^2, g^3]
pk = g^sk

verifier = ProtocolSpec(; g) # Some default parameters. In practice needs to be chosen carefully.

bbord = (;) # A representation of a buletinboard

######## Step 1. ########

bbord = let
    enc = Enc(pk, g)
    ciphertexts_in = ElGamal([enc(options[rand(1:3)], rand(1:10)) for i in 1:10])
    (; bbord..., ciphertexts_in)
end

######## Step 2. ########

bbord = let
    enc = Enc(pk, g)
    simulator = shuffle(bbord.ciphertexts_in, enc, verifier) # now returns simulator as verifier is added
    ciphertexts_out = simulator.proposition.𝐞′
    (; bbord..., ciphertexts_out, simulator.proof)
end

####### Step 3. ########

bbord = let
    dec = Dec(sk)
    votes = dec(bbord.ciphertexts_out)
    (; bbord..., votes)
end

######### As last step we can tally the votes ############

# Before counting we perform an audit that reencryption have happened honestly
proposition = Shuffle(g, pk, bbord.ciphertexts_in, bbord.ciphertexts_out)
@assert verify(proposition, bbord.proof, verifier)

for (i, o) in enumerate(options)
    s = 0
    for v in bbord.votes
        if v == o
            s += 1
        end
    end

    println("Choice $i received $s votes")
end


text:
module ExaTronCUDAKernels
    using ..ExaTronKernels
    using ..CUDA
    include("daxpy.jl")
    include("dcopy.jl")
    include("ddot.jl")
    include("dmid.jl")
    include("dnrm2.jl")
    include("dgpnorm.jl")
    include("dscal.jl")
    include("dssyax.jl")
    include("dnsol.jl")
    include("dtsol.jl")
    include("dtrqsol.jl")
    include("dbreakpt.jl")
    include("dgpstep.jl")
    include("dicf.jl")
    include("dicfs.jl")
    include("dprsrch.jl")
    include("dcauchy.jl")
    include("dtrpcg.jl")
    include("dspcg.jl")
    include("dtron.jl")
    include("TronMatrix.jl")
    ExaTronKernels.supports(::ExaTronKernels.CUDAKernels) = true
end


text:
import Base.@deprecate_binding

@deprecate_binding UfixedBase UFixed
@deprecate_binding UfixedConstructor UFixedConstructor
@deprecate_binding Ufixed UFixed
@deprecate_binding Ufixed8 UFixed8
@deprecate_binding Ufixed10 UFixed10
@deprecate_binding Ufixed12 UFixed12
@deprecate_binding Ufixed14 UFixed14
@deprecate_binding Ufixed16 UFixed16

@deprecate_binding Fixed32 Fixed16
@deprecate Fixed(x::Real) convert(Fixed{Int32, 16}, x)

@deprecate ufixed8(x)  UFixed8(x)
@deprecate ufixed10(x) UFixed10(x)
@deprecate ufixed12(x) UFixed12(x)
@deprecate ufixed14(x) UFixed14(x)
@deprecate ufixed16(x) UFixed16(x)

Compat.@dep_vectorize_1arg Real ufixed8
Compat.@dep_vectorize_1arg Real ufixed10
Compat.@dep_vectorize_1arg Real ufixed12
Compat.@dep_vectorize_1arg Real ufixed14
Compat.@dep_vectorize_1arg Real ufixed16


text:
"""
    struct IntegralOperator{T,K,S,V} <: AbstractMatrix{T}

Representation of an integral operator which takes a density ϕ defined on `Y`
and integrates it with `kernel` for all elements `x ∈ X`.
"""
struct IntegralOperator{T,K,S,V} <: AbstractMatrix{T}
    kernel::K
    X::S
    Y::V
end

IntegralOperator{T}(k,X,Y) where {T} = IntegralOperator{T,typeof(k),typeof(X),typeof(Y)}(k,X,Y)

function IntegralOperator(k,X,Y=X)
    T = return_type(k,X,Y)
    assert_concrete_type(T)
    IntegralOperator{T}(k,X,Y)
end

Base.size(iop::IntegralOperator)      = (length(iop.X), length(iop.Y))

kernel_type(iop::IntegralOperator) = kernel_type(iop.kernel)

kernel(iop::IntegralOperator) = iop.kernel

target_surface(iop::IntegralOperator) = iop.X
source_surface(iop::IntegralOperator) = iop.Y

Base.getindex(iop::IntegralOperator,i::Integer,j::Integer)  = getindex(kernel_type(iop),iop,i,j)

function Base.getindex(::SingleLayer,iop::IntegralOperator,i::Integer,j::Integer)
    x,y,w = qnodes(iop.X)[i], qnodes(iop.Y)[j], qweights(iop.Y)[j]
    return iop.kernel(x,y)*w
end
function Base.getindex(::DoubleLayer,iop::IntegralOperator,i::Integer,j::Integer)
    x,y,ny,w = qnodes(iop.X)[i], qnodes(iop.Y)[j], qnormals(iop.Y)[j], qweights(iop.Y)[j]
    return iop.kernel(x,y,ny)*w
end
function Base.getindex(::AdjointDoubleLayer,iop::IntegralOperator,i::Integer,j::Integer)
    x,y,nx,w = qnodes(iop.X)[i], qnodes(iop.Y)[j], qnormals(iop.X)[i], qweights(iop.Y)[j]
    return iop.kernel(x,y,nx)*w
end
function Base.getindex(::HyperSingular,iop::IntegralOperator,i::Integer,j::Integer)
    x,y,nx,ny,w = qnodes(iop.X)[i], qnodes(iop.Y)[j], qnormals(iop.X)[i], qnormals(iop.Y)[j], qweights(iop.Y)[j]
    return iop.kernel(x,y,nx,ny)*w
end

combined_field_coefficients(iop::IntegralOperator) = combined_field_coefficients(iop.kernel)

# convenience constructors
SingleLayerOperator(op::AbstractPDE,X,Y=X)        = IntegralOperator(SingleLayerKernel(op),X,Y)
DoubleLayerOperator(op::AbstractPDE,X,Y=X)        = IntegralOperator(DoubleLayerKernel(op),X,Y)
AdjointDoubleLayerOperator(op::AbstractPDE,X,Y=X) = IntegralOperator(AdjointDoubleLayerKernel(op),X,Y)
HyperSingularOperator(op::AbstractPDE,X,Y=X)      = IntegralOperator(HyperSingularKernel(op),X,Y)

ambient_dimension(iop::IntegralOperator) = ambient_dimension(iop.kernel)

function isinside(x::SVector,mesh::NystromMesh)
    N   = ambient_dimension(mesh)
    pde = Laplace(dim=N)
    K   = DoubleLayerKernel(pde)
    y   = qnodes(mesh)
    ν   = qnormals(mesh)
    w   = qweights(mesh)
    u   = sum(zip(y,ν,w)) do (y,ν,w)
         K(x,y,ν)*w
    end
    u + 0.5 < 0
end
isinside(x::Tuple,mesh::NystromMesh) = isinside(SVector(x),mesh)

qnodes(vec::Array{<:SVector}) = vec

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
## check that signal's and slots work

using PySide

## test it with a Qt example
w = Widget()
setWindowTitle(w, "Example")

lcd = Qt[:QLCDNumber](w.w)
sld = Slider(w)
btn = Button(w)
setText(btn, "click me")

vbox = VBoxLayout(w)
map(u -> addWidget(vbox, u), (lcd, sld, btn))

setLayout(w, vbox)

## connect valueChanged signal of sld with display slot of lcd
qconnect(sld, :valueChanged, lcd[:display])
qconnect(btn, :clicked) do
  println("clicked button. Value is:", get_value(sld))
end

raise(w)


text:
begin
    using CSV
    using Revise
    using BNPRegressionGGA2021
    using DataFrames
    using DataStructures
    using Distributions
    using Random
    using RCall
    using Statistics
    using StatsBase
    using LinearAlgebra
    const BNP = BNPRegressionGGA2021
    import Distributions: pdf
    import Base: rand
end;

# Helper functions ------------------------------------------------------------

function run_experiment_bnp(dy; N0, Nrep, Niter, id, filename)
    println(filename)
    Random.seed!(1) # seed
    samples = [generate_sample(dy; N0, Nrep) for _ in 1:Niter] # data
    gammas = [get_gamma_bnp(samples[id], id) for id in 1:Niter] # map gamma
    reduced_gammas = reduce_gammas(gammas; method = "bnp", id = id, N0 = N0)
    CSV.write(filename, reduced_gammas) # csv file
end

function run_experiment_freq(dy; N0, Nrep, Niter, id, filename)
    Random.seed!(1) # seed
    samples = [generate_sample(dy; N0, Nrep) for _ in 1:Niter] # data
    gammas = [get_gamma_freq(samples[id]) for id in 1:Niter] # map gamma
    reduced_gammas = reduce_gammas(gammas; method = "freq", id = id, N0 = N0)
    CSV.write(filename, reduced_gammas) # csv file
end

function generate_sample(dy; N0 = 50, Nrep = 10)
    # Simulate the data
    X0d = rand([0, 1], N0, 4)
    x0 = repeat(LinRange(-2, 2, N0 ÷ Nrep), Nrep)
    y0 = @. rand(dy(X0d[:, 4]))
    event0 = y0 .< 5

    # Generate the grid
    N1 = 0
    x1 = LinRange(-1, 1, N1)
    y1 = LinRange(-3, 3, 50) |> collect

    # Expand the grid
    grid = rcopy(R"""expand.grid(y1 = $y1, x1 = $x1)""")
    y1 = grid[!, :y1]
    x1 = grid[!, :x1]

    # Expand x0-x1 using splines
    R"""
    x0min <- min($x0)
    x0max <- max($x0)
    X0 <- splines::bs($x0, df = 6, Boundary.knots = c(x0min, x0max))
    """
    @rget X0

    # Standardise responses and predictors
    mX0, sX0 = mean_and_std(X0, 1)
    for col in 1:size(X0, 2)
        if sX0[col] > 0
            X0[:, col] = (X0[:, col] .- mX0[col]) ./ sX0[col]
        end
    end
    my0, sy0 = mean_and_std(y0, 1)
    y0 = (y0 .- my0) ./ sy0

    # Add a constant term to the design matrices
    X0d = ((X0d .- mean(X0d, dims = 1)) ./ std(X0d, dims = 1))
    X0 = [ones(size(X0, 1)) X0 X0d]
    X1 = zeros(0, size(X0, 2))

    # Generate the true (standardized) density over the grid
    f1 = pdf.(dy.(x1), my0 .+ y1 .* sy0) .* sy0

    # Set the mapping
    mapping = [[1], collect(2:7), [8], [9], [10], [11]]

    # Return the preprocessed results
    return y0, event0, x0, y1, x1, X0, X1, f1, mapping
end;

# Get the MAP estimator of gamma using BNP
function get_gamma_bnp(sample, id)
    println(id)
    y0, event0, x0, y1, x1, X0, X1, f1, mapping = sample
    smpl = BNP.DGSBPNormalDependent(; y0, event0, X0, y1, X1, mapping)
    _, _, chaing = BNP.sample!(smpl; mcmcsize = 10000)
    gamma = begin
        chaing |>
        DataStructures.counter |>
        collect |>
        x -> sort(x, by = x -> x[2], rev = true) |>
        x -> getindex(x, 1) |>
        x -> getindex(x, 1)
    end
    return gamma
end

# Estimate gamma using a frequentist alternative
function get_gamma_freq(sample)
    y0, event0, x0, y1, x1, X0, X1, f1, mapping = sample
    R"""
    library(dplyr)
    data = data.frame(
        status = $event0,
        y = $y0,
        x2 = $x0,
        x3 = $X0[, 8],
        x4 = $X0[, 9],
        x5 = $X0[, 10],
        x6 = $X0[, 11]
    )
    fitted_model <- survival::coxph(survival::Surv(y, status) ~ ., data = data)
    fitted_bestmodel <- MASS::stepAIC(fitted_model)
    best_model_trues <-
        names(fitted_bestmodel$coefficients) %>%
        gsub("x", "", .) %>%
        as.numeric()
    best_gamma <- rep(FALSE, 6)
    best_gamma[best_model_trues] <- TRUE
    best_gamma[1] <- TRUE
    """
    @rget best_gamma
    return best_gamma
end

# Reduce the MAP estimators after their estimation
function reduce_gammas(gammas; method, id, N0)
    names = "g" .* string.(1:6)
    mat = reduce(vcat, gammas')
    df = DataFrame(mat, names)
    df[!, :method] .= method
    df[!, :id] .= id
    df[!, :N0] .= N0
    return df
end

# Experiment 1, N0 = 50 (bnp)
begin
    filename = "data/final/survival-gamma-loggamma-bnp-50.csv"
    dy(x) = MixtureModel(Normal, [(3 + x, 0.8 + 0.2x), (3 - x, 0.8)], [.4, .6])
    # dy(x) = Normal(2 + x, 1)
    run_experiment_bnp(dy; N0 = 50, Nrep = 5, Niter = 100, id = 1, filename)
end

# Experiment 1, N0 = 50 (freq)
begin
    filename = "data/final/survival-gamma-loggamma-freq-50.csv"
    dy(x) = MixtureModel(Normal, [(3 + x, 0.8 + 0.2x), (3 - x, 0.8)], [.4, .6])
    # dy(x) = Normal(2 + x, 1)
    run_experiment_freq(dy; N0 = 50, Nrep = 5, Niter = 100, id = 1, filename)
end

# Experiment 1, N0 = 100 (bnp)
begin
    filename = "data/final/survival-gamma-loggamma-bnp-100.csv"
    dy(x) = MixtureModel(Normal, [(3 + x, 0.8 + 0.2x), (3 - x, 0.8)], [.4, .6])
    # dy(x) = Normal(2 + x, 1)
    run_experiment_bnp(dy; N0 = 100, Nrep = 5, Niter = 100, id = 1, filename)
end

# Experiment 1, N0 = 100 (freq)
begin
    filename = "data/final/survival-gamma-loggamma-freq-100.csv"
    dy(x) = MixtureModel(Normal, [(3 + x, 0.8 + 0.2x), (3 - x, 0.8)], [.4, .6])
    # dy(x) = Normal(2 + x, 1)
    run_experiment_freq(dy; N0 = 100, Nrep = 5, Niter = 100, id = 1, filename)
end

# Experiment 1, N0 = 200 (bnp)
begin
    filename = "data/final/survival-gamma-loggamma-bnp-200.csv"
    dy(x) = MixtureModel(Normal, [(3 + x, 0.8 + 0.2x), (3 - x, 0.8)], [.4, .6])
    # dy(x) = Normal(2 + x, 1)
    run_experiment_bnp(dy; N0 = 200, Nrep = 5, Niter = 100, id = 1, filename)
end

# Experiment 1, N0 = 200 (freq)
begin
    filename = "data/final/survival-gamma-loggamma-freq-200.csv"
    dy(x) = MixtureModel(Normal, [(3 + x, 0.8 + 0.2x), (3 - x, 0.8)], [.4, .6])
    # dy(x) = Normal(2 + x, 1)
    run_experiment_freq(dy; N0 = 1000, Nrep = 5, Niter = 100, id = 1, filename)
end

# Experiment 1, N0 = 500 (bnp)
begin
    filename = "data/final/survival-gamma-loggamma-bnp-200.csv"
    dy(x) = MixtureModel(Normal, [(3 + x, 0.8 + 0.2x), (3 - x, 0.8)], [.4, .6])
    # dy(x) = Normal(2 + x, 1)
    run_experiment_bnp(dy; N0 = 500, Nrep = 5, Niter = 100, id = 1, filename)
end

# Experiment 1, N0 = 500 (freq)
begin
    filename = "data/final/survival-gamma-loggamma-freq-200.csv"
    dy(x) = MixtureModel(Normal, [(3 + x, 0.8 + 0.2x), (3 - x, 0.8)], [.4, .6])
    # dy(x) = Normal(2 + x, 1)
    run_experiment_freq(dy; N0 = 5000, Nrep = 5, Niter = 100, id = 1, filename)
end

# Experiment 1, N0 = 1000 (bnp)
begin
    filename = "data/final/survival-gamma-loggamma-bnp-1000.csv"
    dy(x) = MixtureModel(Normal, [(3 + x, 0.8 + 0.2x), (3 - x, 0.8)], [.4, .6])
    # dy(x) = Normal(2 + x, 1)
    run_experiment_bnp(dy; N0 = 1000, Nrep = 5, Niter = 100, id = 1, filename)
end

# Experiment 1, N0 = 1000 (freq)
begin
    filename = "data/final/survival-gamma-loggamma-freq-1000.csv"
    dy(x) = MixtureModel(Normal, [(3 + x, 0.8 + 0.2x), (3 - x, 0.8)], [.4, .6])
    # dy(x) = Normal(2 + x, 1)
    run_experiment_freq(dy; N0 = 1000, Nrep = 5, Niter = 100, id = 1, filename)
end

# Summary of the results
R"""
data <-
    list.files(
        path = "data/final",
        pattern = "survival-gamma*",
        full.names = TRUE
    ) |>
    purrr::map(readr::read_csv, show_col_types = FALSE) |>
    dplyr::bind_rows() |>
    dplyr::group_by(method, id, N0, g1, g2, g3, g4, g5, g6) |>
    dplyr::count(name = "frequency") |>
    dplyr::ungroup() |>
    dplyr::mutate(gamma = paste0(g1, g2, g3, g4, g5, g6)) |>
    dplyr::select(gamma, frequency, id, method, N = N0)
p <-
    data |>
    dplyr::mutate(
        method = method |>
            dplyr::recode(bnp = "BNP", freq = "BSS (Best Subset Selection)"),
        distribution = id |>
        dplyr::recode(`1` = "Normal", `2` = "Skew Normal"),
        N = factor(N, ordered = TRUE)
    ) |>
    ggplot2::ggplot(ggplot2::aes(y = gamma, x = N, fill = frequency)) +
    ggplot2::geom_tile(color = "grey90") +
    ggplot2::geom_text(ggplot2::aes(label = frequency), color = "red") +
    ggplot2::facet_grid(
        distribution ~ method,
        labeller = ggplot2::labeller(
            distribution = ggplot2::label_both,
            method = ggplot2::label_both
        )
    ) +
    ggplot2::theme_classic() +
    ggplot2::theme(legend.position = "top") +
    ggplot2::scale_fill_gradient(low = "white", high = "black") +
    ggplot2::labs(
        x = "N (sample size)",
        y = "gamma (hypothesis vector), as a string of 0s and 1s"
    )
ggplot2::ggsave("figures/final/heatmap-survival.png")
"""

# TODO:
# Save the number of clusters in each iteration
# Try a model with 3 variables and examine the results
# Plot the density marginal and conditional


text:
#验证所推导的求傅里叶系数公式正确性
#比较variance大小
using LinearAlgebra
using StatsBase
using Plots
using Notifier

function pt(arg...)
    flush(stdout)
    println(arg...)
end
pt()

M=4
m=9
p=normalize!(rand(M),1)#态对应的概率幅list
function A(x)
    sum(p .* cos.(Array(1:M) .* 2x))/2+0.5
end
xk=Array(1:m-1) .* pi ./ m
function est2(i)
    4/m*(1+sum(A.(xk) .* cos.(2i .* xk)))
end
function on(a,b)
    a%b==0 ? b-1 : -1
end
function est3(i)
    #跟est2本质一样，但舍弃est2中的浮点运算
    oo=0
    for j=1:M
        oo+=p[j]*(on(i+j,m)+on(i-j,m))
    end
    4/m*(1+1/2*on(i,m)+1/4*oo)
end
for i=1:M
    pt()
    est2(i)|>pt
    p[i]|>pt
    est3(i)|>pt
end
function D1(N)
    o=0
    for i=1:M
        o+=16/m^2/N*sum(A.(xk) .* (1 .- A.(xk)) .* cos.(2i .* xk).^2)
    end
    o
end
function D2(N)
    o=0
    for i=1:M
        o+=p[i]*(1-p[i])/N/(m-1)
    end
    o
end
pt("variance")
D1(3)|>pt
D2(3)|>pt
# plot(D1.(1:100))
# plot!(D2.(1:100))
# for i=1:m,j=1:m
#     pt(i," ",j)
#     sum(cos.(2i .* Array(1:m-1) .* pi ./ m) .* cos.(2j .* Array(1:m-1) .* pi ./ m))|>pt
# end

max_stars_repo_name: helloinrm/readout

max_issues_repo_name: helloinrm/readout

max_forks_repo_name: helloinrm/readout


text:
using Soss
using Test
using MeasureTheory

include("examples-list.jl")

@testset "Soss.jl" begin
    @testset "Unit tests" begin
        @testset "Linear model" begin
            # include("linear-model.jl")
        end

        @testset "Transforms" begin
            include("transforms.jl")
        end
    end

    @testset "Examples" begin
        for example in EXAMPLES
            @testset "Run example: $(example[1])" begin
                # example_file = joinpath(EXAMPLESROOT, "example-$(example[2]).jl")
                # extra_example_tests = joinpath(TESTROOT, "extra-example-tests", "$(example[2]).jl")
                # @info("Running $(example_file)")
                # include(example_file)
                # if isfile(extra_example_tests)
                #     @info("Running $(extra_example_tests)")
                #     include(extra_example_tests)
                # end
            end
        end
    end

    @testset "Nested models" begin
        inner = @model a, b begin
            p ~ Beta(a, b)
            x ~ Normal(p, 1.0) |> iid(3)
            return x
        end

        outer = @model sub begin
            a ~ Beta(0.5, 0.5)
            b ~ Beta(1, 0.5)
            m ~ sub(a = a, b = b)
        end

        x = rand(outer(sub=inner)).m
        post = outer(sub=inner) | (;m=  (; x))
        t = xform(post)
        @test logdensity(post, t(randn(3))) isa Float64
    end

    @testset "Doctests" begin
        include("doctests.jl")
    end
end


text:
module Dynamic

    using ProjAssistant
    @gen_sub_proj

    using Plots
    import MathProgBase.HighLevelInterface: linprog
    import Clp: ClpSolver
    using ProgressMeter
    using Base.Threads
    using ExtractMacro
    using Statistics
    using Random
    import SimTools
    import SimTools: grad_desc, gd_value

    include("BoxGrid.jl")
    include("LP.jl")
    include("MetNets.jl")
    include("ToyModelD2.jl")
    include("ToyModelD3.jl")
    include("PlotUtils.jl")
    include("Space.jl")
    include("Container.jl")
    include("utils.jl")
    include("SimD2.jl")
    include("SimD3.jl")
    include("subspace.jl")
    include("SimTools.jl")
    include("STST.jl")
    include("maxent.jl")

end


text:
using LinearAlgebra, LowRankApprox, Formatting, ImageFiltering, Images, Plots, ImageMagick, Colors, TestImages, ImageView, ImageTransformations

I=load("rgb.jpg");
X=imresize(I, ratio=1/5);
(p,q)=size(X);
A=channelview(X);
R=Array{Float64}(A[1,:,:]);
G=Array{Float64}(A[2,:,:]);
B=Array{Float64}(A[3,:,:]);

RU, RS, RV=psvd(R);
GU, GS, GV=psvd(G);
BU, BS, BV=psvd(B);

R010=sum(RS[n]*RU[1:p,n]*(RV[1:q,n])' for n=1:10);
R020=sum(RS[n]*RU[1:p,n]*(RV[1:q,n])' for n=1:20);
R030=sum(RS[n]*RU[1:p,n]*(RV[1:q,n])' for n=1:30);
R040=sum(RS[n]*RU[1:p,n]*(RV[1:q,n])' for n=1:40);
R050=sum(RS[n]*RU[1:p,n]*(RV[1:q,n])' for n=1:50);
G010=sum(GS[n]*GU[1:p,n]*(GV[1:q,n])' for n=1:10);
G020=sum(GS[n]*GU[1:p,n]*(GV[1:q,n])' for n=1:20);
G030=sum(GS[n]*GU[1:p,n]*(GV[1:q,n])' for n=1:30);
G040=sum(GS[n]*GU[1:p,n]*(GV[1:q,n])' for n=1:40);
G050=sum(GS[n]*GU[1:p,n]*(GV[1:q,n])' for n=1:50);
B010=sum(BS[n]*BU[1:p,n]*(BV[1:q,n])' for n=1:10);
B020=sum(BS[n]*BU[1:p,n]*(BV[1:q,n])' for n=1:20);
B030=sum(BS[n]*BU[1:p,n]*(BV[1:q,n])' for n=1:30);
B040=sum(BS[n]*BU[1:p,n]*(BV[1:q,n])' for n=1:40);
B050=sum(BS[n]*BU[1:p,n]*(BV[1:q,n])' for n=1:50);

A010=zeros(3,p,q);
A020=zeros(3,p,q);
A030=zeros(3,p,q);
A040=zeros(3,p,q);
A050=zeros(3,p,q);

A010[1,:,:]=R010;
A010[2,:,:]=G010;
A010[3,:,:]=B010;
A020[1,:,:]=R020;
A020[2,:,:]=G020;
A020[3,:,:]=B020;
A030[1,:,:]=R030;
A030[2,:,:]=G030;
A030[3,:,:]=B030;
A040[1,:,:]=R040;
A040[2,:,:]=G040;
A040[3,:,:]=B040;
A050[1,:,:]=R050;
A050[2,:,:]=G050;
A050[3,:,:]=B050;

X010=colorview(RGB, A010);
X020=colorview(RGB, A020);
X030=colorview(RGB, A030);
X040=colorview(RGB, A040);
X050=colorview(RGB, A050);

plot(X,
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="A=A(365) 365 by 548 full rank matrix")
savefig("rgb_svd_365.png")

plot(X010,
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="A(10) rank 10 approximation")
savefig("rgb_svd_010.png")

plot(X020,
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="A(20) rank 20 approximation")
savefig("rgb_svd_020.png")

plot(X030,
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="A(30) rank 30 approximation")
savefig("rgb_svd_030.png")

plot(X040,
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="A(40) rank 40 approximation")
savefig("rgb_svd_040.png")

plot(X050,
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="A(50) rank 50 approximation")
savefig("rgb_svd_050.png")


text:
__precompile__()

module FEBOLPlots

using FEBOL

using PyPlot:
    imshow,
    xlabel,
    ylabel,
    contour,
    figure,
    pause,
    hold,
    axis,
    title,
    scatter,
    gcf,
    savefig,
    matplotlib,
    rc,
    tick_params,
    clf,
    cla

using PyPlot
import PyPlot.plot

using Reel      # for gifs

#rc("font", family="serif")
#rc("text", usetex=true)
rc("font", family="Times New Roman", size=16)
println("rc_context = ", rc_context())

export
    visualize,
    gif,
    plot_vehicle,
    plot,
    hold,
    title,
    pause,
    figure,
    cla,
    gcf,
    savefig,
    Frames

const LocTuple = NTuple{2, Float64}

include("visualize.jl")
include("plotting.jl")
include("gif.jl")
include("special.jl")

end # module


text:
using Test

@testset "run-FlightSims-simulations" begin
    include("environments/environments.jl")
end