text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb04vd(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    MMAX = 20
    NMAX = 20
    LDA  = MMAX
    LDE = MMAX
    LDQ = MMAX
    LDZ = NMAX
    LINUK = max( NMAX,MMAX+1 )
    LIWORK = NMAX
    LDWORK = max( NMAX,MMAX )
    ZERO = 0.0e0
    ONE = 1.0e0
    A = Array{Float64,2}(undef, LDA,NMAX)
    E = Array{Float64,2}(undef, LDE,NMAX)
    Q = Array{Float64,2}(undef, LDQ,MMAX)
    Z = Array{Float64,2}(undef, LDZ,NMAX)
    ISTAIR = Array{BlasInt,1}(undef, MMAX)
# WARNING: desperate attempt to initialize NBLCKS
    NBLCKS = 0
# WARNING: desperate attempt to initialize NBLCKI
    NBLCKI = 0
    IMUK = Array{BlasInt,1}(undef, LINUK)
    INUK = Array{BlasInt,1}(undef, LINUK)
    IMUK0 = Array{BlasInt,1}(undef, NMAX)
    MNEI = Array{BlasInt,1}(undef, 3)
    DWORK = Array{Float64,1}(undef, LDWORK)
    IWORK = Array{BlasInt,1}(undef, LIWORK)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    M = parse(BlasInt, vs[1])
    N = parse(BlasInt, vs[2])
    TOL = parse(Float64, replace(vs[3],'D'=>'E'))
    MODE =  vs[4][1]
    if ( M<0 || M>MMAX )
        @error "Illegal M=$M"
    end
    if ( N<0 || N>NMAX )
        @error "Illegal N=$N"
    end

    vs = String[]
    _isz,_jsz = (M,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       A[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (M,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       E[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
    close(f)
    JOBQ = 'I'
    JOBZ = 'I'
# interp call 1

    RANKE, INFO = SLICOT.mb04ud!(JOBQ, JOBZ, M, N, A, E, Q, Z, ISTAIR, TOL)
    println(io, "RANKE = $RANKE")
    @test INFO == 0
    INFO == 0 || return

    JOBQ = 'U'
    JOBZ = 'U'

# interp call 2

    NBLCKS, NBLCKI, INFO = SLICOT.mb04vd!(MODE, JOBQ, JOBZ, M, N, RANKE, A, E, Q, Z, ISTAIR,
                                          IMUK, INUK, IMUK0, MNEI, TOL)
    @test INFO == 0
    INFO == 0 || return

# interp output 1
    println(io, "Q:")
    _nc = M
    _nr = M
    show(io, "text/plain", Q[1:_nr,1:_nc])
    println(io)

# interp output 2
    println(io, "E:")
    _nc = N
    _nr = M
    show(io, "text/plain", E[1:_nr,1:_nc])
    println(io)

# interp output 3
    println(io, "A:")
    _nc = N
    _nr = M
    show(io, "text/plain", A[1:_nr,1:_nc])
    println(io)

# interp output 4
    println(io, "Z:")
    _nc = N
    _nr = N
    show(io, "text/plain", Z[1:_nr,1:_nc])
    println(io)

if ( ! LSAME( MODE, 'S' ) )
# interp output 5
    println(io, "IMUK:")
    _nr = NBLCKS
    show(io, "text/plain", IMUK[1:_nr])
    println(io)

# interp output 6
    println(io, "INUK:")
    _nr = NBLCKS
    show(io, "text/plain", INUK[1:_nr])
    println(io)

else
# interp output 7
    println(io, "IMUK:")
    _nr = NBLCKS
    show(io, "text/plain", IMUK[1:_nr])
    println(io)

# interp output 8
    println(io, "INUK:")
    _nr = NBLCKS
    show(io, "text/plain", INUK[1:_nr])
    println(io)

# interp output 9
    println(io, "IMUK0:")
    _nr = NBLCKI
    show(io, "text/plain", IMUK0[1:_nr])
    println(io)

# interp output 10
    println(io, "MNEI:")
    _nr = 3
    show(io, "text/plain", MNEI[1:_nr])
    println(io)
end # if
end # run_mb04vd()


text:
using MPI
using ScaLapack
using ScaLapack: BLACS, ScaLapackLite

# for check correction
using LinearAlgebra

# for DEBUG
using Printf
const DEBUG = true
function prtf(x)
    if x < 0; s = @sprintf("%.7f", x);
    else s = @sprintf(" %.7f", x); end;
end

# problem size
const nrows = 20
const ncols = 20
const nrows_block = 6
const ncols_block = 6
const nprocrows = 2
const nproccols = 2

# test type
const testtype = Float64

# finalizer
function mpi_finalizer(comm)
    MPI.Finalize()
end

function mutiple_test(root, comm)

    rank = MPI.Comm_rank(comm)
    nproc = MPI.Comm_size(comm)
    MPI.Barrier(comm);

    # generate matrix
    if rank == 0
        A = Matrix{testtype}(undef, nrows, ncols)
        B = Matrix{testtype}(undef, nrows, ncols)
        for ia::Integer = 1 : nrows
            for ja::Integer = 1 : ncols
                A[ia, ja] = convert(Float64, ia+ja)
                B[ia, ja] = convert(Float64, ia*ja)
            end
        end
    else
        A = Matrix{testtype}(undef, 0, 0)
        B = Matrix{testtype}(undef, 0, 0)
    end
    MPI.Barrier(comm)

    if DEBUG
        if rank == root
            print("\ninput matrix:\n")
        end
        MPI.Barrier(comm)
        print("\nrank[$rank];\nA = \n$A\nB = \n$B\n")
    end
    MPI.Barrier(comm)

    # prepare ScaLapackLite params
    params = ScaLapackLite.ScaLapackLiteParams(nrows_block, ncols_block, nprocrows, nproccols, root)
    slm_A = ScaLapackLite.ScaLapackLiteMatrix(params, A)
    slm_B = ScaLapackLite.ScaLapackLiteMatrix(params, B)

    # perform C = A*B
    slm_C = slm_A * slm_B'
    if DEBUG
        if rank == root
            print("\noutput matrix:\n")
        end
        MPI.Barrier(comm)
        print("\nrank[$rank];\nC = \n$(slm_C.X)\n")
    end

    MPI.Barrier(comm)
    # clean up
    BLACS.exit()

end

function hessenberg_test(root, comm)
    rank = MPI.Comm_rank(comm)
    nproc = MPI.Comm_size(comm)
    MPI.Barrier(comm);

    A = Matrix{testtype}(undef, nrows, ncols)
    for ia::Integer = 1 : nrows
        for ja::Integer = 1 : ncols
            A[ia, ja] = convert(testtype, rand())
        end
    end
    params = ScaLapackLite.ScaLapackLiteParams(nrows_block, ncols_block, nprocrows, nproccols, root)
    slm_A = ScaLapackLite.ScaLapackLiteMatrix(params, A)

    slm_hA = ScaLapackLite.hessenberg(slm_A, true)
    print("[$rank] A: $(slm_hA.X)\n")
end

function schur_test(root, comm)
    rank = MPI.Comm_rank(comm)
    nproc = MPI.Comm_size(comm)
    MPI.Barrier(comm);

    A = Matrix{testtype}(undef, nrows, ncols)
    for ia::Integer = 1 : nrows
        for ja::Integer = 1 : ncols
            A[ia, ja] = convert(testtype, rand())
        end
    end
    params = ScaLapackLite.ScaLapackLiteParams(nrows_block, ncols_block, nprocrows, nproccols, root)
    slm_A = ScaLapackLite.ScaLapackLiteMatrix(params, A)

    if rank == root
        # reference: original matrix A
        A_ref = deepcopy(A)
        # reference: LAPACK routine
        (eig_ref,eigvec_ref)=LinearAlgebra.eigen(A)
    end

    # perform Schur decomposition and en-return T
    (eig, slm_Z, slm_T) = ScaLapackLite.schur(slm_A, true)

    # CHECK: recompose A from Schur vectors Z; A = Z * T * Z^H
    slm_A_ = slm_Z * slm_T * slm_Z'
    if rank == root
        frobenius_norm = norm(slm_A_.X-A_ref, 2)
    end

    if rank == root
        eidx=sortperm(real(eig))
        eidx_ref=sortperm(real(eig_ref))
        diff = [ abs(eig[eidx[j]]-eig_ref[eidx_ref[j]]) for j=1:nrows ]

        print("\n--- eigenvalues ---\n")
        wlines = "\nindex / eig(ScaLapackLite) / eig(LAPACK) / | eig(ScaLapackLite)-eig(LAPACK) |\n"
        for idx = 1:nrows
            wlines *= @sprintf("  %s / %s + i %s / %s + i %s / %s\n",
                          lpad(idx, 3, "0"),
                          prtf(real(eig[eidx[idx]])), prtf(imag(eig[eidx[idx]])),
                          prtf(real(eig_ref[eidx_ref[idx]])), prtf(imag(eig_ref[eidx_ref[idx]])),
                          prtf(diff[idx]))
        end
        wlines *= "\n"; print(wlines);

        print("\n--- Schur vectors / eigen vectors ---\n")

        wline = @sprintf("\nFrobenius norm: || Z * T * Z^H - A ||_{2} = %s\n", prtf(frobenius_norm))
        print(wline)

        for jdx = 1:ncols
            print("\n[$jdx-th vec]:\n")
            wlines = "\nindex / Schur vector (ScaLapackLite) / eig vector (LAPACK)\n"
            for idx = 1:nrows
                wlines *= @sprintf("  %s / %s + i %s / %s + i %s\n",
                            lpad(idx, 3, "0"),
                            prtf(real(slm_Z.X[idx, eidx[jdx]])),
                            prtf(imag(slm_Z.X[idx, eidx[jdx]])),
                            prtf(real(eigvec_ref[idx,eidx_ref[jdx]])),
                            prtf(imag(eigvec_ref[idx,eidx_ref[jdx]])))
            end
            wlines *= "\n"; print(wlines);
        end

    end

end

function eig_test(root, comm)
    rank = MPI.Comm_rank(comm)
    nproc = MPI.Comm_size(comm)
    MPI.Barrier(comm);

    A = Matrix{testtype}(undef, nrows, ncols)
    for ia::Integer = 1 : nrows
        for ja::Integer = 1 : ncols
            A[ia, ja] = convert(testtype, rand())
            # A[ia, ja] = convert(testtype, rand()+im*rand())
        end
    end
    A_copied = deepcopy(A)
    params = ScaLapackLite.ScaLapackLiteParams(nrows_block, ncols_block, nprocrows, nproccols, root)
    slm_A = ScaLapackLite.ScaLapackLiteMatrix(params, A)

    if rank == root
        # reference: LAPACK routine
        (eigvals_ref,eigvecs_ref)=LinearAlgebra.eigen(A)
    end

    # find eigenvalues and eigenvectors
    eigidx = Vector{Int64}([])
    (eigvals, eigvecs) = ScaLapackLite.eigen(slm_A, eigidx, false)

    testtype_e = testtype
    if testtype == Float32; testtype_e = ComplexF32;
    elseif testtype == Float64; testtype_e = ComplexF64; end

    if rank == root
        eidx = sortperm(real(eigvals))
        eidx_ref = sortperm(real(eigvals_ref))

        print("\n--- eigenvalues ---\n")
        wlines = "\nindex / eig(ScaLapackLite) / eig(LAPACK)\n"
        for idx = 1:nrows
            wlines *= @sprintf("  %s / %s + i %s / %s + i %s\n",
                          lpad(idx, 3, "0"),
                          prtf(real(eigvals[eidx[idx]])), prtf(imag(eigvals[eidx[idx]])),
                          prtf(real(eigvals_ref[eidx_ref[idx]])), prtf(imag(eigvals_ref[eidx_ref[idx]])))
        end
        wlines *= "\n"; print(wlines);

        print("\n--- eigv(ScaLapackLite) / eigv(LAPACK) ---\n")
        for jdx = 1:ncols

            # calc error
            x = Vector{testtype_e}(eigvecs[:,eidx[jdx]])
            x_ref = Vector{testtype_e}(eigvecs_ref[:,eidx_ref[jdx]])
            ej = eigvals[eidx[jdx]]
            ej_ref = eigvals_ref[eidx_ref[jdx]]
            error = norm(A_copied*x - ej*x, 2)
            error_ref = norm(A_copied*x_ref - ej_ref*x_ref, 2)
            wlines = @sprintf("\nlog10(|| Ax - λx ||)[ScaLapackLite] = %s\nlog10(|| Ax - λx ||)[LAPACK]        = %s\n",
                               prtf(log10(error)), prtf(log10(error_ref)))

            print("\n[$jdx-th vec]:\n")
            print(wlines)
            wlines = "\nindex / eigv(ScaLapackLite) / eigv(LAPACK)\n"
            for idx = 1:nrows
                wlines *= @sprintf("  %s / %s + i %s / %s + i %s\n",
                            lpad(idx, 3, "0"),
                            prtf(real(eigvecs[idx, eidx[jdx]])), prtf(imag(eigvecs[idx, eidx[jdx]])),
                            prtf(real(eigvecs_ref[idx, eidx_ref[jdx]])), prtf(imag(eigvecs_ref[idx, eidx_ref[jdx]])))
            end
            wlines *= "\n"; print(wlines);
        end

    end

end

function main()
    # MPI initilaize
    MPI.Init()
    comm = MPI.COMM_WORLD
    finalizer(mpi_finalizer, comm)
    # test
    # mutiple_test(0, comm)
    # hessenberg_test(0, comm)
    # schur_test(0, comm)
    eig_test(0, comm)
end

main()

ext: jl


text:
samples = [0.0, 1.0, 2.0, 3.0, 4.0]
h5file  = "temp.h5"

@testset "Parameter" begin

    @test_throws AssertionError Parameter(:μ, 1.0, 4.0, samples)
    @test_throws AssertionError Parameter(:μ, 0.0, 3.0, samples)
    @test_throws AssertionError Parameter(:μ, 1.0, 0.0, 5)
    @test_throws AssertionError Parameter(:μ, 1.0, 0.0)

    p1 = Parameter(:μ, 0.0, 4.0, samples)
    p2 = Parameter(:μ, 0.0, 4.0, samples)
    p3 = Parameter(:μ, 0.0, 4.0, 5)
    p4 = Parameter(:μ, 0.0, 4.0)
    p5 = Parameter(:μ, 0.0, 4.0, nothing)

    @test p1 == p2
    @test p1 == p3
    @test p1 != p4
    @test p4 == p5

    @test hash(p1) == hash(p2)
    @test hash(p1) == hash(p3)
    @test hash(p1) != hash(p4)
    @test hash(p4) == hash(p5)

    @test length(p1) == length(samples)
    @test length(p2) == length(samples)
    @test length(p3) == length(samples)
    @test length(p4) == 0
    @test length(p5) == 0

    @test size(p1) == size(samples)
    @test size(p2) == size(samples)
    @test size(p3) == size(samples)
    @test size(p4) == (0,)
    @test size(p5) == (0,)

    @test hassamples(p1) == true
    @test hassamples(p2) == true
    @test hassamples(p3) == true
    @test hassamples(p4) == false
    @test hassamples(p5) == false

    @test collect(p1) == samples
    @test collect(p2) == samples
    @test collect(p3) == samples
    @test collect(p4) === nothing
    @test collect(p5) === nothing

    @test minimum(p1) == 0.0
    @test minimum(p2) == 0.0
    @test minimum(p3) == 0.0
    @test minimum(p4) == 0.0
    @test minimum(p5) == 0.0

    @test maximum(p1) == 4.0
    @test maximum(p2) == 4.0
    @test maximum(p3) == 4.0
    @test maximum(p4) == 4.0
    @test maximum(p5) == 4.0

    p1 = Parameter(:μ, 0.0, 1.0, 3)
    p2 = Parameter(:ν, 1.0, 1.0, 1)
    p3 = Parameter(:σ, 0.0, 4.0, 2)

    @test NamedTuple(p1, p2, p3) == NamedTuple{(:μ, :ν, :σ)}((p1, p2, p3))

    h5save(p1, h5file; mode="w")
    @test isfile(h5file)

    p2 = h5load(Parameter, h5file, "μ")
    rm(h5file)
    @test p1 == p2

end


text:
using FourierFlows: CPU, TwoDGrid
using FFTW
using Random: seed!, rand

dev = CPU()           # Device.
N = 256               # Number of grid points.
L = 2.0               # Length of the square domain.
stepper = "RK4"       # Type of time stepper (e.g., RK4, AB3).
dt = 2e-3             # Size of timestep.
effort = FFTW.PATIENT # FFTW effort flag (e.g., FFTW.MEASURE, FFTW.PATIENT).
nlog = 500            # Display log and save Hovmöller slices every `nlog` timesteps.
β = 0.0               # Planetary background potential vorticity gradient.
U₀ = 1.0              # Amplitude of basic jet.
Lj = L/10             # Half-width of basic jet.
qmag = 0.004          # Amplitude of the random initial condition.
Rthresh = 0.01        # Threshold energy ratio for renormalization.
cycles = 3            # Total number of renormalization cycles.

cyclesm = cycles - 1
gr = TwoDGrid(N, L, N, L, x0=0.0, aliased_fraction=0.0)
U, Uyy = NondimLinearOneHalfLayerQG.get_bickley(gr, U₀, Lj)
seed!(1234) # Fix the seed of the random number generator for reproducibility.
q₀ = qmag * gr.nx * gr.ny * rand(gr.nx, gr.ny)

max_stars_repo_name: apaloczy/RoughOceanQGEnergetics

max_issues_repo_name: apaloczy/RoughOceanQGEnergetics

max_forks_repo_name: apaloczy/RoughOceanQGEnergetics


text:
"""
    complement(g)

Return the [graph complement](https://en.wikipedia.org/wiki/Complement_graph)
of a graph

### Implementation Notes
Preserves the eltype of the input graph.
"""
function complement(g::Graph)
    gnv = nv(g)
    h = Graph(gnv)
    for i=1:gnv
        for j=i+1:gnv
            if !has_edge(g, i, j)
                add_edge!(h,i,j)
            end
        end
    end
    return h
end

function complement(g::DiGraph)
    gnv = nv(g)
    h = DiGraph(gnv)
    for i in vertices(g), j in vertices(g)
      if i != j && !has_edge(g,i,j)
        add_edge!(h,i,j)
      end
    end
    return h
end

"""
    reverse(g)

Return a directed graph where all edges are reversed from the
original directed graph.

### Implementation Notes
Preserves the eltype of the input graph.
"""
function reverse end
@traitfn function reverse(g::::IsDirected)
    gnv = nv(g)
    gne = ne(g)
    h = DiGraph(gnv)
    h.fadjlist = deepcopy(g.badjlist)
    h.badjlist = deepcopy(g.fadjlist)
    h.ne = gne
    return h
end

"""
    reverse!(g)

In-place reverse of a directed graph (modifies the original graph).
"""
function reverse! end
@traitfn function reverse!(g::::IsDirected)
    g.fadjlist, g.badjlist = g.badjlist, g.fadjlist
    return g
end

doc"""
    blkdiag(g, h)

Return a graph with ``|V(g)| + |V(h)|`` vertices and ``|E(g)| + |E(h)|``
edges where the vertices an edges from graph `h` are appended to graph `g`.

### Implementation Notes
Preserves the eltype of the input graph. Will error if the
number of vertices in the generated graph exceeds the eltype.
"""
function blkdiag(g::T, h::T) where T<:AbstractGraph
    gnv = nv(g)
    r = T(gnv + nv(h))
    for e in edges(g)
        add_edge!(r, e)
    end
    for e in edges(h)
        add_edge!(r, gnv+src(e), gnv+dst(e))
    end
    return r
end

"""
    intersect(g, h)

Return a graph with edges that are only in both graph `g` and graph `h`.

### Implementation Notes
This function may produce a graph with 0-degree vertices.
Preserves the eltype of the input graph.
"""
function intersect(g::T, h::T) where T<:AbstractGraph
    gnv = nv(g)
    hnv = nv(h)

    r = T(min(gnv, hnv))
    for e in intersect(edges(g),edges(h))
        add_edge!(r,e)
    end
    return r
end

"""
    difference(g, h)

Return a graph with edges in graph `g` that are not in graph `h`.

### Implementation Notes
Note that this function may produce a graph with 0-degree vertices.
Preserves the eltype of the input graph.
"""
function difference(g::T, h::T) where T<:AbstractGraph
    gnv = nv(g)
    hnv = nv(h)

    r = T(gnv)
    for e in edges(g)
        !has_edge(h, e) && add_edge!(r,e)
    end
    return r
end

"""
    symmetric_difference(g, h)

Return a graph with edges from graph `g` that do not exist in graph `h`,
and vice versa.

### Implementation Notes
Note that this function may produce a graph with 0-degree vertices.
Preserves the eltype of the input graph. Will error if the
number of vertices in the generated graph exceeds the eltype.
"""
function symmetric_difference(g::T, h::T) where T<:AbstractGraph
    gnv = nv(g)
    hnv = nv(h)

    r = T(max(gnv, hnv))
    for e in edges(g)
        !has_edge(h, e) && add_edge!(r, e)
    end
    for e in edges(h)
        !has_edge(g, e) && add_edge!(r, e)
    end
    return r
end

"""
    union(g, h)

Return a graph that combines graphs `g` and `h` by taking the set union
of all vertices and edges.

### Implementation Notes
Preserves the eltype of the input graph. Will error if the
number of vertices in the generated graph exceeds the eltype.
"""
function union(g::T, h::T) where T<:AbstractGraph
    gnv = nv(g)
    hnv = nv(h)

    r = T(max(gnv, hnv))
    r.ne = ne(g)
    for i in vertices(g)
        r.fadjlist[i] = deepcopy(g.fadjlist[i])
        if is_directed(g)
            r.badjlist[i] = deepcopy(g.badjlist[i])
        end
    end
    for e in edges(h)
        add_edge!(r, e)
    end
    return r
end

"""
    join(g, h)

Return a graph that combines graphs `g` and `h` using `blkdiag` and then
adds all the edges between the vertices in `g` and those in `h`.

### Implementation Notes
Preserves the eltype of the input graph. Will error if the number of vertices
in the generated graph exceeds the eltype.
"""
function join(g::T, h::T) where T<:AbstractGraph
    r = blkdiag(g, h)
    for i in vertices(g)
        for j=nv(g)+1:nv(g)+nv(h)
            add_edge!(r, i, j)
        end
    end
    return r
end

"""
    crosspath(len::Integer, g::Graph)

Return a graph that duplicates `g` `len` times and connects each vertex
with its copies in a path.

### Implementation Notes
Preserves the eltype of the input graph. Will error if the number of vertices
in the generated graph exceeds the eltype.
"""
function crosspath end
@traitfn function crosspath(len::Integer, g::::(!IsDirected))
    T = eltype(g)
    p = PathGraph(len)
    h = Graph{T}(p)
    return cartesian_product(h, g)
end

# The following operators allow one to use a LightGraphs.Graph as a matrix in eigensolvers for spectral ranking and partitioning.
# """Provides multiplication of a graph `g` by a vector `v` such that spectral
# graph functions in [GraphMatrices.jl](https://github.com/jpfairbanks/GraphMatrices.jl) can utilize LightGraphs natively.
# """
function *(g::Graph, v::Vector{T}) where T<:Real
    length(v) == nv(g) || error("Vector size must equal number of vertices")
    y = zeros(T, nv(g))
    for e in edges(g)
        i = src(e)
        j = dst(e)
        y[i] += v[j]
        y[j] += v[i]
    end
    return y
end

function *(g::DiGraph, v::Vector{T}) where T<:Real
    length(v) == nv(g) || error("Vector size must equal number of vertices")
    y = zeros(T, nv(g))
    for e in edges(g)
        i = src(e)
        j = dst(e)
        y[i] += v[j]
    end
    return y
end

"""
    sum(g, i)

Return a vector of indegree (`i`=1) or outdegree (`i`=2) values for graph `g`.
"""
function sum(g::AbstractGraph, dim::Int)
    dim == 1 && return indegree(g, vertices(g))
    dim == 2 && return outdegree(g, vertices(g))
    error("Graphs are only two dimensional")
end

size(g::AbstractGraph) = (nv(g), nv(g))
"""
    size(g, i)

Return the number of vertices in `g` if `i`=1 or `i`=2, or `1` otherwise.
"""
size(g::Graph,dim::Int) = (dim == 1 || dim == 2)? nv(g) : 1

"""
    sum(g)

Return the number of edges in `g`
"""
sum(g::AbstractGraph) = ne(g)

"""
    sparse(g)

Return the default adjacency matrix of `g`.
"""
sparse(g::AbstractGraph) = adjacency_matrix(g)

#arrayfunctions = (:eltype, :length, :ndims, :size, :strides, :issymmetric)
# eltype(g::AbstractGraph) = Float64
length(g::AbstractGraph) = nv(g)*nv(g)
ndims(g::AbstractGraph) = 2
issymmetric(g::AbstractGraph) = !is_directed(g)

"""
    cartesian_product(g, h)

Return the (cartesian product)[https://en.wikipedia.org/wiki/Tensor_product_of_graphs]
of `g` and `h`.

### Implementation Notes
Preserves the eltype of the input graph. Will error if the number of vertices
in the generated graph exceeds the eltype.
"""
function cartesian_product(g::G, h::G) where G<:AbstractGraph
    z = G(nv(g)*nv(h))
    id(i, j) = (i-1)*nv(h) + j
    for e in edges(g)
        i1, i2 = Tuple(e)
        for j=1:nv(h)
            add_edge!(z, id(i1,j), id(i2,j))
        end
    end

    for e in edges(h)
        j1, j2 = Tuple(e)
        for i in vertices(g)
            add_edge!(z, id(i,j1), id(i,j2))
        end
    end
    return z
end

"""
    tensor_product(g, h)

Return the (tensor product)[https://en.wikipedia.org/wiki/Tensor_product_of_graphs]
of `g` and `h`.

### Implementation Notes
Preserves the eltype of the input graph. Will error if the number of vertices
in the generated graph exceeds the eltype.
"""
function tensor_product(g::G, h::G) where G<:AbstractGraph
    z = G(nv(g)*nv(h))
    id(i, j) = (i-1)*nv(h) + j
    for e1 in edges(g)
        i1, i2 = Tuple(e1)
        for e2 in edges(h)
            j1, j2 = Tuple(e2)
            add_edge!(z, id(i1, j1), id(i2, j2))
        end
    end
    return z
end

## subgraphs ###

"""
    induced_subgraph(g, vlist)
    induced_subgraph(g, elist)

Return the subgraph of `g` induced by the vertices in  `vlist` or edges in `elist`
along with a vector mapping the new vertices to the old ones
(the  vertex `i` in the subgraph corresponds to the vertex `vmap[i]` in `g`.)

The returned graph has `length(vlist)` vertices, with the new vertex `i`
corresponding to the vertex of the original graph in the `i`-th position
of `vlist`.

### Usage Examples
```doctestjl
julia> g = CompleteGraph(10)

julia> sg, vmap = induced_subgraph(g, 5:8)

julia> @assert g[5:8] == sg

julia> @assert nv(sg) == 4

julia> @assert ne(sg) == 6

julia> @assert vm[4] == 8

julia> sg, vmap = induced_subgraph(g, [2,8,3,4])

julia> @assert sg == g[[2,8,3,4]]

julia> elist = [Edge(1,2), Edge(3,4), Edge(4,8)]

julia> sg, vmap = induced_subgraph(g, elist)

julia> @assert sg == g[elist]
```
"""
function induced_subgraph(g::T, vlist::AbstractVector{U}) where T<:AbstractGraph where U<:Integer
    allunique(vlist) || error("Vertices in subgraph list must be unique")
    h = T(length(vlist))
    newvid = Dict{U, U}()
    vmap =Vector{U}(length(vlist))
    for (i,v) in enumerate(vlist)
        newvid[v] = U(i)
        vmap[i] = v
    end

    vset = Set(vlist)
    for s in vlist
        for d in out_neighbors(g, s)
            # println("s = $s, d = $d")
            if d in vset && has_edge(g, s, d)
                newe = Edge(newvid[s], newvid[d])
                add_edge!(h, newe)
            end
        end
    end
    return h, vmap
end

function induced_subgraph(g::T, elist::AbstractVector{U}) where T<:AbstractGraph where U<:AbstractEdge
    h = zero(g)
    et = eltype(h)
    newvid = Dict{et, et}()
    vmap = Vector{et}()

    for e in elist
        u, v = Tuple(e)
        for i in (u,v)
            if !haskey(newvid, i)
                add_vertex!(h)
                newvid[i] = nv(h)
                push!(vmap, i)
            end
        end
        add_edge!(h, newvid[u], newvid[v])
    end
    return h, vmap
end

"""
    g[iter]

Return the subgraph induced by `iter`.
Equivalent to [`induced_subgraph`](@ref)`(g, iter)[1]`.
"""
getindex(g::AbstractGraph, iter) = induced_subgraph(g, iter)[1]

"""
    egonet(g, v:, d)

Return the subgraph of `g` induced by the neighbors of `v` up to distance
`d`.
This is equivalent to [`induced_subgraph`](@ref)`(g, neighborhood(g, v, d, dir=dir))[1].`

### Optional Arguments
- `dir=:out`: if `g` is directed, this argument specifies the edge direction
with respect to `v` (i.e. `:in` or `:out`).
"""
egonet(g::AbstractGraph, v::Integer, d::Integer; dir=:out) = g[neighborhood(g, v, d, dir=dir)]

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
import Base: @deprecate

@deprecate by(d::AbstractDataFrame, cols, s::Vector{Symbol}) aggregate(d, cols, map(eval, s))
@deprecate by(d::AbstractDataFrame, cols, s::Symbol) aggregate(d, cols, eval(s))
@deprecate nullable!(colnames::Array{Symbol,1}, df::AbstractDataFrame) nullable!(df, colnames)
@deprecate nullable!(colnums::Array{Int,1}, df::AbstractDataFrame) nullable!(df, colnums)

import Base: keys, values, insert!
@deprecate keys(df::AbstractDataFrame) names(df)
@deprecate values(df::AbstractDataFrame) DataFrames.columns(df)
@deprecate insert!(df::DataFrame, df2::AbstractDataFrame) merge!(df, df2)

@deprecate read_rda(args...) FileIO.load(args...)

@deprecate pool categorical
@deprecate pool! categorical!


text:
using BinaryBuilder
using LibGit2
import SHA: sha256

function filehash(path)
    open(path, "r") do f
        bytes2hex(sha256(f))
    end
end

function installed_julia_dir(ver)
     jp = julia_path(ver)
     jp_contents = readdir(jp)
     # Allow the unpacked directory to either be insider another directory (as produced by
     # the buildbots) or directly inside the mapped directory (as produced by the BB script)
     if length(jp_contents) == 1
         jp = joinpath(jp, first(jp_contents))
     end
     jp
end

"""
    prepare_julia(the_ver)

Download the specified version of Julia using the information provided in `Versions.toml`.
"""
function prepare_julia(the_ver::VersionNumber)
    vers = read_versions()
    for (ver, data) in vers
        ver == string(the_ver) || continue
        dir = julia_path(ver)
        mkpath(dirname(dir))
        if haskey(data, "url")
            url = data["url"]

            file = get(data, "file", "julia-$ver.tar.gz")
            @assert !isabspath(file)
            file = downloads_dir(file)
            mkpath(dirname(file))

            Pkg.PlatformEngines.download_verify_unpack(url, data["sha"], dir;
                                                        tarball_path=file, force=true)
        else
            file = data["file"]
            !isabspath(file) && (file = downloads_dir(file))
            Pkg.PlatformEngines.verify(file, data["sha"])
            isdir(dir) || Pkg.PlatformEngines.unpack(file, dir)
        end
        return
    end
    error("Requested Julia version $the_ver not found")
end

"""
    version = download_julia(name::String)

Download Julia from an on-line source listed in Builds.toml as identified by `name`.
Returns the `version` (what other functions use to identify this build).
This version will be added to Versions.toml.
"""

function download_julia(name::String)
    builds = read_builds()
    @assert haskey(builds, name) "Julia build $name is not registered in Builds.toml"
    data = builds[name]

    # get the filename and extension from the url
    url = data["url"]
    filename = basename(url)
    if endswith(filename, ".tar.gz")
        ext = ".tar.gz"
        base = filename[1:end-7]
    else
        ext, base = splitext(filename)
    end

    # download
    temp_file = downloads_dir(base)
    mkpath(dirname(temp_file))
    ispath(temp_file) && rm(temp_file)
    Pkg.PlatformEngines.download(url, temp_file)

    # unpack
    temp_dir = julia_path(base)
    ispath(temp_dir) && rm(temp_dir; recursive=true)
    Pkg.PlatformEngines.unpack(temp_file, temp_dir)

    # figure out stuff from the downloaded binary
    version = VersionNumber(read(`$(installed_julia_dir(base))/bin/julia -e 'print(Base.VERSION_STRING)'`, String))
    if version.prerelease != ()
        commit_short = read(`$(installed_julia_dir(base))/bin/julia -e 'print(Base.GIT_VERSION_INFO.commit_short)'`, String)
        version = VersionNumber(string(version) * string("-", commit_short))
    end
    rm(temp_dir; recursive=true) # let `prepare_julia` unpack; keeps code simpler here

    versions = read_versions()
    if haskey(versions, string(version))
        @info "Julia $name (version $version) already available"
        rm(temp_file)
    else
        # always use the hash of the downloaded file to force a check during `prepare_julia`
        hash = filehash(temp_file)

        # move to its final location
        tarball = "julia-$version$ext"
        tarball_path = downloads_dir(tarball)
        if ispath(tarball_path)
            @warn "Destination file $tarball already exists, assuming it matches"
            rm(temp_file)
        else
            mv(temp_file, tarball_path)
        end

        # Update Versions.toml
        version_stanza = """
            ["$version"]
            file = "$tarball"
            sha = "$hash"
            """
        open(versions_file(); append=true) do f
            println(f, version_stanza)
        end
    end

    return version
end

"""
    version = build_julia(ref::String="master"; binarybuilder_args::Vector{String}=String["--verbose"])

Check-out and build Julia at git reference `ref` using BinaryBuilder.
Returns the `version` (what other functions use to identify this build).
This version will be added to Versions.toml.
"""
function build_julia(ref::String="master"; binarybuilder_args::Vector{String}=String["--verbose"])
    # get the Julia repo
    repo_path = downloads_dir("julia")
    if !isdir(repo_path)
        @info "Cloning Julia repository..."
        repo = LibGit2.clone("https://github.com/JuliaLang/julia", repo_path)
    else
        repo = LibGit2.GitRepo(repo_path)
        LibGit2.fetch(repo)
    end

    # lookup the version number and commit hash
    reference = LibGit2.GitCommit(repo, ref)
    tree = LibGit2.peel(LibGit2.GitTree, reference)
    version = VersionNumber(chomp(LibGit2.content(tree["VERSION"])))
    commit = string(LibGit2.GitHash(reference))
    if version.prerelease != ()
        commit_short = commit[1:10]
        version = VersionNumber(string(version) * string("-", commit_short))
    end

    versions = read_versions()
    if haskey(versions, string(version))
        @info "Julia $ref (version $version) already available"
        return version
    end

    # Collection of sources required to build julia
    sources = [
        "https://github.com/JuliaLang/julia.git" => commit,
    ]

    # Bash recipe for building across all platforms
    script = raw"""
    cd $WORKSPACE/srcdir
    mount -t devpts -o newinstance jrunpts /dev/pts
    mount -o bind /dev/pts/ptmx /dev/ptmx

    cd julia
    cat > Make.user <<EOF
    JULIA_CPU_TARGET=generic;sandybridge,-xsaveopt,clone_all;haswell,-rdrnd,base(1)
    EOF
    make -j${nproc}

    make install
    cp LICENSE.md ${prefix}
    contrib/fixup-libgfortran.sh ${prefix}/lib/julia
    contrib/fixup-libstdc++.sh ${prefix}/lib ${prefix}/lib/julia
    """

    # These are the platforms we will build for by default, unless further
    # platforms are passed in on the command line
    platforms = [
        Linux(:x86_64, libc=:glibc)
    ]

    # The products that we will ensure are always built
    products = Product[
        ExecutableProduct("julia", :julia)
    ]

    # Dependencies that must be installed before this package can be built
    dependencies = []

    # Build the tarballs
    product_hashes = cd(joinpath(@__DIR__, "..", "deps")) do
        build_tarballs(binarybuilder_args, "julia", version, sources, script, platforms, products, dependencies, preferred_gcc_version=v"7", skip_audit=true)
    end
    temp_file, hash = product_hashes[platforms[1]]
    tarball = basename(temp_file)

    # Update Versions.toml
    version_stanza = """
        ["$version"]
        file = "$tarball"
        sha = "$hash"
        """
    open(versions_file(); append=true) do f
        println(f, version_stanza)
    end

    # Copy the generated tarball to the downloads folder
    tarball_path = downloads_dir(tarball)
    if ispath(tarball_path)
        # NOTE: we can't use the previous file here (like in `download_julia`)
        #       because the hash will most certainly be different
        @warn "Destination file $tarball already exists, overwriting"
    end
    mv(temp_file, tarball_path)

    return version
end


text:
# kmeans implementation for the ospa distance.

function ospa_kmeans(measurements::Vector{Pointcloud},
                     c::S
                     ;
                     nclusters::Integer=2,
                     niter::Integer=100,
                     eps::Float64=0.05,
                     verbose=false
                    ) where {S <: Real}
    # This has to go into the for loop
    clusters = init_barycenter_clusters(measurements, nclusters)
    costs = fill(Inf, nclusters)
    max_barycenter_rounds = 1
    old_indices = fill(Inf, length(measurements))
    for i=1:niter
        if verbose
            println("Iteration: $(i)")
        end
        dmat = ospa_dist(clusters, measurements, c)
        indices = map(x -> findmin(dmat[:,x])[2], 1:size(dmat)[2])
        for j =1:nclusters
            inds = indices.==j
            # this should be provided by the user or automatically derived
            weights = fill(1/sum(inds), sum(inds))
            clusters[j], costs[j] = ospa_barycenter(measurements[inds], c, weights, :binned)
            if verbose
                println("Cost $(j): $(costs[j])")
            end
        end
        # Mechanism for ending prematurely.
        if old_indices == indices
            break
        end
        old_indices = indices
    end
    clusters, costs
end

function init_barycenter_clusters(datapoints::Vector{Pointcloud}, nclusters::Integer)
    datapoints[[1,end]]
end


text:
using SQLite
co = SQLite.connect(Pkg.dir() * "/SQLite/test/Chinook_Sqlite.sqlite")

df = query("SELECT * FROM Employee;")
createtable(df; name="test")
query("select * from test")
droptable("test")

query("SELECT * FROM sqlite_master WHERE type='table' ORDER BY name;")
query("SELECT * FROM Album;")
query("SELECT a.*, b.AlbumId
	FROM Artist a
	LEFT OUTER JOIN Album b ON b.ArtistId = a.ArtistId
	ORDER BY name;")

using DataFrames
df2 = DataFrame(ones(1000000,5))
@time createtable(df2;name="test2")
@time query("SELECT * FROM test2;")
droptable("test2")

df3 = DataFrame(ones(1000,5))
@time sqldf("select * from df3")

@time readdlmsql(Pkg.dir() * "/SQLite/test/sales.csv";sql="select * from sales",name="sales")
@time query("select typeof(f_year), typeof(base_lines), typeof(dollars) from sales")
droptable("sales")


text:
# Quate Displace
function BioAFMize_QD_FFT(observedAfm, pdbList, quateList, sigma, config)::BioAFMResult
    imgH, imgW = size(observedAfm)
    pdbNum = size(pdbList, 1)
    quateNum = size(quateList, 1)
    posteriors = Float64[]
    mostPropability = Float64(0)
    mostPropabilityPDB = Int(0)
    mostPropabilityQuate = Int(0)

    for pdbId in 1:pdbNum
        sumPropability = Float64(0)
        for quateId in 1:quateNum
            pdbRotated = MDToolbox.rotate(pdbList[pdbId], quateList[quateId, :])
            calculatedAfm = afmize_beta(pdbRotated, config)

            nowPropability = bayesEstimationWithFFTConvolution(observedAfm, calculatedAfm, sigma)

            sumPropability += nowPropability
            if mostPropability < nowPropability
                mostPropability = nowPropability
                mostPropabilityPDB = pdbId
                mostPropabilityQuate = quateId
            end
        end
        push!(posteriors, sumPropability)
    end
    println(posteriors)
    posteriors ./= sum(posteriors)
    return BioAFMResult(posteriors, mostPropabilityPDB, mostPropabilityQuate)
end

max_stars_repo_name: matsunagalab/BioAFM

max_issues_repo_name: matsunagalab/BioAFM

max_forks_repo_name: matsunagalab/BioAFM


text:
module Ledger

export Transaction, Block, Blockchain, gethash, appendblock, checkdifficulty

using Dates
using SHA

struct Transaction
	sender_id::UInt32
	sender_name::String
	recipient_id::UInt32
	recipient_name::String
	quantity::Float32
	message::String
	signature::String
	fee::Float32

end

struct Block
    block_id::UInt32
    timestamp_created::DateTime
    transactions::Vector{Transaction}
    parent_hash::String
    hash::String
    proof_id::String

    function Block(block_id, timestamp_created, transactions, parent_hash, proof_id)
        hash = sha2_256(string(transactions, parent_hash, proof_id))
        new(block_id, timestamp_created, transactions, parent_hash, bytes2hex(hash), proof_id)
    end
end

mutable struct Blockchain
    name::String
    timestamp_created::DateTime
    timestamp_updated::DateTime

    blocks::Vector{Block}

    max_block_size::UInt32
    difficulty_level::UInt32
    mining_reward::Float32

	function Blockchain(name, timestamp_created, timestamp_updated, blocks, max_block_size, difficulty_level, mining_reward)
		new(name, timestamp_created, timestamp_updated, blocks, max_block_size, difficulty_level, mining_reward)
	end

	# Construct genesis block
	function Blockchain(name, max_block_size=1, difficulty_level=4, mining_reward=0)
		timestamp_created = Dates.now()
		first_block = Block(42, timestamp_created, Transaction[], "Genesis block", "69")

		new(name, timestamp_created, timestamp_created, Block[first_block], max_block_size, difficulty_level, mining_reward)
	end

end

# Functions belonging to Blockchain struct

function verifyblock(blockchain::Blockchain, block::Block)
    valid = true
    # verify block size

    # verify hash
    hash = gethash(block)
    valid = checkdifficulty(blockchain, block)

    return valid
end

function gethash(block::Block)
    tx_string = join([x.signature for x in block.transactions])
    hash = bytes2hex(sha2_256(string(block.proof_id) * tx_string * block.parent_hash))
    return hash
end

function checkdifficulty(blockchain::Blockchain, block::Block)
    return startswith(block.hash, "0"^blockchain.difficulty_level)
end

function appendblock(blockchain::Blockchain, block::Block)
    if verifyblock(blockchain, block)
        push!(blockchain.blocks, block)
	return true
    else
	println("WTF dude..")
	return false
    end
end

end

max_stars_repo_name: marchunter/litcoin

max_issues_repo_name: marchunter/litcoin

max_forks_repo_name: marchunter/litcoin


text:
postagg = FieldAccess("name", "field")

@testset "Arithmetic" begin
    @test_throws ErrorException("Unknown arithmetic function") Arithmetic("name", "%", [postagg])
    @test_throws ErrorException("Invalid ordering") Arithmetic("name", "+", [postagg], ordering="a")
end

@testset "Greatest" begin
    @test_throws ErrorException("Invalid data type") Greatest("name", [postagg], "float")
end

@testset "Least" begin
    @test_throws ErrorException("Invalid data type") Least("name", [postagg], "float")
end


text:
#output.jl
outputCount = 0
#Example

function checkEndCon()
	if outputCount > 10
		println("CheckEndCon: End!")
		println(" ")
		return true
	else
		println("CheckEndCon: Nope ", outputCount)
		println(" ")
		return false
	end
end

max_stars_repo_name: voutasaurus/metalbolicalgorithm

max_issues_repo_name: voutasaurus/metalbolicAlgorithm

max_forks_repo_name: voutasaurus/metalbolicAlgorithm


text:
using Test, Colors

@testset "Algorithms" begin
    @test isconcretetype(eltype(colormap("Grays")))

    @test_throws ArgumentError colormap("Grays", N=10)

    col = distinguishable_colors(10)
    @test isconcretetype(eltype(col))
    local mindiff
    mindiff = Inf
    for i = 1:10
        for j = i+1:10
            mindiff = min(mindiff, colordiff(col[i], col[j]))
        end
    end
    @test mindiff > 8
end


text:
# MIT License

# Copyright (c) 2021 Fraunhofer-Chalmers Centre

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

########################################
# Reads the results from the "result" directory (assuming that the CPLEX results from http://soa.iti.es/files/RCmax.7z has been formatted and placed there also)
# output results the data profiles and the result_table.md
# The data tables are generated from the results in results.7z
########################################
using DelimitedFiles, Plots, Statistics, Latexify, Printf
cd(@__DIR__ )
folders = [name  for  name in readdir("results") if isnothing(findfirst('.',name))]
##
function data_profile(n_task=nothing)
    T_lr = []
    T_gur = []
    for folder in folders
        res_path = "results/"*folder*"/res.txt"
        res,header = readdlm(res_path,',',header=true)
        t = res[:,(header[:].=="time")]
        if !isnothing(n_task)
            t = t[res[:,(header[:].=="n")].<=n_task]
        end
        T_lr = vcat(T_lr,t)

        res_path = "results/"*folder*"/res_gur_cut.txt"
        res,header = readdlm(res_path,',',header=true)
        t = res[:,(header[:].=="time")]
        if !isnothing(n_task)
            t = t[res[:,(header[:].=="n")].<=n_task]
        end
        T_gur = vcat(T_gur,t)
    end

    T = [T_lr T_gur]

    τ = 10 .^ range(-2.1, log10(120), length = 100)
    n, n_solv = size(T)
    ρ =  [[count(T[:,i] .< t) for t in τ] ./ n for i in 1:n_solv ]
    plotly()
    title_label = "Data profiles"
    if !isnothing(n_task)
        title_label *= ", n ≤ "*string(n_task)
    end
    plot(τ, 100*ρ, ylims=[0,100], title=title_label, label = ["LR" "Gur"], xlabel="CPU time [s]", ylabel ="Solved [%]", xaxis=:log, xticks = 10.0 .^(-2:2),legend=:bottomright)
end
plt = data_profile()
savefig(plt, "data_profile.png")
display(plt)

plt = data_profile(200)
savefig(plt, "data_profile200.png")
display(plt)
##
function getvalue(res,header,instance,label)
    row = findfirst(res[:,1].==instance)
    if isnothing(row)
        return NaN
    end
    return res[row, header[:].==label]
end
function markbest(z,zmin)
    if z == zmin
        return "**"*string(z)*"**"
    else
        return z
    end
end
function calcmean(res, header, m, n, label)
    rows = (res[:,(header[:].=="m")] .== m) .& (res[:,(header[:].=="n")] .==n)
    if !any(rows)
        return 0
    end
    return  mean(res[rows[:],(header[:].==label)])
end
function n_unsolved(res, header, m, n)
    rows = (res[:,(header[:].=="m")] .== m) .& (res[:,(header[:].=="n")] .==n)
    if !any(rows)
        return 0
    end
    return Int(sum(1 .-res[rows[:],(header[:].=="solved")])  )
end

function data_tables()
    # for each size compute a mean time, iter, nodes, solved, obj
    machines = [10, 20,30,40,50]
    tasks = [100, 200,500,1000]
    outfile = "result_tables.md"
    io = open(outfile,"w")
    write(io, "#Result summary tables\n")
    write(io, "f <- solution not proved to be optimal\n")
    for folder in folders
        res_path = "results/"*folder*"/res.txt"

        out_table =["m" "n" "LR [s]" "LR #f" "LR iter" "LR nodes" "LR z" "LR root gap [%]" "Gur [s]" "Gur #f" "Gur iter" "Gur nodes" "Gur z" "Gur root gap [%]" "Gur gap [%]" "Cplex z"]

        LRres, LRheader = readdlm(res_path,',',header=true)
        res_path = "results/"*folder*"/res_gur_cut.txt"
        GURres, GURheader = readdlm(res_path,',',header=true)
        CplexRes, CplexHeader = readdlm("results/"*folder*"/Cplex2horas"*folder*".csv",';',header=true)
        CplexHeader = [CplexHeader "n" "m"]
        CplexRes = Any[CplexRes Int.(floor.(CplexRes[:,1] .-1,sigdigits=1)) Int.(floor.(CplexRes[:,1].-1;digits=-1)) .% 100]

        for m in machines
            for n in tasks
                LRtime = @sprintf("%0.2f", calcmean(LRres, LRheader, m, n, "time"))
                LRunsolved = n_unsolved(LRres, LRheader, m, n)
                LRiter = round(Int,calcmean(LRres, LRheader, m, n, "iter"))
                LRnodes = round(Int,calcmean(LRres, LRheader, m, n, "nodes"))
                LRz = calcmean(LRres, LRheader, m, n, "z")
                LRrootgap = @sprintf("%0.2f", calcmean(LRres, LRheader, m, n, "lbgaproot")*100)
                Gurtime =  @sprintf("%0.2f",calcmean(GURres, GURheader, m, n, "time"))
                Gurunsolved = n_unsolved(GURres, GURheader, m, n)
                Guriter =round(Int, calcmean(GURres, GURheader, m, n, "iter"))
                Gurnodes =round(Int, calcmean(GURres, GURheader, m, n, "nodes"))
                Gurz = round(calcmean(GURres, GURheader, m, n, "z"), digits=1)
                Gurrootgap = @sprintf("%0.2f", calcmean(GURres, GURheader, m, n, "lbrootgap")*100)
                Gurrgap =@sprintf("%0.2f",  calcmean(GURres, GURheader, m, n, "gap")*100  )
                Cplexz = calcmean(CplexRes, CplexHeader, m, n, "Objective")
                out_table = vcat(out_table, [m n LRtime LRunsolved LRiter LRnodes LRz LRrootgap Gurtime Gurunsolved Guriter Gurnodes Gurz Gurrootgap Gurrgap Cplexz])
            end
        end
        write(io, "## "*folder*"\n\n")
        write(io,  string(md(out_table, latex=false)))
        write(io, "\n")
    end

    out_table =["type" "instance" "LR value" "LR opt" "Gurobi value" "Gurobi opt" "Cplex value" "Cplex opt"]
    for folder in folders
        res_path = "results/"*folder*"/res.txt"
        LRres, LRheader = readdlm(res_path,',',header=true)
        res_path = "results/"*folder*"/res_gur_cut.txt"
        GURres, GURheader = readdlm(res_path,',',header=true)
        GURres = round.(Int,GURres)
        CplexRes, CplexHeader = readdlm("results/"*folder*"/Cplex2horas"*folder*".csv",';',header=true)
        CplexRes = round.(Int,CplexRes)
        for instance in LRres[:,1]
            LRz = getvalue(LRres,LRheader,instance,"z")
            LRopt = getvalue(LRres,LRheader,instance,"solved")
            GURz = getvalue(GURres,GURheader,instance,"z")
            Guropt = getvalue(GURres,GURheader,instance,"solved")
            Cplexz = getvalue(CplexRes,CplexHeader,instance,"Objective")
            Cplexopt = getvalue(CplexRes,CplexHeader,instance,"Proven opt.")
            out_table = vcat(out_table,[folder instance LRz LRopt GURz Guropt Cplexz Cplexopt] )
        end
    end

    write(io, "## All instances summary \n\n")
    A = copy(out_table[2:end,3:end])
    A[isnan.(A)].=Inf
    z = A[:,1:2:end]
    best_z = Int.(minimum(z, dims=2))
    has_best =  z .<= best_z
    has_strict_best = reduce(hcat,[z[:,j] .< minimum(z[:,setdiff(1:3,j)],dims=2) for j in 1:3])
    has_opt = A[:,2:2:end] .== 1
    unique_opt = [has_opt[i,j] && sum(has_opt[i,:])==1  for i in 1:size(has_opt,1),j in 1:3]

    write(io,  string(md( [["" "LR" "Gurobi" "Cplex"];
     hcat("Best sol [%]", round.(100*sum(has_best,dims=1)/size(has_best,1),digits=2));
     hcat("Strict best sol [%]", round.(100*sum(has_strict_best,dims=1)/size(has_best,1),digits=2));
     hcat("Proved opt [%]",round.(100*sum(has_opt,dims=1)/size(has_opt,1),digits=2));
     hcat("Uniquely proved opt [%]",round.(100*sum(unique_opt,dims=1)/size(has_opt,1),digits=2))], latex=false)))

    writedlm("res_all_instances.csv", out_table,",")

    write(io, "## All instances \n\n")
    for (i, j) in Tuple.(findall(has_best))
        out_table[i+1,1+2*j] = "**"*string( out_table[i+1,1+2*j] )*"**"

    end

    write(io,  string(md(out_table, latex=false)))

    close(io)

end
data_tables()

max_stars_repo_name: Fraunhofer-Chalmers-Centre/RCmax

max_issues_repo_name: Fraunhofer-Chalmers-Centre/RCmax

max_forks_repo_name: Fraunhofer-Chalmers-Centre/RCmax


text:
module ComputationalBiologyCourse

using DataDeps

end  # module ComputationalBiologyCourse


text:
#=
    Figure S3
=#

using Model # To get colours

using CSV
using ColorSchemes
using DataFrames
using DataFramesMeta
using Plots
using StatsPlots

include("../figure_defaults.jl")

##############################################################
## Data
##############################################################

data = @subset(CSV.read("Data/CSV/DualMarker.csv",DataFrame), :Temp .== 37.0, :Time .> 0.0)

##############################################################
## Univariate
##############################################################

T = sort(unique(data.Time))

plts_q = [plot() for t in T]
plts_u = [plot() for t in T]

for (i,t) in enumerate(T)

    @df @subset(data,:Time .== t,(!).(:Quenched)) density!(plts_q[i],:Signal_Cy5 / 1e4,c=col_Q,lw=2.0,frange=0.0,fα=0.4,label="Not quenched")
    @df @subset(data,:Time .== t,:Quenched) density!(plts_q[i],:Signal_Cy5 / 1e4,c=col_Q̄,lw=2.0,frange=0.0,fα=0.4,ls=:dash,label="Quenched")
    plot!(plts_q[i],xlim=(0.0,3.0),widen=:false,yticks=[],xlabel="Cy5 (10⁴)",axis=:x,title="$(Int(t)) min")

    @df @subset(data,:Time .== t,(!).(:Quenched)) density!(plts_u[i],:Signal_BDP,c=col_U,lw=2.0,frange=0.0,fα=0.4,label="Not quenched")
    @df @subset(data,:Time .== t,:Quenched) density!(plts_u[i],:Signal_BDP,c=col_Ū,lw=2.0,frange=0.0,fα=0.4,ls=:dash,label="Quenched")
    plot!(plts_u[i],xlim=(0.0,120.0),widen=:false,yticks=[],xlabel="BDP",axis=:x)

    if i > 1
        plot!(plts_q[i],legend=:none)
        plot!(plts_u[i],legend=:none)
    end

end

figS3 = plot(plts_q...,plts_u...,layout=grid(2,7),size=(800,250))

##############################################################
## Save...
##############################################################

savefig("$(@__DIR__)/FigS3.svg")

max_stars_repo_name: ap-browning/internalisation

max_issues_repo_name: ap-browning/internalisation

max_forks_repo_name: ap-browning/internalisation


text:
# classify the nodes in a tree as constant, linear, or nonlinear with respect to the input

@enum Linearity CONSTANT LINEAR PIECEWISE_LINEAR NONLINEAR

export CONSTANT, LINEAR, PIECEWISE_LINEAR, NONLINEAR

function classify_linearity(nd::Vector{NodeData},adj,subexpression_linearity,fixed_variable = [])

    linearity = Array{Linearity}(undef, length(nd))

    # do a forward pass through the graph, which is reverse order of nd

    children_arr = rowvals(adj)

    for k in length(nd):-1:1
        # compute the value of node k
        nod = nd[k]
        if nod.nodetype == VARIABLE
            if length(fixed_variable) > 0 && fixed_variable[nod.index]
                linearity[k] = CONSTANT
            else
                linearity[k] = LINEAR
            end
        elseif nod.nodetype == VALUE || nod.nodetype == PARAMETER
            linearity[k] = CONSTANT
        elseif nod.nodetype == SUBEXPRESSION
            linearity[k] = subexpression_linearity[nod.index]
        else
            children_idx = nzrange(adj,k)
            num_constant_children = 0
            any_nonlinear = false
            for r in children_idx
                if linearity[children_arr[r]] == NONLINEAR
                    any_nonlinear = true
                    break
                elseif linearity[children_arr[r]] == CONSTANT
                    num_constant_children += 1
                end
            end
            if any_nonlinear
                # if any children are nonlinear, then we're nonlinear
                linearity[k] = NONLINEAR
                # except in the case of ifelse. if the operands are linear then we're piecewise linear
                if nod.nodetype == CALL && operators[nod.index] == :ifelse
                    if linearity[children_arr[children_idx[2]]] == LINEAR && linearity[children_arr[children_idx[3]]] == LINEAR
                        linearity[k] = PIECEWISE_LINEAR
                    end
                end
                continue
            elseif num_constant_children == length(children_idx)
                # if all children are constant, then we're constant
                linearity[k] = CONSTANT
                continue
            end
            # some children are constant and some are linear
            # if operator is nonlinear, then we're nonlinear
            op = nod.index
            if nod.nodetype == CALLUNIVAR
                if op < USER_UNIVAR_OPERATOR_ID_START && (univariate_operators[op] == :+ || univariate_operators[op] == :-)
                    linearity[k] = LINEAR
                else
                    linearity[k] = NONLINEAR
                end
            elseif nod.nodetype == CALL
                # operator with more than 1 argument
                if op >= USER_OPERATOR_ID_START
                    linearity[k] = NONLINEAR
                elseif operators[op] == :+ || operators[op] == :-
                    linearity[k] = LINEAR
                elseif operators[op] == :* && num_constant_children == length(children_idx) - 1
                    linearity[k] = LINEAR
                elseif operators[op] == :/
                    if linearity[children_arr[last(children_idx)]] == CONSTANT
                        # denominator is constant, we're linear
                        linearity[k] = LINEAR
                    else
                        # denominator is linear
                        linearity[k] = NONLINEAR
                    end
                else # all other operators are nonlinear
                    linearity[k] = NONLINEAR
                end
            elseif nod.nodetype == LOGIC || nod.nodetype == COMPARISON
                linearity[k] = NONLINEAR
            else
                error("Unknown node type")
            end
        end
    end

    return linearity

end

export Linearity,classify_linearity


text:
alpha = 0;
m = 3;

# ground truth solution
function u(x1, x2)
    return sin(pi*x1)*sin(pi*x2) + 2*sin(4*pi*x1)*sin(4*pi*x2)
end

# right hand side
function f(x1, x2)
    return 2*pi^2*sin(pi*x1)*sin(pi*x2)+64*pi^2*sin(4*pi*x1)*sin(4*pi*x2)+alpha*u(x1,x2)^m
end

# boundary value
function g(x1, x2)
    return u(x1, x2)
end


text:
module Day05

using OffsetArrays

export day05

function parseLine(line)
    raw = split.(split(line, " -> "), ",")
    return complex(parse.(Int64, raw[1])...), complex(parse.(Int64, raw[2])...)
end

function readInput(inputfile="inputs/input05.txt")
    return [parseLine(line) for line in eachline(inputfile)]
end

function day05(ventdata = readInput("inputs/input05.txt"))
    minx = minimum((a) -> min(real(a[1]), real(a[2])), ventdata)
    miny = minimum((a) -> min(imag(a[1]), imag(a[2])), ventdata)
    maxx = maximum((a) -> max(real(a[1]), real(a[2])), ventdata)
    maxy = maximum((a) -> max(imag(a[1]), imag(a[2])), ventdata)

    ## previous version used a Dict here
    ## that's ~25x slower
    grid = OffsetArray(zeros(UInt8, maxx-minx+1, maxy-miny+1), minx:maxx, miny:maxy)

    mask = ishorizontal.(ventdata)
    mapvents!(grid, ventdata[mask])
    part1 = count(>=(2), grid)

    mapvents!(grid, ventdata[.!(mask)])
    part2 = count(>=(2), grid)
    return part1, part2
end

function mapvents!(grid, data)
    for (start, stop) in data
        Δ = stop-start
        dist = distance(Δ)
        δ = intdiv(Δ, dist)
        for l in 0:dist
            pos = start+δ*l
            grid[real(pos), imag(pos)] += 1
        end
    end
    return grid
end

ishorizontal((from, to)) = real(from) == real(to) || imag(from) == imag(to)
distance(x::Complex) = let r = abs(real(x)); r == 0 ? abs(imag(x)) : r end
intdiv(x::Complex, by::Int) = complex(div(real(x), by), div(imag(x), by))

count_crossings(grid) = count(>=(2), values(grid))

end #module


text:
handle("updateeditor") do data
    @destruct [
        text || "",
        mod || "Main",
        path || nothing,
        updateSymbols || true
    ] = data

    try
        updateeditor(text, mod, path, updateSymbols)
    catch err
        []
    end
end

# NOTE: update outline and symbols cache all in one go
function updateeditor(text, mod = "Main", path = nothing, updateSymbols = true)
    # update symbols cache
    # ref: https://github.com/JunoLab/Juno.jl/issues/407
    updateSymbols && updatesymbols(mod, path, text)

    # return outline
    outline(toplevelitems(text))
end

function outline(items)
    filter!(map(outlineitem, items)) do item
        item !== nothing
    end
end

function outlineitem(binding::ToplevelBinding)
    expr = binding.expr
    bind = binding.bind
    lines = binding.lines

    name = bind.name
    if CSTParser.has_sig(expr)
        name = str_value(CSTParser.get_sig(expr))
    end
    type = static_type(bind)
    icon = static_icon(bind)
    Dict(
        :name  => name,
        :type  => type,
        :icon  => icon,
        :lines => [first(lines), last(lines)]
    )
end
function outlineitem(call::ToplevelCall)
    expr = call.expr
    lines = call.lines

    # don't show doc strings
    isdoc(expr) && return nothing

    # show first line verbatim of macro calls
    if ismacrocall(expr)
        return Dict(
            :name  => strlimit(first(split(call.callstr, '\n')), 100),
            :type  => "snippet",
            :icon  => "icon-mention",
            :lines => [first(lines), last(lines)],
        )
    end

    # show includes
    if isinclude(expr)
        return Dict(
            :name  => call.callstr,
            :type  => "module",
            :icon  => "icon-file-code",
            :lines => [first(lines), last(lines)],
        )
    end

    return nothing
end
function outlineitem(tupleh::ToplevelTupleH)
    expr = tupleh.expr
    lines = tupleh.lines

    # `expr.parent` is always `CSTParser.EXPR`
    type = isconstexpr(expr.parent) ? "constant" : "variable"
    icon = isconstexpr(expr.parent) ? "c" : "v"

    Dict(
        :name  => str_value(expr),
        :type  => type,
        :icon  => icon,
        :lines => [first(lines), last(lines)]
    )
end


text:
# This is the core functionality for manipulating dimension names as express as a tuple
# It is central to the performance of this named dimensions as a zero cost abstraction
# For constant input all functions should constant propagate to constant output.
# Each function in this file is annotated with a comment and a benchmark
# running that benchmark should return 0 allocations,
# and you can swap its `@btime` for `@code_typed` to check if it is constant propergating
# This constant propagate lets us write faily naive code elsewhere in the package
# knowing the creation and destruction of NamedDimsArrays will be optimised away.

"""
    compile_time_return_hack(x)

This is a cunning hack to strongly encourage the compiler to use constant propagation,
when returning a `Tuple` of `Symbols`.
Under normal circumstances, constant propergation will not work to fully compute
returned tuples of non-bits types.
But because `Symbol`s are a allowed in type-parameters (unlike other nonbits types),
we can put them into a type-parameter, which triggers (??different??) constant propergation,
(for ??Reasons??)
and once something is a type-parameter, then extracting it and returning it is non-allocating.

See https://discourse.julialang.org/t/zero-allocation-tuple-subsetting/23122/8

The short version of this is:
if a function returns a tuple of symbols and is allocating when it looks like it shouldn't
then try wrapping its returned value in `compile_time_return_hack`.
and then look at the `@code_lowered` again.
"""
compile_time_return_hack(x::Tuple{Vararg{Symbol}}) = _compile_time_return_hack(Val{x}())
_compile_time_return_hack(::Val{X}) where X = X

"""
    dim(dimnames, name)

For `dimnames` being a tuple of names (symbols) for the dimensions.
and `name` being a one of those names
This returns the dimension corresponding to that `name`,
e.g. `dim((:a, :b), :b) == 2`
If that `name` is not one of the given `dimnames` then an error is thrown.
"""
function dim(dimnames::Tuple, name::Symbol)::Int
    # 0-Allocations see: `@btime  (()->dim((:a, :b), :b))()`
    dimnum = dim_noerror(dimnames, name)
    if dimnum === 0
        throw(ArgumentError(
            "Specified name ($(repr(name))) does not match any dimension name ($dimnames)"
        ))
    end
    return dimnum
end

function dim(dimnames::Tuple, names)
    # 0-Allocations see: `@btime (()->dim((:a,:b), (:a,:b)))()`
    return map(name->dim(dimnames, name), names)
end

function dim(dimnames::Tuple, d::Union{Integer, Colon})
    # This is the fallback that allows `NamedDimsArray`'s to be have dimensions
    # referred to by number. This is required to allow functions on `AbstractArray`s
    # and that use function like `sum(xs; dims=2)` to continue to work without changes
    # `:` is the default for most methods that take `dims`
    return d
end

Base.@pure function dim_noerror(dimnames::Tuple{Vararg{Symbol, N}}, name::Symbol) where N
    # 0-Allocations see: @btime  (()->dim_noerror((:a, :b, :c), :c))()
    for ii in 1:N
        getfield(dimnames, ii) === name && return ii
    end
    return 0
end

"""
    permute_dimnames(dimnames, perm)

Reorder `dimnames` according `perm`.
`perm` should be ordered set of numerical indexs for the new position of the name.
Note: this does not throw errors if you give it a permutation that skips some positions
and duplicates others.
"""
function permute_dimnames(dimnames::NTuple{N, Symbol}, perm) where N
    # 0-Allocations, but does not seem to fully calculate at compile time
    # even with the `compile_time_return_hack`, though that is still required to
    # prevent allocations. `@code_typed permute_dimnames((:a,:b,:c), (1,3,2))`

    new_dimnames = ntuple(length(perm)) do ii
        ind = perm[ii]
        dimnames[ind]
    end
    return compile_time_return_hack(new_dimnames)
end

"""
    tuple_issubset
A version of `is_subset` sepecifically for `Tuple`s of `Symbol`s, that is `@pure`.
This helps it get optimised out of existance. It is less of an abuse of `@pure` than
most of the stuff for making `NamedTuples` work.
"""
Base.@pure function tuple_issubset(lhs::Tuple{Vararg{Symbol,N}}, rhs::Tuple{Vararg{Symbol,M}}) where {N,M}
    N <= M || return false
    for a in lhs
        found = false
        for b in rhs
            found |= a === b
        end
        found || return false
    end
    return true
end

"""
    order_named_inds(dimnames::Tuple; named_inds...)

Returns the values of the `named_inds`, sorted as per the order they appear in `dimnames`,
with any missing dimnames, having there value set to `:`.
An error is thrown if any dimnames are given in `named_inds` that do not occur in `dimnames`.
"""
order_named_inds(dimnames::Tuple; named_inds...) = order_named_inds(dimnames, named_inds.data)

function order_named_inds(dimnames::Tuple{Vararg{Symbol,N}}, named_inds::NamedTuple) where {N}
    if !tuple_issubset(keys(named_inds), dimnames)
        throw(DimensionMismatch("Expected subset of $(dimnames), got $(keys(named_inds))"))
    end

    full_inds = ntuple(N) do ii
        name_ii = dimnames[ii]
        get(named_inds, name_ii, :)
    end

    return full_inds
end

"""
    incompatible_dimension_error(names_a, names_b)
Throws a `DimensionMismatch` explaining that these dimension names are not compatible.
"""
function incompatible_dimension_error(names_a, names_b)
    throw(DimensionMismatch("Incompatible dimension names: $names_a ≠ $names_b"))
end

"""
    unify_names(a, b)

Produces the merged set of names for tuples of names `a` and `b`,
or an error if it is not possibly to unify them.
Two tuples of names can be unified they are the same length
and if for each position the names are either the same, or one is a wildcard (`:_`).
When combining wildcard with non-wildcard the resulting name is the non-wildcard.
(This is somewhat like the very simplest case of unification in e.g prolog).

For example:
 - `(:a, :b)` and `(:a, :b)` can be unified to give `(:a, :b)`
 - similarly: `(:a, :_)` and `(:a, :b)` can also be unified to give `(:a, :b)`
 - `(:a, :b)` and `(:b, :c)` cannot be unified as the names at each position to not match.
 - `(:a, :b)` and `(:a, :b, :c)` cannot be unified as they have different lengths

This is the type of name combination used for binary array operations.
Where the dimensions of both arrays must be the same.
"""
function unify_names(names_a, names_b)
    # 0-Allocations if inputs are the same
    # 0-Allocation, if has a `:_` see  `@btime (()->unify_names((:a, :b), (:a, :_)))()`

    names_a === names_b && return names_a

    # Error message should not include names until it is thrown, as othrwise
    # the interpolation allocates and slows everything down a lot.
    length(names_a) != length(names_b) && incompatible_dimension_error(names_a, names_b)

    ret = ntuple(length(names_a)) do ii  # remove :_ wildcards
        a = getfield(names_a, ii)
        b = names_b[ii]
        a === :_ && return b
        b === :_ && return a
        a === b && return a

        return false  # mismatch occured, we mark this with a nonSymbol result
    end
    ret isa Tuple{Vararg{Symbol}} || incompatible_dimension_error(names_a, names_b)
    return compile_time_return_hack(ret)
end

unify_names(a) = a
unify_names(a, b, cs...) = unify_names(unify_names(a,b), cs...)
# @btime (()->unify_names((:a, :b), (:a, :_), (:_, :b)))()

"""
    unify_names_longest(a, b)

This is the same as [`unify_names`](@ref), but with the equal length requirement removed.
It unifies the names up to the length of the shortest, and takes the named from the longest
for the remainder.
It can also be considered as padding the shorter of the two given tuples of names
with trailing wildcards (`:_`).
This is the type of name combination used for broadcating array operations.
Where the smaller (dimensionally) array is broadcast against the longer, repeating it for
all entries in the missing trailing dimensions.
"""
unify_names_longest(names, ::Tuple{}) = names
unify_names_longest(::Tuple{}, names) = names
unify_names_longest(::Tuple{}, ::Tuple{}) = tuple()
function unify_names_longest(names_a, names_b)
    # 0 Allocations: @btime (()-> unify_names_longest((:a,:b), (:a,)))()

    length(names_a) == length(names_b) && return unify_names(names_a, names_b)
    long, short = length(names_a) > length(names_b) ? (names_a, names_b) : (names_b, names_a)
    ret = ntuple(length(long)) do ii
        a = getfield(long, ii)
        b = ii <= length(short) ? short[ii] : :_
        a === :_ && return b
        b === :_ && return a
        a === b && return a
        return false  # mismatch occured, we mark this with a nonSymbol result
    end
    ret isa Tuple{Vararg{Symbol}} || incompatible_dimension_error(names_a, names_b)
    return compile_time_return_hack(ret)
end

unify_names_shortest(names, ::Tuple{}) = ()
unify_names_shortest(::Tuple{}, names) = ()
unify_names_shortest(::Tuple{}, ::Tuple{}) = ()
function unify_names_shortest(names_a, names_b)
    # 0 Allocations: @btime (()-> unify_names_shortest((:a,:b), (:a,)))()

    length(names_a) == length(names_b) && return unify_names(names_a, names_b)
    long, short = length(names_a) > length(names_b) ? (names_a, names_b) : (names_b, names_a)
    ret = ntuple(length(short)) do ii
        a = getfield(long, ii)
        b = getfield(short, ii)
        a === :_ && return b
        b === :_ && return a
        a === b && return a
        return false  # mismatch occured, we mark this with a nonSymbol result
    end
    ret isa Tuple{Vararg{Symbol}} || incompatible_dimension_error(names_a, names_b)
    return compile_time_return_hack(ret)
end

# The following are helpers for remaining_dimnames_from_indexing
# as a generated function it can get unhappy if asked to use anon functions
# and it can only call function declared before it. So we declare them explictly here.
is_noninteger_type(::Type{<:Integer}) = false
is_noninteger_type(::Any) = true

"""
    remaining_dimnames_from_indexing(dimnames::Tuple, inds)
Given a tuple of dimension names
and a set of index expressesion e.g `1, :, 1:3, [true, false]`,
determine which are not dropped.
Dimensions indexed with scalars are dropped
"""
@generated function remaining_dimnames_from_indexing(dimnames::Tuple, inds)
    # 0-Allocation see:
    # `@btime (()->remaining_dimnames_from_indexing((:a, :b, :c), (:,390,:)))()``
    keep_names = []
    dim_num = 0
    for type in inds.parameters
        if type <: Integer
            dim_num += 1
        elseif type <: CartesianIndex
            dim_num += type.parameters[1]
        elseif type == Array{CartesianIndex{0},1}
            push!(keep_names, QuoteNode(:_))
        else
            dim_num += 1
            push!(keep_names, :(getfield(dimnames, $dim_num)) )
        end
    end
    return Expr(:call, :compile_time_return_hack, Expr(:tuple, keep_names...))
end

"""
    remaining_dimnames_after_dropping(dimnames::Tuple, dropped_dims)
Given a tuple of dimension names, and either a collection of dimensions,
or a single dimension, expressed as a number,
Returns the dimension names with those dimensions dropped.
"""
function remaining_dimnames_after_dropping(dimnames::Tuple, dropped_dim::Int)
    # 0 allocations. See `@btime remaining_dimnames_after_dropping((:a,:b,:c,:d,:e), 4)`
    return _remaining_dimnames_after_dropping(dimnames, Tuple{dropped_dim})
end

function remaining_dimnames_after_dropping(dimnames::NTuple{N, Symbol}, dropped_dims::Tuple{Vararg{Int}}) where N
    # 0-Allocations see:
    # `@code_typed (()->remaining_dimnames_after_dropping((:a,:b,:c,:d,:e), (1,3)))()`
    return _remaining_dimnames_after_dropping(dimnames, Tuple{dropped_dims...})
end

# dropped_dims must be a Tuple type where the values are Int literal for the dimensions being dropped
@generated function _remaining_dimnames_after_dropping(dimnames::NTuple{N,Symbol}, dropped_dims::Type) where N
    dropped_dims_vals = dropped_dims.parameters[1].parameters
    keep_names = [:(getfield(dimnames, $ii)) for ii in 1:N if ii ∉ dropped_dims_vals]
    return Expr(:call, :compile_time_return_hack, Expr(:tuple, keep_names...))
end

"""
    tuple_cat(x, y, zs...)

This is like `vcat` for tuples, it splats everything into one long tuple.
"""
tuple_cat(x::Tuple, ys::Tuple...) = (x..., tuple_cat(ys...)...)
tuple_cat() = ()
# @btime tuple_cat((1, 2), (3, 4, 5), (6,)) # 0 allocations


text:
struct CoordMode{D} <: IndexMode
    dims::D
end

dims(m::CoordMode) = m.dims
order(m::CoordMode) = Unordered()

"""
    Coord <: Dimension

A coordinate dimension itself holds dimensions.

This allows combining point data with other dimensions, such as time.

# Example

```julia
julia> using DimensionalData

julia> dim = Coord([(1.0,1.0,1.0), (1.0,2.0,2.0), (3.0,4.0,4.0), (1.0,3.0,4.0)], (X(), Y(), Z()))
Coord ::
  val: Tuple{Float64, Float64, Float64}[(1.0, 1.0, 1.0), (1.0, 2.0, 2.0), (3.0, 4.0, 4.0), (1.0, 3.0,
4.0)]
  mode: CoordMode
Coord{Vector{Tuple{Float64, Float64, Float64}}, DimensionalData.CoordMode{Tuple{X{Colon, AutoMode{Auto
Order}, NoMetadata}, Y{Colon, AutoMode{AutoOrder}, NoMetadata}, Z{Colon, AutoMode{AutoOrder}, NoMetada
ta}}}, NoMetadata}

julia> da = DimArray(0.1:0.1:0.4, dim)
4-element DimArray{Float64,1} with dimensions:
  Coord (): Tuple{Float64, Float64, Float64}[(1.0, 1.0, 1.0), (1.0, 2.0, 2.0), (3.0, 4.0, 4.0), (1.0,
3.0, 4.0)]
    CoordMode
 0.1
 0.2
 0.3
 0.4

julia> da[Coord(Z(At(1.0)), Y(Between(1, 3)))]
1-element DimArray{Float64,1} with dimensions:
  Coord (): Tuple{Float64, Float64, Float64}[(1.0, 1.0, 1.0)] CoordMode
 0.1

julia> da[Coord(4)] == 0.4
true

julia> da[Coord(Between(1, 5), :, At(4.0))]
2-element DimArray{Float64,1} with dimensions:
  Coord (): Tuple{Float64, Float64, Float64}[(3.0, 4.0, 4.0), (1.0, 3.0, 4.0)] CoordMode
 0.3
 0.4
```
"""
struct Coord{T,Mo<:IndexMode,Me<:AllMetadata} <: Dimension{T,Mo,Me}
    val::T
    mode::Mo
    metadata::Me
end
function Coord(val::T, dims::Tuple, metadata::Me=NoMetadata()) where {T<:AbstractVector,Me<:AllMetadata}
    length(dims) = length(first(val))
    mode = CoordMode(key2dim(dims))
    Coord{T,typeof(mode),Me}(val, mode, metadata)
end
const SelOrStandard = Union{Selector,StandardIndices}
Coord(s1::SelOrStandard, s2::SelOrStandard, sels::SelOrStandard...) = Coord((s1, s2, sels...))
Coord(sel::Selector) = Coord((sel,))
Coord(d1::Dimension, dims::Dimension...) = Coord((d1, dims...))
Coord(val::T=:) where T = Coord{T,AutoMode,NoMetadata}(val, AutoMode(), NoMetadata())

dims(d::Coord) = dims(mode(d))
bounds(d::Coord) = ntuple(i -> extrema((x[i] for x in val(d))), length(first(d)))

# Return a Vector{Bool} for matching coordinates
sel2indices(dim::Coord, sel::DimTuple) = sel2indices(dim, sortdims(sel, dims(dim)))
sel2indices(dim::Coord, sel::Tuple) = [all(map(_matches, sel, x)) for x in val(dim)]
sel2indices(dim::Coord, sel::StandardIndices) = sel

_matches(sel::Dimension, x) = _matches(val(sel), x)
_matches(sel::Between, x) = (x >= first(sel)) & (x < last(sel))
_matches(sel::At, x) = x == val(sel)
_matches(sel::Colon, x) = true
_matches(sel::Nothing, x) = true

_tozero(xs) = map(x -> zero(x), xs)
@inline _reducedims(mode::CoordMode, dim::Dimension) =
    rebuild(dim, [_tozero(dim.val[1])], dim.mode)


text:
using Test
using NormalHermiteSplines
using NormalHermiteSplines: ElasticCholesky, ElasticNormalSpline

using DoubleFloats
using LinearAlgebra
using Random
using StaticArrays
using UnPack

@testset "NormalHermiteSplines.jl" begin
    include("1D.jl")
    include("2D.jl")
    include("3D.jl")
    include("elastic.jl")
end


text:
function initialise_Q(gp::GPBase)
    # TODO: Use PDMats for the below
    V = cov(gp.kernel, gp.x, gp.data)
    Ω = inv(V)
    K = deepcopy(Ω)
    m = mean(gp.mean, gp.x)
    Q = Approx(m, V)
    return Q, V, K
end

function update_Q!(Q, m, V)
    Q.m = m
    Q.V = V
end

function elbo(y::AbstractArray, μ::AbstractArray, Ω::AbstractMatrix, m::AbstractArray, V::AbstractMatrix, ll::Likelihood)
    @assert length(μ) == length(m)
    @assert size(V) == size(Ω)
    println("y: ", size(y))
    println("μ: ", size(μ))
    println("Ω: ", size(Ω))
    println("m: ", size(m))
    println("V: ", size(V))
    N = length(y)
    VprodΩ = V * Ω
    vexp =  var_exp(ll, y, m, V)
    return 0.5*(logdet(VprodΩ) - tr(VprodΩ) - transpose(m - μ)*Ω*(m - μ) + N) + vexp
end

function push_back(mat::AbstractMatrix, idx::Integer)
    # TODO: Possibly more efficient way to do this.
    ori_size = size(mat)
    target = mat[:, idx]
    temp =  mat[:,setdiff(1:end, idx)]
    mat = cat(temp, target, dims=2)
    return mat
end

function push_back(mat::AbstractVector)
    target = mat[1]
    append!(mat[2:end], target)
    return mat
end

function vi(gp::GPBase; nits::Int64=100)
    function objective_func(params)
        n = gp.nobs
        m = params[1:n]
        V = params[(n + 1):end]
        prod_term = V .* ω
        obj = -0.5*(logdet(prod_term .* Array{Float64}(I, n, n)) - sum(prod_term) - ((m - gp.μ)' * Ω * (m - gp.μ)) + n + var_exp(gp.lik, gp.y, m, V))
        return obj
    end

    # Initialise log-target and log-target's derivative
    mcmc(gp; nIter=1);

    # Initialise variational approximation
    Q, Ω, K = initialise_Q(gp)
    ω = diag(Ω)

    θ = cat(Q.m, diag(Q.V), dims=1)

    # res = optimize(objective_func, θ, Newton())# LBFGS(); autodiff = :forward)
    lower = cat(repeat([-Inf], gp.nobs), repeat([1e-10], gp.nobs), dims = 1)
    upper = cat(repeat([Inf], gp.nobs), repeat([Inf], gp.nobs), dims = 1)
    initial_x = cat(Q.m, diag(Q.V), dims=1)
    inner_optimizer = ConjugateGradient()
    res = optimize(objective_func, lower, upper, initial_x, Fminbox(inner_optimizer))

    optθ = Optim.minimizer(res)
    Q.m = optθ[1:gp.nobs]
    Q.V = (optθ[(gp.nobs + 1):end] .* Array{Float64}(I, gp.nobs, gp.nobs))  + Array{Float64}(I, gp.nobs, gp.nobs)*1e-10

    return Q
end

function predict_f(gp::GPBase, x::AbstractMatrix, Q::Approx)
    size(x,1) == gp.dim || throw(ArgumentError("Gaussian Process object and input observations do not have consistent dimensions"))
    return predict_full(gp, x, Q)
end

# Sample from the variational GP
function Random.rand!(gp::GPBase, x::AbstractMatrix, A::DenseMatrix, Q::Approx)
    nobs = size(x, 2)
    n_sample = size(A, 2)
    μ, Σraw = predict_f(gp, x, Q)
    Σraw, chol = make_posdef!(Σraw)
    Σ = PDMat(Σraw, chol)
    return broadcast!(+, A, μ, unwhiten!(Σ,randn(nobs, n_sample)))
end

# Generate samples from the variational approxiamtion
Random.rand(gp::GPBase, x::AbstractMatrix, n::Int, Q::Approx) = rand!(gp, x, Array{Float64}(undef, size(x, 2), n), Q)
Random.rand(gp::GPBase, x::AbstractVector, Q::Approx) = vec(rand(gp, x', 1, Q))
Random.rand(gp::GPBase, x::AbstractMatrix, Q::Approx) = vec(rand(gp, x, 1, Q))

predict_full(gp::GPA, xpred::AbstractMatrix, Q::Approx) = predictMVN(gp,xpred, gp.x, gp.y, gp.kernel, gp.mean, gp.v, gp.covstrat, Q)

"""
    predict_y(gp::GPA, x::Union{Vector{Float64},Matrix{Float64}}[; full_cov::Bool=false])
Return the predictive mean and variance of Gaussian Process `gp` at specfic points which
are given as columns of matrix `x`. If `full_cov` is `true`, the full covariance matrix is
returned instead of only variances.
"""
function predictMVN(gp::GPBase,xpred::AbstractMatrix, xtrain::AbstractMatrix, ytrain::AbstractVector,
                   kernel::Kernel, meanf::Mean, alpha::AbstractVector,
                   covstrat::CovarianceStrategy, Q::Approx)
    Ktrain = PDMat(Q.V)
    crossdata = KernelData(kernel, xtrain, xpred)
    priordata = KernelData(kernel, xpred, xpred)
    Kcross = cov(kernel, xtrain, xpred, crossdata)
    Kpred = cov(kernel, xpred, xpred, priordata)
    mx = mean(meanf, xpred)
    mu, Sigma_raw = predictMVNvi!(gp,Kpred, Ktrain, Kcross, mx, alpha) #TODO: Handle through multiple dispatch
    return mu, Sigma_raw
end

function predictMVN!(gp::GPA, Kxx, Kff, Kfx, mx, αf)
    Lck = whiten!(Kff, Kfx)
    mu = mx + Lck' * αf
    subtract_Lck!(Kxx, Lck)
    return mu, Kxx
end

function predictMVNvi!(gp::GPA, Kxx, Kff, Kfx, mx, αf)
    Lck = whiten!(Kff, Kfx)
    mu = mx + Lck' * αf
    return mu, Kxx
end


text:
using Documenter, SparseGaussianProcesses

makedocs(
  sitename = "SparseGaussianProcesses.jl",
  modules = [SparseGaussianProcesses],
  pages = [
    "Home"=>"index.md",
    "API"=>"api.md"
  ])

deploydocs(repo = "github.com/aterenin/SparseGaussianProcesses.jl.git")


text:
using CSV, DataFrames, DataFramesMeta, Base.Dates, Plots
using JuMP, Gurobi, Distributions, Clustering, Distances, JLD

@everywhere function f(Sigma_f,c,wf,r)
    1/2*vecdot(wf, Sigma_f*wf) - r*vecdot(c,wf)
end

@everywhere function f_pos(Sigma_f,c,wf,r)
    1/2*vecdot(wf, Sigma_f*wf) + r*vecdot(c,wf)
end

@everywhere function GenerateData(n_f,Sigma_f,c,Af,b,T,ra)
   x = SharedArray{Float64,2}((n_f,T));
   @parallel (+) for i = 1:T
       genData = Model( solver=GurobiSolver( OutputFlag=0, MIPGap = 1e-6, TimeLimit = 600 ) )
       @variable(genData, wf[1:n_f] )
       @objective(genData, Min, f(Sigma_f,c,wf,ra) )
       @constraint(genData, Af*wf .>= b )
       solve(genData)
       x[:,i] = getvalue(wf)
   end
   return x
end

@everywhere function GenerateData_QEP(n_f,Sigma_f,c,Af,b,T,ra)
   x = SharedArray{Float64,2}((n_f,T));
   @parallel (+) for i = 1:T
       genData = Model( solver=GurobiSolver(OutputFlag=0, MIPGap = 1e-6, TimeLimit = 300 ) )
       @variable(genData, wf[1:n_f] )
       @objective(genData, Min, f_pos(Sigma_f,c,wf,ra) )
       @constraint(genData, Af*wf .== b )
       solve(genData)
       x[:,i] = getvalue(wf)
   end
   return x
end

@everywhere function compute_c(A_return)
#    A_return: matrix of size (num of time steps, num of assets)

    num_asset = size(A_return)[2];
    c = SharedArray{Float64,1}(num_asset);
    for ind_asset = 1:num_asset
        A_vec = A_return[:, ind_asset];
        nonzeroInd = find(x->x!=0, A_vec);
        non_zeroA_vec = A_vec[nonzeroInd];
        c[ind_asset] = mean(non_zeroA_vec);
    end
   return c
end

@everywhere function compute_c_max(A_return)
#    A_return: matrix of size (num of time steps, num of assets)

    num_asset = size(A_return)[2];
    c = SharedArray{Float64,1}(num_asset);
    for ind_asset = 1:num_asset
        A_vec = A_return[:, ind_asset];
        nonzeroInd = find(x->x>=0, A_vec);
        non_zeroA_vec = A_vec[nonzeroInd];
        c[ind_asset] = median(non_zeroA_vec);
    end
   return c
end

@everywhere function remove_extreme_negative_return(A_return)
    for i=1:size(A_return)[1]
        for j=1:size(A_return)[2]
           if A_return[i,j]<-10
#                 println(string(i));
#                 println(string(j));
#                 println(A_return[i,j]);
                A_return[i,j] = 0 ;
            end
        end
    end
    A_return[isnan.(A_return)] .= 0;
    return A_return
end

@everywhere function remove_allzeros_rows(A_return)

    if sum(broadcast(abs, A_return[:,1]))>0.0
        R = A_return[:,1]
        for i=2:size(A_return)[2]
            if sum(broadcast(abs, A_return[:,i]))>0.0
                R = hcat(R, A_return[:,i]);
            end
        end
    else
        R = A_return[:,2]
        for i=3:size(A_return)[2]
            if sum(broadcast(abs, A_return[:,i]))>0.0
                R = hcat(R, A_return[:,i]);
            end
        end
    end
    return R
end

@everywhere function remove_allzeros_rows_withobs(A_return, X_obs)
    if (size(A_return,2)!=size(X_obs,2))
        println("Unequal sample size");
        return None, None
    else
        if sum(broadcast(abs, A_return[:,1]))>0.0
            R = A_return[:,1];
            Xr = X_obs[:,1];
            for i=2:size(A_return)[2]
                if sum(broadcast(abs, A_return[:,i]))>0.0
                    R = hcat(R, A_return[:,i]);
                    Xr =hcat(Xr, X_obs[:,i]);
                end
            end
        else
            R = A_return[:,2];
            Xr = X_obs[:,2];
            for i=3:size(A_return)[2]
                if sum(broadcast(abs, A_return[:,i]))>0.0
                    R = hcat(R, A_return[:,i]);
                    Xr =hcat(Xr, X_obs[:,i]);
                end
            end
        end
        return R, Xr
    end

end

@everywhere function decompose(A_return, num_comp)
    # spectral decompositon of the covariance matrix
    cov_mat = cov(A_return);
    Q = cov_mat;

    factor = eigfact(Q);
    eigvalues = factor[:values];
    eigvectors = factor[:vectors];

    large_eigs = eigvalues[length(eigvalues)-num_comp+1:length(eigvalues)];

    F = eigvectors[:,length(eigvalues)-num_comp+1:length(eigvalues)];

    Sigma_f = Diagonal(large_eigs);

   return F, Sigma_f
end

@everywhere function decompose_new(A_return, num_comp)
    # spectral decompositon of the covariance matrix
    cov_mat = cov(A_return);
    Q = cov_mat;
    Q[isnan(Q)] = 0;

    da = diag(Q);
    da[da.==0]=0;
    n=size(Q,1);
    Q[1:n+1:end]=da;

    factor = eigfact(Q);
    eigvalues = factor[:values];
    eigvectors = factor[:vectors];

    large_eigs = eigvalues[length(eigvalues)-num_comp+1:length(eigvalues)];

    F = eigvectors[:,length(eigvalues)-num_comp+1:length(eigvalues)];

    Sigma_f = Diagonal(large_eigs);

   return F, Sigma_f
end

@everywhere function qp_r_online_ver2(m,nf,M,eta,Sigma_f,c,Af,b,wf,ra)
   problem = Model(solver=GurobiSolver(OutputFlag = 0, MIPGap = 1e-6, TimeLimit = 300) )
   @variable(problem, 0.000001<= ra_t<= 100) # risk aversion value at time t
   @variable(problem, x[1:nf])
#    @variable(problem, 0<=x[1:n]<=1 )
   @variable(problem, u[1:m] >= 0 )
   @variable(problem, z[1:m], Bin )
   @objective(problem, Min, vecdot(ra - ra_t,ra - ra_t)/2 + eta*vecdot(wf - x,wf - x) )
   @constraint(problem, Af*x .>= b )
   @constraint(problem, u .<= M*z )
   @constraint(problem, Af*x - b .<= M*(1-z))
   @constraint(problem, Sigma_f*x - ra_t*c - Af'*u .== zeros(nf) )
   solve(problem)
   return getvalue(ra_t), getvalue(x)
end

@everywhere function qp_r_online_ver3(m,nf,M,eta,Sigma_f,c,Af,b,wf,ra)
   problem = Model(solver=GurobiSolver(OutputFlag = 0, MIPGap = 1e-6, TimeLimit = 3000) )
   @variable(problem, 0.000001<= ra_t<= 100) # risk aversion value at time t
   @variable(problem, x[1:nf])
#  @variable(problem, 0<=x[1:n]<=1 )
   @variable(problem, u[1:m] >= 0 )
   @variable(problem, z[1:m], Bin )
   #@variable(problem, e<=0.0001)
   @objective(problem, Min, vecdot(ra - ra_t,ra - ra_t)/2 + eta*vecdot(wf - x,wf - x))
   @constraint(problem, Af*x .>= b)
   @constraint(problem, u .<= M*z )
   @constraint(problem, Af*x - b .<= M*(1-z))
   @constraint(problem, Sigma_f*x - ra_t*c - Af'*u .== zeros(nf) )
   solve(problem)
   return getvalue(ra_t), getvalue(x)
end

@everywhere function fixdiagonal(sig)
    d = diag(sig);
        if d[length(d)]>5*d[length(d)-1]
            println("eigen value is too large");
            d[length(d)] =  d[length(d)-1]+0.02;
        end
    return Diagonal(d)

end

@everywhere function normalize_doublesigns(vec)
    outvec = deepcopy(vec);
    pos_ind=find(x->x>0, outvec);
    neg_ind =find(x->x<0, outvec);
    sum_neg=sum(outvec[neg_ind]);
    sum_pos=sum(outvec[pos_ind]);

    outvec[pos_ind]=outvec[pos_ind]./sum_pos;
    outvec[neg_ind]=-outvec[neg_ind]./sum_neg;
    return outvec;
end

@everywhere function remove_allzeros_rows_withobs(A_return, X_obs)

    if (size(A_return,2)!=size(X_obs,2))
        println("Unequal sample size");
        return None, None
    else
        if sum(broadcast(abs, A_return[:,1]))>0.0
            R = A_return[:,1];
            Xr = X_obs[:,1];
            for i=2:size(A_return)[2]
                if sum(broadcast(abs, A_return[:,i]))>0.0
                    R = hcat(R, A_return[:,i]);
                    Xr =hcat(Xr, X_obs[:,i]);
                end
            end
        else
            R = A_return[:,2];
            Xr = X_obs[:,2];
            for i=3:size(A_return)[2]
                if sum(broadcast(abs, A_return[:,i]))>0.0
                    R = hcat(R, A_return[:,i]);
                    Xr =hcat(Xr, X_obs[:,i]);
                end
            end
        end
        return R, Xr
    end
end

@everywhere function validation(all_configs, Sigma_f, caf, k, n, F)

    T = 200;
    return_x_err = SharedArray{Float64,2}((size(all_configs,1),T));
    return_r_err = SharedArray{Float64,2}((size(all_configs,1),T));

    for i = 1:size(all_configs,1)

        # obtaint the parameters from the all_configs vector
        M=all_configs[i][1];
        lambda=all_configs[i][2];
        err=all_configs[i][3];
        ra_guess=all_configs[i][4];
        ra_true=all_configs[i][5];

        Af = -[eye(n);ones(1,n);-ones(1,n);-eye(n)]*F;

        bf = -[ones(n,1);1;-1;err*ones(n,1)];
        mf = length(bf);

        x_err_vec = SharedArray{Float64}(T);
        r_err_vec = SharedArray{Float64}(T);

        for j = 1:T

            x_recon = GenerateData(k,Sigma_f,caf,Af,bf,1,ra_true);
            ra_guess, x2 = qp_r_online_ver3(mf,k,M,lambda*j^(-1/2),Sigma_f,caf,Af,bf,x_recon,ra_guess);

            x_err_vec[j]=norm(x2-x_recon,2)/norm(x_recon,2)
            r_err_vec[j]=(ra_guess-ra_true).^2./(ra_true.^2)

        end

        return_x_err[i,:]=x_err_vec;
        return_r_err[i,:]=r_err_vec;

    end

    return return_x_err, return_r_err;
end

function plot_asset_eigenspace(outfile)
    c=ColorGradient([:red,:yellow,:blue]);
    plot(ef_riskvec_ineq,ef_profitvec_ineq, linewidth = 4, lc=c, line_z=log(steps), label="inequality constraint directly in eigenspace", color =:red, xtickfontsize=5, ytickfontsize=5, size = (500, 500), legendfontsize = 8, titlefontsize=7,legend=:bottomright, title="inequality constraiant directly in eigenspace", fmt = :png);
    plot!(ef_riskvec_asset_projected,ef_profitvec_asset_projected, seriestype = :scatter, markersize=1, legend = false, color = [:black], size = (500, 500), fmt = :png);
    plot!(ef_riskvec_fund_projected,ef_profitvec_fund_projected, seriestype = :scatter, markersize=5, marker=:cross, legend = false, color = [:cyan], size = (500, 500), fmt = :png);
    png(outfile);
end

function plot_asset_originalspace(outfile)
    c=ColorGradient([:red,:yellow,:blue]);
    plot(ef_riskvec_recon_ineq, ef_profitvec_recon_ineq, linewidth = 4, lc=c, line_z=log(steps), label="inequality constraint mapped to original space", color =:red, xtickfontsize=5, ytickfontsize=5, size = (500, 500), legendfontsize = 8, titlefontsize=7, legend=:bottomright, title="inequality constraint mapped to original space", fmt = :png);
    plot!(ef_riskvec_asset_original,ef_profitvec_asset_original, seriestype = :scatter, markersize=1, legend = false, color = [:black], size = (500, 500), fmt = :png);
    plot!(ef_riskvec_fund_original,ef_profitvec_fund_original, seriestype = :scatter, markersize=5, marker=:cross, legend = false, color = [:cyan], size = (500, 500), fmt = :png);
    png(outfile);
end

function plot_asset_ori(outfile)
    c=ColorGradient([:red,:yellow,:blue]);
    plot(ef_riskvec_ori, ef_profitvec_ori, linewidth = 4, lc=c, line_z=log(steps), label="inequality constraint mapped to original space \n double unnormalized", color = :red, xtickfontsize=5, ytickfontsize=5, size = (500, 500), legendfontsize = 8, titlefontsize=7, legend=:bottomright, title="inequality constraint original", fmt = :png);
    plot!(ef_riskvec_asset_original,ef_profitvec_asset_original, seriestype = :scatter, markersize=1, legend = false, color = [:black], size = (500, 500), fmt = :png);
    plot!(ef_riskvec_fund_original,ef_profitvec_fund_original, seriestype = :scatter, markersize=5, marker=:cross, legend = false, color = [:cyan], size = (500, 500), fmt = :png);
    png(outfile);
end

function plot_asset_sector(outfile)
    c=ColorGradient([:red,:yellow,:blue]);
    plot(ef_riskvec_sector, ef_profitvec_sector, linewidth = 4, lc=c, line_z=log(steps), label="inequality constraint sector space", color = :red, xtickfontsize=5, ytickfontsize=5, size = (500, 500), legendfontsize = 8, titlefontsize=7, legend=:bottomright, title="inequality constraint sector space", fmt = :png);
    plot!(ef_riskvec_fund_sector,ef_profitvec_fund_sector, seriestype = :scatter, markersize=1, legend = false, color = [:black], size = (500, 500), fmt = :png);
    png(outfile);
end

function plot_asset_sector_withRandomPortfolio(outfile)
    x = range(1,1,10);
    labels = string.(collect(x));
    c=ColorGradient([:red,:yellow,:blue]);

    # varying risk aversion value to generate portfolio
    p1=plot(ef_riskvec_sector, ef_profitvec_sector, linewidth=4, lc=c, line_z=log(steps), label="inequality constraint sector space", color = :red, xtickfontsize=10, ytickfontsize=10, size=(800, 800), legendfontsize=8, titlefontsize=7, legend=:bottomright, title="Efficient Frointer in Sector Space", fmt = :png);

    #plot random portfolio <risk,profit> points
    #plot!(ef_riskvec_random, ef_profitvec_random, text=labels, seriestype = :scatter, label="random portfolio", fmt = :png);
    plot!(ef_riskvec_random, ef_profitvec_random, seriestype = :scatter, legend=:bottomright, label="random portfolio", fmt = :png);

    #plot real fund portfolio on sector level
    plot!(ef_riskvec_fund_sector, ef_profitvec_fund_sector, seriestype = :scatter, legend=:bottomright, label="VFINX portfolio",  fmt = :png);

    # guess from random portfolios
    p2=bar(all_r,color="red",alpha=0.4, titlefontsize=7, title="Estimated Risk Aversion \n from Random Portfolios");

    # guess from real fund sector portfolios
    p3=bar(all_fund_r,color="green",alpha=0.4, titlefontsize=7, title="Estimated Risk Aversion \n from Fund Portfolios");

    l = @layout [a{0.8w} [b{0.5h}; b{0.5h}]];

    display(plot(p1, p2, p3, layout = l, legend = false));
    png(outfile);
end

# generate portfolio directly in eigen space
function plot_diagonal_rguess(outfile)
    #i_rand is the
    i_rand = sortperm(r_rand);
    i_benchmark = sortperm(all_benchmark_r);
    plot(i_rand,i_benchmark, seriestype = :scatter, xlabel="True Risk Aversion Order", ylabel="Estimated Risk Aversion Order", xtickfontsize=20, ytickfontsize=20,  titlefontsize=30, title="Order Rank", markersize=r_rand[i_rand]*10, legend = false, color = [:blue], size = (500, 500), fmt = :png);
    plot!([0,30],[0,30], color=[:black], linewidth=0.8);
    png(outfile);
end

# generate portfolio directly in eigen space
function plot_diagonal_rguess_eigen(outfile)
    #i_rand is the
    i_rand = sortperm(r_rand);
    i_benchmark = sortperm(all_benchmark_r);
    plot(i_rand,i_benchmark, seriestype = :scatter, xlabel="True Risk Aversion Order", ylabel="Estimated Risk Aversion Order", xtickfontsize=10, ytickfontsize=10,  titlefontsize=10, title="Order Rank", markersize=r_rand[i_rand]/10, legend = false, color = [:orange], size = (500, 500), fmt = :png);
    plot!([0,30],[0,30], color=[:black], linewidth=0.8);
    png(outfile);
end

# generate portfolio directly in eigen space
function plot_diagonal_rguess_truevalue(outfile)
    plot(r_rand,all_benchmark_r, seriestype = :scatter, xlabel="True Risk Aversion Value", ylabel="Estimated Risk Aversion Value", xtickfontsize=20, ytickfontsize=20,  titlefontsize=30, title="Comparison Rank", markersize=r_rand/2, legend = false, color = [:blue], size = (500, 500), fmt = :png);

    plot!([0,80],[0,80], color=[:black], linewidth=0.8);
    png(outfile);
end

# generate portfolio in the original space, and project it to eigenspace
function plot_diagonal_rguess_projection(outfile)
    i_rand_p = sortperm(r_rand_ori);
    i_benchmark_p = sortperm(all_benchmark_projected_r);
    plot(i_rand_p,i_benchmark_p, seriestype = :scatter, markersize=r_rand_ori[i_rand_p],  xlabel="True Risk Aversion Order", ylabel="Estimated Risk Aversion Order", legend = false, color = [:orange], size = (500, 500), fmt = :png);
    png(outfile);
end

# generate portfolio in the original space, and project it to eigenspace
function plot_diagonal_rguess_projection_truevalue(outfile)
    plot(r_rand, all_benchmark_r, seriestype = :scatter, title="Comparison Rank", markersize=r_rand/10, xtickfontsize=10, ytickfontsize=10,  titlefontsize=10, xlabel="True Risk Aversion Value", ylabel="Estimated Risk Aversion Value", legend = false, color = [:orange], size = (500, 500), fmt = :png);
    plot!([0,80],[0,80], color=[:black], linewidth=0.8);
    png(outfile);
end

# plot the EFs
function plot_rand_portfolio_onef(outfile)
    c=ColorGradient([:red,:yellow,:blue]);
    p1=plot(ef_riskvec_sector, ef_profitvec_sector, linewidth = 4, lc=c, line_z=log(steps), label="efficient frontier", color = :red, xtickfontsize=15, ytickfontsize=15, size = (500, 500), legendfontsize = 8, titlefontsize=30, legend=:bottomright, title="Efficient frontier", fmt = :png);
    plot!(ef_riskvec_random_onef, ef_profitvec_random_onef, lc=c, markersize=r_rand*15, seriestype = :scatter, label="random portfolio", fmt = :png);
    png(outfile);
end

function plot_asset_eigenspace_mono(outfile)
    c=ColorGradient([:red,:yellow,:blue]);
    # plot EF directly generated in eigenspace
    # color is red
    plot(ef_riskvec_ineq,ef_profitvec_ineq, linewidth = 4, lc=c, line_z=log(steps), xlabel="Risk", ylabel="Profit", color =:red, titlefontsize=10, xtickfontsize=10, ytickfontsize=10, size = (500, 500), legendfontsize = 8, legend=:bottomright, title="Efficient frontier for factor-based portfolio (k=5)", fmt = :png);

    plot!(ef_riskvec_random_onef, ef_profitvec_random_onef, seriestype=:scatter, markersize=r_rand./5, marker=:circle, legend=:bottomright,  label="random portfolio generated in factor space", color = [:orange], size = (500, 500), fmt=:png);

    png(outfile);
end

function plot_asset_eigenspace_k(outfile)
    c=ColorGradient([:red,:yellow,:blue]);
    plot(ef_riskvec_ineq_5,ef_profitvec_ineq_5, linewidth = 1, lc=c, line_z=log(steps), xlabel="Risk", ylabel="Profit", label="k=5", color =:red, titlefontsize=10, xtickfontsize=10, ytickfontsize=10, size = (500, 500), legendfontsize = 8, legend=:bottomright, title="inequality constriant directly in eigenspace", fmt = :png);
    plot!(ef_riskvec_ineq_10,ef_profitvec_ineq_10, linewidth = 2, lc=c, line_z=log(steps), xlabel="Risk", ylabel="Profit", label="k=10", color =:red, titlefontsize=10, xtickfontsize=10, ytickfontsize=10, size = (500, 500), legendfontsize = 8, legend=:bottomright, fmt = :png);
    plot!(ef_riskvec_ineq_15,ef_profitvec_ineq_15, linewidth = 2, lc=c, line_z=log(steps), xlabel="Risk", ylabel="Profit", label="k=15", color =:red, titlefontsize=10, xtickfontsize=10, ytickfontsize=10, size = (500, 500), legendfontsize = 8, legend=:bottomright, fmt = :png);
    plot!(ef_riskvec_ineq_20,ef_profitvec_ineq_20, linewidth = 3, lc=c, line_z=log(steps), xlabel="Risk", ylabel="Profit", label="k=20", color =:red, titlefontsize=10, xtickfontsize=10, ytickfontsize=10, size = (500, 500), legendfontsize = 8, legend=:bottomright, fmt = :png);
    plot!(ef_riskvec_ineq_25,ef_profitvec_ineq_25, linewidth = 4, lc=c, line_z=log(steps), xlabel="Risk", ylabel="Profit", label="k=25", color =:red, titlefontsize=10, xtickfontsize=10, ytickfontsize=10, size = (500, 500), legendfontsize = 8, legend=:bottomright, fmt = :png);
    plot!(ef_riskvec_ineq_30,ef_profitvec_ineq_30, linewidth = 5, lc=c, line_z=log(steps), xlabel="Risk", ylabel="Profit", label="k=30", color =:red, titlefontsize=10, xtickfontsize=10, ytickfontsize=10, size = (500, 500), legendfontsize = 8, legend=:bottomright, fmt = :png);
    png(outfile);
end

max_stars_repo_name: cyberyu/ipo

max_issues_repo_name: cyberyu/ipo

max_forks_repo_name: cyberyu/ipo