text:
function methods!(KS::KineticSolver, psi::Array{Float64,3}, psiOld::Array{Float64,3}, k::Int64)

    if (KS.ThisQuadratureType == "octa" || KS.ThisQuadratureType == "ico" ) && KS.rotationmagnitude > 0 # rSN rotation step
        u = randn(3);
        u = u / norm(u);
        theta = KS.rotationmagnitude*KS.dt / KS.norder;
        rotateandinterpolate!(u,theta,KS.Q,psiOld,psi,"forward") # writes into psiOld
        @inbounds for i=1:length(psi)
            psi[i] = psiOld[i]
        end

    elseif KS.ThisQuadratureType == "tensorized" && KS.rotationmagnitude > 0 # Rotation just around the z-axis
        if k % 2 == 0
            theta = KS.rotationmagnitude*KS.dt / KS.norder;       # and ensure it lies in Delta psi
            alpha = theta / (pi / KS.norder) % 1.0
            rotateOnly!([0.0; 0.0; 1.0],theta,KS.Q,"forward")
            psiInterp = zeros(Int(KS.norder / 2),2 * KS.norder);
            for i = 1:size(psi,2)
                for j = 1:size(psi,3)
                    psiMat = reshape(psi[:, i, j],Int(KS.norder / 2),2 * KS.norder);
                    for q = 1:(2 * KS.norder - 1)
                        psiInterp[:, q] = (1 - alpha) * psiMat[:, q] + alpha*psiMat[:,q + 1];
                    end
                    psiInterp[:, 2 * KS.norder] = (1 - alpha) * psiMat[:, 2 * KS.norder] + alpha*psiMat[:,1];
                    psi[:, i, j] = psiInterp[:];
                end
            end
        else
            theta = KS.rotationmagnitude*KS.dt / KS.norder;       # and ensure it lies in Delta psi
            alpha = theta % (pi / KS.norder)
            rotateOnly!([0.0; 0.0; -1.0],theta,KS.Q,"forward")
            psiInterp = zeros(Int(KS.norder / 2),2 * KS.norder);
            for i = 1:size(psi,2)
                for j = 1:size(psi,3)
                    psiMat = reshape(psi[:, i, j],Int(KS.norder / 2),2 * KS.norder);
                    for q = 2:2 * KS.norder
                        psiInterp[:, q] = (1 - alpha) * psiMat[:, q] + alpha*psiMat[:,q - 1];
                    end
                    psiInterp[:, 1] = (1 - alpha) * psiMat[:, 1] + alpha*psiMat[:,2 * KS.norder];
                    psi[:, i, j] = psiInterp[:];
                end
            end
        end
    end

    if KS.convolutionmagnitude > 0 # Idea of Cory Hauck, perform convolution in the angular domain

        @inbounds @views        for j = 1:size(psi,3)
            for i = 1:size(psi,2)
                if norm(psi[:,i,j])>1e-12
                    mul!(KS.tmp, KS.ScatteringKernelConvolution, KS.Q.weights.*psi[:,i,j])
                    psi[:,i,j] .+= KS.dt .* KS.convolutionmagnitude.*(KS.tmp - KS.OutscatteringVector.*KS.Q.weights.*psi[:,i,j]);
                end
            end
        end
    end
end

function solveFluxUpwind!(KS::KineticSolver, psi::Array{Float64,3}, flux::Array{Float64,3})
    # computes the numerical flux over cell boundaries for each ordinate
    # for faster computation, we split the iteration over quadrature points
    # into four different blocks: North West, Nort East, Sout West, South East
    # this corresponds to the direction the ordinates point to
    idxPosPos = findall((KS.Q.pointsxyz[:,1].>=0.0) .&(KS.Q.pointsxyz[:,2].>=0.0))
    idxPosNeg = findall((KS.Q.pointsxyz[:,1].>=0.0) .&(KS.Q.pointsxyz[:,2].<0.0))
    idxNegPos = findall((KS.Q.pointsxyz[:,1].<0.0)  .&(KS.Q.pointsxyz[:,2].>=0.0))
    idxNegNeg = findall((KS.Q.pointsxyz[:,1].<0.0)  .&(KS.Q.pointsxyz[:,2].<0.0))

    # PosPos
    @inbounds for j=KS.rangex,i=KS.rangey, q = idxPosPos
        s0 = psi[q,i,j-3]
        s1 = psi[q,i,j-2]
        s2 = psi[q,i,j-1]
        s3 = psi[q,i,j]
        northflux = 7/4 * s3 - s2 + 1/4 * s1
        southflux = 7/4 * s2 - s1 + 1/4 * s0

        s0 = psi[q,i-3,j]
        s1 = psi[q,i-2,j]
        s2 = psi[q,i-1,j]
        s3 = psi[q,i,j]
        eastflux = 7/4 * s3 - s2 + 1/4 * s1
        westflux = 7/4 * s2 - s1 + 1/4 * s0

        flux[q,i,j] = KS.Q.pointsxyz[q,1] ./KS.dx .* (eastflux-westflux) +
        KS.Q.pointsxyz[q,2]./KS.dy .* (northflux-southflux)
    end
    #PosNeg
    @inbounds for j=KS.rangex,i=KS.rangey,q = idxPosNeg
        s0 = psi[q,i,j+3]
        s1 = psi[q,i,j+2]
        s2 = psi[q,i,j+1]
        s3 = psi[q,i,j]
        northflux = 7/4 * s2 - s1 + 1/4 * s0
        southflux = 7/4 * s3 - s2 + 1/4 * s1

        s0 = psi[q,i-3,j]
        s1 = psi[q,i-2,j]
        s2 = psi[q,i-1,j]
        s3 = psi[q,i,j]
        eastflux = 7/4 * s3 - s2 + 1/4 * s1
        westflux = 7/4 * s2 - s1 + 1/4 * s0

        flux[q,i,j] = KS.Q.pointsxyz[q,1] ./KS.dx .*(eastflux-westflux) +
        KS.Q.pointsxyz[q,2] ./KS.dy .*(northflux-southflux)
    end

    # NegPos
    @inbounds for j=KS.rangex,i=KS.rangey,q = idxNegPos
        s0 = psi[q,i,j-3]
        s1 = psi[q,i,j-2]
        s2 = psi[q,i,j-1]
        s3 = psi[q,i,j]
        northflux = 7/4 * s3 - s2 + 1/4 * s1
        southflux = 7/4 * s2 - s1 + 1/4 * s0

        s0 = psi[q,i+3,j]
        s1 = psi[q,i+2,j]
        s2 = psi[q,i+1,j]
        s3 = psi[q,i,j]
        westflux = 7/4 * s3 - s2 + 1/4 * s1
        eastflux = 7/4 * s2 - s1 + 1/4 * s0

        flux[q,i,j] = KS.Q.pointsxyz[q,1]./KS.dx .*(eastflux-westflux) +
        KS.Q.pointsxyz[q,2] ./KS.dy .*(northflux-southflux)
    end

    # NegNeg
    @inbounds for j=KS.rangex,i=KS.rangey,q = idxNegNeg
        s0 = psi[q,i,j+3]
        s1 = psi[q,i,j+2]
        s2 = psi[q,i,j+1]
        s3 = psi[q,i,j]
        northflux = 7/4 * s2 - s1 + 1/4 * s0
        southflux = 7/4 * s3 - s2 + 1/4 * s1

        s0 = psi[q,i+3,j]
        s1 = psi[q,i+2,j]
        s2 = psi[q,i+1,j]
        s3 = psi[q,i,j]
        westflux = 7/4 * s3 - s2 + 1/4 * s1
        eastflux = 7/4 * s2 - s1 + 1/4 * s0

        flux[q,i,j] = KS.Q.pointsxyz[q,1] ./KS.dx .*(eastflux-westflux) +
        KS.Q.pointsxyz[q,2] ./KS.dy .*(northflux-southflux)
    end
end

function Euler!(KS::KineticSolver, psi::Array{Float64,3}, flux::Array{Float64,3}, dt::Float64)
    # computes the Euler step
    @inbounds for j=KS.rangex,i=KS.rangey
        integral = 0.0;
        @inbounds for q=KS.rangequad
            integral += psi[q, i, j] .*KS.Q.weights[q]
        end
        integral = integral*KS.SigmaS[i,j] / 4.0 / pi;
        sigmaT = KS.SigmaT[i, j]
        @inbounds for q=KS.rangequad
            psi[q, i, j]  += dt .*(-psi[q, i, j] .*sigmaT - flux[q, i, j] + integral)
        end
    end
end

# the first minmod code is the fast version of the second minmod  below that is commented
@inline minmod(x::Float64, y::Float64) = ifelse(x < 0, clamp(y, x, 0.0), clamp(y, 0.0, x))
#@inline function minmod(x::Float64, y::Float64)
#    return sign(x) * max(0.0, min(abs(x),y*sign(x) ) );
#end

@inline function slopefit(left::Float64, center::Float64, right::Float64)
    tmp = minmod(0.5 * (right - left),2.0 * (center - left));
    return minmod(2.0 * (right - center),tmp);
end

function computeMoments(KS::KineticSolver, psi::Array{Float64,3})
    rho = zeros(KS.ny,KS.nx)
    for q = KS.rangequad, j=1:KS.nx, i = 1:KS.ny#KS.rangex, i=KS.rangey
        rho[i, j] += psi[q, i, j] * KS.Q.weights[q]
    end
    return rho
end

function forcePeriodicity!(KS::KineticSolver, psi::Array{Float64,3})
    if KS.periodicXflag
        for q=1:KS.nquadpoints
            for i=1:KS.ny
                psi[q, i, 1] = psi[q,i,end - 3]
                psi[q, i, 2] = psi[q,i,end - 2]
                psi[q,i,end - 1] = psi[q, i, 3]
                psi[q,i,end] = psi[q, i, 4]
            end
        end
    end
    if KS.periodicYflag
        for q=1:KS.nquadpoints
            for j=1:KS.nx
                psi[q, 1, j] = psi[q,end - 3,j]
                psi[q, 2, j] = psi[q,end - 2,j]
                psi[q,end - 1,j] = psi[q, 3, j]
                psi[q,end,j] = psi[q, 4, j]
            end
        end
    end
end

function copyFromToNoGhosts!(psiFrom::Array{Float64,3}, psiTo::Array{Float64,2}, KS::KineticSolver)
    for q=KS.rangequad
        k = 0;
        for i =3:KS.ny - 2
            for j = 3:KS.nx - 3
                k = k + 1
                psiTo[q, k] = psiFrom[q, i, j]
            end
        end
    end
end

function copyFromToWithGhosts!(psiFrom::Array{Float64,2}, psiTo::Array{Float64,3}, KS::KineticSolver)
    for q=KS.rangequad
        k = 0;
        for i =1:KS.ny - 4
            for j = 1:KS.nx - 4
                k = k + 1
                psiTo[q, i + 2, j + 2] = psiFrom[q, k]
            end
        end
    end
end

function solveFlux!(KS::KineticSolver, phi::Array{Float64,3}, flux::Array{Float64,3})
    # computes the numerical flux over cell boundaries for each ordinate
    # for faster computation, we split the iteration over quadrature points
    # into four different blocks: North West, Nort East, Sout West, South East
    # this corresponds to the direction the ordinates point to
    idxPosPos = findall((KS.Q.pointsxyz[:,1].>=0.0) .&(KS.Q.pointsxyz[:,2].>=0.0))
    idxPosNeg = findall((KS.Q.pointsxyz[:,1].>=0.0) .&(KS.Q.pointsxyz[:,2].<0.0))
    idxNegPos = findall((KS.Q.pointsxyz[:,1].<0.0)  .&(KS.Q.pointsxyz[:,2].>=0.0))
    idxNegNeg = findall((KS.Q.pointsxyz[:,1].<0.0)  .&(KS.Q.pointsxyz[:,2].<0.0))

    # PosPos
    @inbounds for j=KS.rangex,i=KS.rangey, q = idxPosPos
        s1 = phi[q,i,j-2]
        s2 = phi[q,i,j-1]
        s3 = phi[q,i,j]
        s4 = phi[q,i,j+1]
        northflux = s3+0.5 .*slopefit(s2,s3,s4)
        southflux = s2+0.5 .*slopefit(s1,s2,s3)

        s1 = phi[q,i-2,j]
        s2 = phi[q,i-1,j]
        s3 = phi[q,i,j]
        s4 = phi[q,i+1,j]
        eastflux = s3+0.5 .*slopefit(s2,s3,s4)
        westflux = s2+0.5 .*slopefit(s1,s2,s3)

        flux[q,i,j] = KS.Q.pointsxyz[q,1] ./KS.dx .* (eastflux-westflux) +
        KS.Q.pointsxyz[q,2]./KS.dy .* (northflux-southflux)
    end
    #PosNeg
    @inbounds for j=KS.rangex,i=KS.rangey,q = idxPosNeg
        s1 = phi[q,i,j-1]
        s2 = phi[q,i,j]
        s3 = phi[q,i,j+1]
        s4 = phi[q,i,j+2]
        northflux = s3-0.5 .* slopefit(s2,s3,s4)
        southflux = s2-0.5 .*slopefit(s1,s2,s3)

        s1 = phi[q,i-2,j]
        s2 = phi[q,i-1,j]
        s3 = phi[q,i,j]
        s4 = phi[q,i+1,j]
        eastflux = s3+0.5 .*slopefit(s2,s3,s4)
        westflux = s2+0.5 .*slopefit(s1,s2,s3)

        flux[q,i,j] = KS.Q.pointsxyz[q,1] ./KS.dx .*(eastflux-westflux) +
        KS.Q.pointsxyz[q,2] ./KS.dy .*(northflux-southflux)
    end

    # NegPos
    @inbounds for j=KS.rangex,i=KS.rangey,q = idxNegPos
        s1 = phi[q,i,j-2]
        s2 = phi[q,i,j-1]
        s3 = phi[q,i,j]
        s4 = phi[q,i,j+1]
        northflux = s3+0.5 .*slopefit(s2,s3,s4)
        southflux = s2+0.5 .*slopefit(s1,s2,s3)

        s1 = phi[q,i-1,j]
        s2 = phi[q,i,j]
        s3 = phi[q,i+1,j]
        s4 = phi[q,i+2,j]
        eastflux = s3-0.5 .*slopefit(s2,s3,s4)
        westflux = s2-0.5 .*slopefit(s1,s2,s3)

        flux[q,i,j] = KS.Q.pointsxyz[q,1]./KS.dx .*(eastflux-westflux) +
        KS.Q.pointsxyz[q,2] ./KS.dy .*(northflux-southflux)
    end

    # NegNeg
    @inbounds for j=KS.rangex,i=KS.rangey,q = idxNegNeg
        s1 = phi[q,i,j-1]
        s2 = phi[q,i,j]
        s3 = phi[q,i,j+1]
        s4 = phi[q,i,j+2]
        northflux = s3-0.5 .* slopefit(s2,s3,s4)
        southflux = s2-0.5 .* slopefit(s1,s2,s3)

        s1 = phi[q,i-1,j]
        s2 = phi[q,i,j]
        s3 = phi[q,i+1,j]
        s4 = phi[q,i+2,j]
        eastflux = s3-0.5 .* slopefit(s2,s3,s4)
        westflux = s2-0.5 .* slopefit(s1,s2,s3)

        flux[q,i,j] = KS.Q.pointsxyz[q,1] ./KS.dx .*(eastflux-westflux) +
        KS.Q.pointsxyz[q,2] ./KS.dy .*(northflux-southflux)
    end
end

max_stars_repo_name: camminady/SN

max_issues_repo_name: camminady/SN

max_forks_repo_name: camminady/SN


text:
using Flux
using Flux: train!
using Flux: throttle
using IterTools

actual(x) = 4x .+ 2

x_train = [1 2]
y_train = actual.(x_train)

opt = Descent(0.001)

data = [(x_train, y_train)]

predict = Dense(1, 1)

parameters = params(predict)

loss(x, y) = Flux.Losses.mse(predict(x), y)

#for epoch in 1:2000
train!(loss, parameters, ncycle(data, 100000), opt, cb = throttle(() -> println("treinando"), 0.002))
#end

loss(x_train, y_train)

parameters

print("fool")

ext: jl


text:
using Random
SansSerifFonts = ["Trebuchet MS", "Heiti TC", "微軟正黑體", "Arial Unicode MS", "Droid Fallback Sans", "sans-serif", "Helvetica", "Verdana", "Hei",
"Arial", "Tahoma", "Trebuchet MS", "Microsoft Yahei", "Comic Sans MS", "Impact", "Segoe Script", "STHeiti", "Apple LiGothic", "MingLiU", "Ubuntu",
"Segoe UI", "DejaVu Sans", "DejaVu Sans Mono", "Noto Sans CJK"]
SerifFonts = ["Baskerville", "Times New Roman", "華康儷金黑 Std", "華康儷宋 Std",  "DFLiKingHeiStd-W8", "DFLiSongStd-W5", "DejaVu Serif", "SimSun",
    "Hiragino Mincho Pro", "LiSong Pro", "新細明體", "serif", "Georgia", "STSong", "FangSong", "KaiTi", "STKaiti", "Courier New"]
CandiFonts = union(SansSerifFonts, SerifFonts)
CandiWeights = ["", " Regular", " Normal", " Medium", " Bold", " Light"]
function checkfonts(fonts::AbstractVector)
    fname = tempname()
    r = Bool[]
    open(fname, "w") do f
        redirect_stderr(f) do
            p = position(f)
            for font in fonts
                rendertext("a", 1 + rand(), font=font) # 相同字体相同字号仅warning一次，故首次执行最准
                # flush(f) #https://en.cppreference.com/w/cpp/io/c/fseek The standard C++ file streams guarantee both flushing and unshifting
                seekend(f)
                p2 = position(f)
                push!(r, p2 == p)
                p = p2
            end
        end
    end
    return r
end
checkfonts(f) = checkfonts([f]) |> only
function filterfonts(;fonts=CandiFonts, weights=CandiWeights)
    candi = ["$f$w" for w in weights, f in fonts] |> vec
    candi[checkfonts(candi)]
end
AvailableFonts = filterfonts()
push!(AvailableFonts, "")

Schemes_colorbrewer = filter(s -> occursin("colorbrewer", colorschemes[s].category), collect(keys(colorschemes)))
Schemes_colorbrewer =  filter(s -> (occursin("Accent", String(s))
        || occursin("Dark", String(s))
        || occursin("Paired", String(s))
        || occursin("Pastel", String(s))
        || occursin("Set", String(s))
        || occursin("Spectral", String(s))
        || occursin("BuPu", String(s))
        || occursin("PuOr", String(s))
        || occursin("RdPu", String(s))
        ), Schemes_colorbrewer)
Schemes_seaborn = filter(s -> occursin("seaborn", colorschemes[s].category), collect(keys(colorschemes)))
Schemes_tableau = filter(s -> occursin("tableau", colorschemes[s].category), collect(keys(colorschemes)))
Schemes_cvd = filter(s -> occursin("cvd", colorschemes[s].category), collect(keys(colorschemes)))
Schemes = [Schemes_colorbrewer..., Schemes_seaborn..., Schemes_tableau..., Schemes_cvd...]

function displayschemes()
    for scheme in Schemes
        display(scheme)
        colors = Render.colorschemes[scheme].colors
        display(colors)
    end
end
function randomscheme()
    if rand() < 0.95
        scheme = rand(Schemes)
        c = Render.colorschemes[scheme].colors
        colors = randsubseq(c, rand())
        isempty(colors) &&  (colors = c)
        println("color scheme: ", repr(scheme), ", length: ", length(colors))
    else
        colors = rand((0, 1, 0, 1, 0, 1, (0, 1), rand(), (rand(), rand())))
        @show colors
    end
    (colors...,)
end
function randomwh(sz::Number=800)
    s = sz * sz
    ratio = (9/16 + rand()*7/16)
    ratio > 0.9 && (ratio = 1.0)
    h = round(Int, sqrt(s * ratio))
    w = round(Int, h / ratio)
    w, h
end
randomwh(sz::Tuple) = sz
randomwh(arg...) = arg
equalwh(sz::Number=800) = sz, sz
equalwh(sz::Tuple) = sz
equalwh(arg...) = arg
function randommask(args...; maskshape=:rand, kargs...)
    rd = Dict(squircle => 0.4, box => 0.6, ellipse => 0.8, ngon => 0.9, star => 1)
    ran = get(rd, maskshape, rand())
    if ran <= 0.4
        return randomsquircle(randomwh(args...)...; kargs...)
    elseif ran <= 0.6
        return randombox(randomwh(args...)...; kargs...)
    elseif ran <= 0.8
        return randomellipse(randomwh(args...)...; kargs...)
    elseif ran <= 0.9
        return randomngon(equalwh(args...)...; kargs...)
    else
        return randomstar(equalwh(args...)...; kargs...)
    end
end
function showcallshape(args...; kargs...)
    ags = [string(args[1]), repr.(args[2:end])..., ("$k=$(repr(v))" for (k, v) in kargs)...]
    println("shape(", join(ags, ", "), ")")
end
function randombox(w, h; cornerradius=:rand, keeparea=false, kargs...)
    if cornerradius == :rand
        r = rand() * 0.5 - 0.05 # up to 0.45
        r < 0. && (r = 0.) # 10% for 0.
        r = round(Int, h * r)
    else
        r = cornerradius
    end
    sc = keeparea ? sqrt(w*h/box_area(w, h, cornerradius=r)) : 1
    w = round(Int, w*sc); h = round(Int, h*sc); r = round(Int, r*sc)
    showcallshape(box, w, h; cornerradius=r, kargs...)
    return shape(box, w, h; cornerradius=r, kargs...)
end
function randomsquircle(w, h; rt=:rand, keeparea=false, kargs...)
    if rt == :rand
        if rand() < 0.8
            rt = rand()
        else
            ran = rand()
            if ran < 0.5
                rt = 2
            else
                rt = 1 + 1.5rand()
            end
        end
        rt = round(rt, digits=3)
    end
    sc = keeparea ? sqrt(w*h/squircle_area(w, h, rt=rt)) : 1
    w = round(Int, w*sc); h = round(Int, h*sc)
    showcallshape(squircle, w, h, rt=rt; kargs...)
    return shape(squircle, w, h, rt=rt; kargs...)
end
function randomellipse(w, h; keeparea=false, kargs...)
    sc = keeparea ? sqrt(w*h/ellipse_area(w, h)) : 1
    w = round(Int, w*sc); h = round(Int, h*sc)
    showcallshape(ellipse, w, h; kargs...)
    return shape(ellipse, w, h; kargs...)
end
function randomorientation(n)
    if n == 3
        ori = rand((0, π/2, π/3))
    elseif n % 2 == 0
        ori = rand((0, π/n))
    else
        ori = 0
    end
    return ori
end
function randomngon(w, h; npoints=:rand, orientation=:rand, keeparea=false, kargs...)
    npoints == :rand && (npoints = rand(3:12))
    orientation == :rand && (orientation = randomorientation(npoints))
    sc = keeparea ? sqrt(w*h/ngon_area(w, h, npoints=npoints)) : 1
    w = round(Int, w*sc); h = round(Int, h*sc)
    showcallshape(ngon, w, h; npoints=npoints, orientation=orientation, kargs...)
    return shape(ngon, w, h; npoints=npoints, orientation=orientation, kargs...)
end
function randomstar(w, h; npoints=:rand, starratio=:rand, orientation=:rand, keeparea=false, kargs...)
    npoints == :rand && (npoints = rand(5:12))
    orientation == :rand && (orientation = randomorientation(npoints))
    if starratio == :rand
        starratio = cos(π/npoints) * (0.7 + 0.25rand())
        starratio = round(starratio, digits=3)
    end
    sc = keeparea ? sqrt(w*h/star_area(w, h, npoints=npoints, starratio=starratio)) : 1
    w = round(Int, w*sc); h = round(Int, h*sc)
    showcallshape(star, w, h; npoints=npoints, starratio=starratio, orientation=orientation, kargs...)
    return shape(star, w, h; npoints=npoints, starratio=starratio, orientation=orientation, kargs...)
end
function randomangles()
    a = rand((-1, 1)) .* rand((0, (0, 90), (0, 90, 45), (0, 90, 45, -45), (0, 45, -45), (45, -45), -90:90, 0:90))
    println("angles = ", a)
    a
end
function randommaskcolor(colors)
    colors = parsecolor.(unique(colors))
    try
        g = Gray.(colors) |> sort
        m = g[1]
        M = g[end]
        if length(g) > 1
            d = diff(g)
            I = maximum(d)
            i = findlast(isequal(I), d)
        else
            I = 0
            i = -1
        end
        # @show I, m, M
        if I > 3(1 - M) && I > 3m
            middle = (g[i] + g[i + 1]) / 2
            th1 = clamp(max(g[i] + 0.15, middle - rand(0:0.001:0.05)), 0, middle)
            th2 = clamp(min(g[i + 1] - 0.15, middle + rand(0:0.001:0.05)), middle, 1)
            default = middle
        elseif sum(g) / length(g) < 0.7 && (m + M) / 2 < 0.7 && !(m > 2(1 - M))# 明亮
            th1 = clamp(max(M + 0.15, rand(0.85:0.001:1.0)), 0, 1)
            th2 = clamp(th1 + 0.1, 0, 1)
            default = 1.0
        else    # 黑暗
            th2 = clamp(min(m - 0.15, rand(0.0:0.001:0.3)), 0, 1) # 对深色不敏感，+0.15
            th1 = clamp(th2 - 0.15, 0, 1)
            default = 0.0
        end
        maskcolor = rand((default,
        (round(rand(th1:0.001:th2), digits=3),
        round(rand(th1:0.001:th2), digits=3),
        round(rand(th1:0.001:th2), digits=3))))
        # @show maskcolor
        return maskcolor
    catch e
        @show e
    @show "colors sum failed", colors
        return "white"
    end
end
function randomlinecolor(colors)
    if rand() < 0.8
        linecolor = rand((colors[1], colors[1], rand(colors)))
    else
        linecolor = (
            round(rand(), digits=3),
            round(rand(), digits=3),
            round(rand(), digits=3),
            min(1., round(0.5 + rand()/2, digits=3)))
    end
    linecolor
end
randomoutline() = rand((0, 0, 0, rand(2:10)))
function randomfont()
    font = rand(AvailableFonts)
    @show font
    font
end


text:
# # 生物資訊基礎教學系列：Julia程式語言基礎
# #### Date: 2020-06-28
# #### Speaker：米夏糯

# ---

# ## 本日介紹內容包含了
# + 學習程式語言的目的及意義
# + Julia程式語言簡介
# + Julia程式語言基礎
# + 學完基礎語法之後

# ---

# # 學習程式語言的目的及意義
# ## 學了程式語言能做些什麼事情
#

# ---

# # 學習程式語言的目的及意義
# ## __沒有!__
#
# >> 寫程式絕對是適合男女老幼的，「寫程式的都是肥宅」只是大家對它的誤解。寫程式更加是一種藝術，一種無中生有的精神。

# ---

# ## 初學程式語言時所需要具備的心態: __Grit__
#

# ---

# # Julia程式語言簡介
# ## Julia程式語言的發展脈絡
# + 想要開源的語言、自由的版權、C語言的速度、如Lisp一般的一致性、支援像Matlab一樣的數學符號、如Python一般的通用、處理統計能像R一樣簡單、處理字串如Perl般的自然，同時能像shell一樣串接程式。
#

# ---

# ## Julia語言的特性 (1/3)
# + Write like Python, run like C
# + 高效的運算速度
# + 直接簡單的語法
# + 同時支援動態、靜態類型
# + 支援多種程式設計風格

# ---

# ## Julia語言的特性 (2/3)
# #### Benchmark
#

# ---

# ## Julia語言的特性 (3/3)
#

# ---

# ## 學習Julia語言的管道途徑以及學習資源 (1/2)
# + 書籍
#   - [Julia 程式設計：新世代資料科學與數值運算語言](https://www.tenlong.com.tw/products/9789577633828?list_name=srh)
#   - [Python 接班人出世：最新科學專用語言 Julia 入門實戰](https://www.tenlong.com.tw/products/9789865501006?list_name=srh)
# + 影片
#   - [Julia Taiwan Meetup](https://www.youtube.com/channel/UCHS7SCM_GDDMreNQfW5LFzQ/videos)
#   - [阿杜的YouTube頻道](https://www.youtube.com/c/%E6%9D%9C%E5%B2%B3%E8%8F%AF/videos)
# + 線上課程
#   - [Julia Academy](https://juliaacademy.com/)
#   - [Udemy](https://www.udemy.com/topic/julia-programming-language/)

# ---

# ## 學習Julia語言的管道途徑以及學習資源 (2/2)
# + [Julia Taiwan社群所翻譯之文件](https://hackmd.io/@7WeiUEuJSBKp7WCRouAWVg/r106dSpkb/%2FKwBhCMBYA4GYE4C0B2e4BsiYEYCGjwAmZaRbddbcEbSEeAExCA%3D%3D?type=book)
#

# ---

# # Julia程式語言基礎
# ## 安裝Julia (1/5)
#

# ---

# ## 安裝Julia (2/5)
#

# ---

# ## 安裝Julia (3/5)
#

# ---

# ## 安裝Julia (4/5)
#

# ---

# ## 安裝Julia (5/5)
#

# ---

# ## 撰寫Julia時使用的工具平台 - VS Code (1/6)
#

# ---

# ## 撰寫Julia時使用的工具平台 - VS Code (2/6)
#

# ---

# ## 撰寫Julia時使用的工具平台 - VS Code (3/6)
#

# ---img/

# ## 撰寫Julia時使用的工具平台 - VS Code (4/6)
#

# ---

# ## 撰寫Julia時使用的工具平台 - VS Code (5/6)
#

# ---

# ## 撰寫Julia時使用的工具平台 - VS Code (6/6)
#

# ---

# ## 撰寫Julia時使用的工具平台 - Jupyter Notebook (1/6)
# + 安裝IJulia
#

# ---

# ## 撰寫Julia時使用的工具平台 - Jupyter Notebook (2/6)
# + 安裝Anaconda (1/2)
#

# ---

# ## 撰寫Julia時使用的工具平台 - Jupyter Notebook (3/6)
# + 安裝Anaconda (2/2)
#

# ---

# ## 撰寫Julia時使用的工具平台 - Jupyter Notebook (4/6)
# + 安裝Jupyter Notebook (1/3)
#

# ---

# ## 撰寫Julia時使用的工具平台 - Jupyter Notebook (5/6)
# + 安裝Jupyter Notebook (2/3)
#

# ---

# ## 撰寫Julia時使用的工具平台 - Jupyter Notebook (6/6)
# + 安裝Jupyter Notebook (3/3)
#

# ---

# ## 基礎語法與概念 - 寫出你的第一支hello world

# + 進入Julia環境
#

# + 輸入以下程式碼，並按Enter送出
println("Hello world")

# ---

# ## 基礎語法與概念 - 註解

# + 單行註解
## 井字號開頭的是註解
println("上面是註解")

# + 多行註解
#=
如果想要使用多行註解,
可以用#=及=#夾住所需要註解的內容
=#

# ---

# ## 基礎語法與概念 - 變數 (1/4)
# > 在Julia裡頭，一個變數是指與一個值所存在的記憶體空間__有關__（或稱為__綁定__）的一個名字。當我們想在後續使用目前運算後所得到的值，透過變數存取就會比直接對記憶體位置存取來的簡單易懂。
#
## 將5這個值指定到一個名為x的變數, 這邊的=所表示的是"指派"的意思而非"等於"
## 這個過程我們也稱之為"宣告"
x = 5
## 輸入變數名稱x, 我們就可以獲得這個變數所存取之值的內容
x

# ---

# ## 基礎語法與概念 - 變數 (2/4)
## 將前一張投影片中所宣告的變數x加上6後所得運算結果重新指派給x這個變數
x = x + 6
## 此時x的值就會等於11
x

# ---

# ## 基礎語法與概念 - 變數 (3/4)
# #### 變數命名規則
# + 開頭僅能使用字母(A-Z, a-z)、底線或是__Unicode__(須大於00A0)的字元
#   - _123 (valid)
#   - Abcd (valid)
#   - 5566 (invalid)
#   - δ (valid)
# + 不能使用保留字作為變數名稱
#   - if
#   - while
#   - local
#   - end
#   - true

# ---

# ## 基礎語法與概念 - 變數 (4/4)
# #### 支援[Unicode](https://docs.julialang.org/en/v1/manual/unicode-input/) & LaTex語法作為變數名稱

# + Example 1:
## 輸入\delta 再按tab鍵，就會產生一個希臘字母 δ
## 宣告一個名為α₂的變數，並綁定到0.23這個值
α₂ = 0.23

# + Example 2:
🌱 = 0
🚰 = 5
🌅 = 5
🌳 =  🌱 + 🚰 + 🌅

# ---

# ## 基礎語法與概念 - 加減乘除
# + 加
2 + 5
# + 減
2 - 5
# + 乘
3 * 8
# + 除
8 / 3

# ---

# ## 基礎語法與概念 - 浮點數
# > 有小數點的數字在程式語言中就被稱為__浮點數__
# + 純整數相加
5 + 5
# + 浮點數相加
5.0 + 5

# ---

# ## 基礎語法與概念 - 型別
# + 使用`typeof()`來檢查型別
typeof(5)

# + 檢查浮點數
typeof(5.0)

# + 在Julia裡頭使用到的型別主要有：
#   - Int (Int8, Int16, Int32, Int64, UInt etc.)
#   - Float (Float16, Float32, Float64 etc.)
#   - Char
#   - String
#   - Boolean

# ---

# ## 基礎語法與概念 - 數字與運算子 (1/4)
# #### 整數 (1/2)
# > 如果沒有特別宣告，會依系統位元數決定
#

# ---

# ## 基礎語法與概念 - 數字與運算子 (1/4)
# #### 整數 (2/2)
# + 宣告非預設型別的整數
Int8(10)
typeof(ans)
# + 取得不同位元數的上下限
typemax(Int8) # 取得上限
typemin(Int8) # 取得下限
# + 溢位
typemax(Int8) + Int8(1)

# ---

# ## 基礎語法與概念 - 數字與運算子 (2/4)
# #### 浮點數
# + 含小數點浮點數
#
# + 特殊浮點數
#

# ---

# ## 基礎語法與概念 - 數字與運算子 (3/4)
# #### 常數及係數
# + 預設常數
pi
# + 宣告常數
const C = 1.23456
# + 字面係數
x = 3
2x^2 - 3x + 1

# ---

# ## 基礎語法與概念 - 數字與運算子 (4/4)
# #### 運算子 (1/3)
# + 算術運算子
## -x: 變號
## x + y, x - y, x * y, x / y: 四則運算
## div(x,y): 求商
## x % y: 求餘
## x \ y: 與 x / y 相同
## x ^ y: 次方

# + 比較運算子
## x == y: 等於
## x != y, x ≠ y: 不等於
## x < y: 小於
## x > y: 大於
## x <= y, x ≤ y: 小於等於
## y >=, x ≥ y: 大於等於

# ---

# ## 基礎語法與概念 - 數字與運算子 (4/4)
# #### 運算子 (2/3)

## 1.更新運算子
## 舊寫法

x = 5
y = 0
y = y + 2x
## 使用更新運算子的寫法
y += 2x
## 其他更新運算子：-=, *=, /=, %=等

## 2. 布林運算子
## x = true; typeof(x)
## !x, ~x: 非 (negation)
## x & y: 且
## x | y: 或

# ---

# ## 基礎語法與概念 - 控制流程 (1/3)

# #### if 條件判斷 (1/2)
# + if條件判斷語法結構：
## if <判斷式>
##    <運算式>
## end

x = 0
y = 5
if x < y
    println("x is less than y")
end

# ---

# ## 基礎語法與概念 - 控制流程 (1/3)

# #### if 條件判斷 (2/2)
# + if-else條件判斷語法結構：
## if <判斷式>
##    <運算式-1>
## else
##    <運算式-2>
## end

# + if-elseif-else條件判斷語法結構：
## if <判斷式-1>
##    <運算式-1>
## elseif <判斷式-2>
##    <運算式-2>
## else
##    <運算式-3>
## end

# ---

# ## 基礎語法與概念 - 控制流程 (2/3)

# #### for 迴圈 (1/2)
# + for迴圈語法結構：
## for i = <起始>:<結束>
##     <運算式程式碼區塊>
## end

a = 0
for i = 1:100
    global a
    a += i
end
a

# ---

# ## 基礎語法與概念 - 控制流程 (2/3)

# #### for 迴圈 (2/2)
# + 巢狀for迴圈，以九九乘法表為例：
for i = 1:9
    for j = 1:9
        print(i, " x ", j, " = ", i*j, "\t")
    end
    println()
end

# ---

# ## 基礎語法與概念 - 控制流程 (3/3)

# #### while 迴圈
# + while 迴圈語法結構：
## while <持續條件>
##     <運算式程式碼區塊>
## end

## Example:
a = 1
while a <= 100
    global a
    a += 2
end
a

# ---

# ## 基礎語法與概念 - 函式

# > 函式就像是一段特地的程式碼區塊，專門用來完成特定的工作，為了簡化同時又不用重複寫這段程式碼，我們就可以將這段程式碼定義成一個__函式__。隨後當需要執行這項工作時，只需要呼叫此函式即可。
# + 定義一個加法
function add(x, y)
    return x + y
end
## 呼叫這個定義好的函式
add(1, 2)
# + 函式也可以這樣定義
add2(x, y) = x + y
add2(1,2)

# ---

# ## 基礎語法與概念 - 集合容器 (1/4)

# #### 陣列 (1/4)
# > 把同性質的資料按順序排在一起，像是一整排的格子那樣的排列，就稱為__陣列__

# + 創建一個空陣列
x = []

# + 也可以把元素型別放在方括號前創建
Any[]

# + 創建一個Int64的空陣列
Int64[]

# ---

# ## 基礎語法與概念 - 集合容器 (1/4)

# #### 陣列 (2/4)
# + Example 1:
x = [1, 2, 3]
# + Example 2:
x = [1, 1.2]
# + Example 3:
Int8[1, 2, 3, 4]

# ---

# ## 基礎語法與概念 - 集合容器 (1/4)

# #### 陣列 (3/4)
# + 索引 (indexing) / 取值
x = [1, 1.2, 2.4]
x[1]
# + 求陣列大小
length(x)
# + 取得某區段的元素
x[2:length(x)]

# ---

# ## 基礎語法與概念 - 集合容器 (1/4)

# #### 陣列 (4/4)
# + 不得使用超出陣列大小的索引值
x = [6.0, 3.2, 7.6, 0.9, 2.3]
# x這個陣列一共5個元素，要是這樣想去取第六個元素如: x[6]，就會出現錯誤
#

# + 其他操作
## push!: 從陣列後面加入新的值， push!(x, 9.0)
## append!: 串接兩個陣列, y = [10.0, 3.4]; append!(x, y)
## pop!: 從陣列最後移出一個元素, pop!(x)

# ---

# ## 基礎語法與概念 - 集合容器 (2/4)

# #### 集合 (Set) (1/4)
# > 如同數學中的集合一樣，不允許重複的元素，每個元素都是唯一的，且元素間也沒有次序
## 把一個陣列放入Set()中來建構一個集合
x = Set([1, 2, 3, 4])
## 增加集合中的元素
push!(x, 5)

# ---

# ## 基礎語法與概念 - 集合容器 (2/4)

# #### 集合 (Set) (2/4)
# + 判斷集合間的關係
y = Set([1, 3, 5])
issubset(y, x) #確認y是否為x的子集, 或輸入 y ⊆ x, 也可得到相同結果

# + 交集
intersect(x, y) # 或使用 x ∩ y

# ---

# ## 基礎語法與概念 - 集合容器 (2/4)

# #### 集合 (Set) (3/4)

# + 聯集
union(x, y) # 或使用 x ∪ y

# 差集
setdiff(x, y)

# ---

# ## 基礎語法與概念 - 集合容器 (2/4)

# #### 集合 (Set) (4/4)

# + 判斷集合是否相等
x == Set([3, 1, 2, 5, 4])

# + 迭代
for i in x
    println(i)
end

# ---

# ## 基礎語法與概念 - 集合容器 (3/4)

# #### 字串 (Strings) (1/4)
# > 一個__字串__可以視為多個__字元 (characters)__的序列組合，但這個組合內容不可變 (immutable)。
# + 字元用__單__引號夾住
typeof('A')
# + 見其ASCII編碼
'A'
# + 字串用__雙__引號夾住
typeof("A")

# ---

# ## 基礎語法與概念 - 集合容器 (3/4)

# #### 字串 (Strings) (2/4)
# + 字元有其對應的整數
Int('A')
# + 字元也可以進行加減法運算
'A' + 1
# + 字元也可以比大小
'A' > 'a'

# ---

# ## 基礎語法與概念 - 集合容器 (3/4)

# #### 字串 (Strings) (3/4)
# > 字串具有不可變的特性，因此一旦字串被創造出來，就不能被修改，既不能變長也不能變短。因此對其進行的任何操作，實際上是會產出新的字串。
x = "Hello world!"

# > 也可以用3個雙引號刮起來，就可使用多行字串
"""
Hello
world
!
"""

# ---

# ## 基礎語法與概念 - 集合容器 (3/4)

# #### 字串 (Strings) (4/4)

# + 字串的串接
## 使用 * 將多個字串串接起來
x="1 + "
x * "2 = " * "3"

# + 字串的內插
x = "Today"
y = "Sunday"
string(x, "is", y)

# ---

# ## 基礎語法與概念 - 集合容器 (4/4)

# #### 字典 (Dictionary) (1/3)
# + 所謂__字典__是由一組鍵vs.值 (key vs. value) 對應關係所組成的資料結構。鍵必須是唯一的、不能重複，但值可以重複出現。
x = Dict("A" => 1, "B" => 2, "C" => 3)
## 透過鍵這個索引來取得值, 輸入x["A"]可以得到1

# + 要是給的鍵不存在，就會發生錯誤
#

# ---

# ## 基礎語法與概念 - 集合容器 (4/4)

# #### 字典 (Dictionary) (2/3)

# + 增加新的鍵-值對
x["D"] = 4
x

# + 複寫
x["A"] = 5
x

# ---

# ## 基礎語法與概念 - 集合容器 (4/4)

# #### 字典 (Dictionary) (3/3)

# + 迭代

for (k, v) in x
    println(k, "->", v)
end

# ---

# ## 學完基礎語法之後 -- 進階

# + 型別 (Types)
#   - 複合型別, 參數化型別, 抽象型別等
# + 巨集 (Macro)
#   - metaprogramming
# + 方法 (Methods)
#   - 多重分派
#   - 參數化方法
#   - 泛型

# ---

# ## 學完基礎語法之後 -- 實做

# + 加入Julia Language Study Group並完成其作業
# + 上線上課 (Udemy等)，寫作業
# + 參與別人的研究，透過協助他人分析資料來練功
# + 拿Paper所提供的資料重複跑一次分析

# ---

# # Acknowledgement
# + Taipei Bioinformatics Omnibus
#
# + 杜岳華
#


text:
"""
    kstest(class, var; rev = true)

Calculate empirical 2 sample Kolmogorov-Smirnov statistic and its location.\\
`class` is a 2 level categorical variable, `var` is the distribution to analyze.

Returns:
- n, total number of observations
- n1, number of observations of class 1
- n0, number of observations of class 0
- baserate, incidence rate of class 1
- ks, the maximum separation between the two cumulative distributions
- ksarg, the value of `var` at which maximum separation is achieved
- ksdep, depth of ksarg in the sorted values of `var`
`rev` = true counts depth from high value towards low value
"""
function kstest(class::BitArray{1}, var::Vector; rev = true)
    n = length(class)
    n == length(var) || error(ArgumentError("class and var should have the same length"))
    n1 = sum(class)
    n1 == 0 && error(ArgumentError("there are no class 1"))
    n0 = n - n1
    n0 == 0 && error(ArgumentError("there are no class 0"))

    baserate = n1 / n
    idx = sortperm(var, rev = rev)
    tgt = class[idx]
    cdf1 = cumsum(tgt) ./ n1            # TPR: True Positive Rate
    cdf0 = cumsum(.!tgt) ./ n0          # FPR: False Positive Rate
    sep = cdf1 .- cdf0
    ks, ksidx = findmax(sep)
    ksarg = var[idx[ksidx]]             # var @ max sep
    ksdep = ksidx / n                   # depth @ max sep

    (n = n, n1 = n1, n0 = n0, baserate = baserate, ks = ks, ksarg = ksarg, ksdep = ksdep)
end

function kstest(class::Vector, var::Vector; rev = true)
    uc = sort!(unique(class))
    length(uc) == 2 || error(ArgumentError("class must be 2 levels"))

    kstest(class .== uc[2], var; rev = rev)
end


text:
@testset "disjointSet.jl" begin
    @testset "union(disjointSet::DisjointSet, representative1::Int, representative2::Int)" begin
        disjointSet = SeaPearl.DisjointSet(5)
        SeaPearl.setUnion!(disjointSet, 2, 3)
        @test disjointSet.parent[2] == 3
        @test disjointSet.parent[3] == -2

        SeaPearl.setUnion!(disjointSet, 3, 1)
        @test disjointSet.parent[3] == -2
        @test disjointSet.parent[1] == 3

        SeaPearl.setUnion!(disjointSet, 4, 5)
        @test disjointSet.parent[4] == 5
        @test disjointSet.parent[5] == -2

        SeaPearl.setUnion!(disjointSet, 5, 3)
        @test disjointSet.parent[1] == 3
        @test disjointSet.parent[2] == 3
        @test disjointSet.parent[3] == -3
        @test disjointSet.parent[4] == 5
        @test disjointSet.parent[5] == 3
    end

    @testset "find_representative(disjointSet::DisjointSet, element)::Int" begin
        disjointSet = SeaPearl.DisjointSet(6)
        SeaPearl.setUnion!(disjointSet, 3, 4)
        @test disjointSet.parent[3] == 4
        @test disjointSet.parent[4] == -2

        SeaPearl.setUnion!(disjointSet, 5, 6)
        @test disjointSet.parent[5] == 6
        @test disjointSet.parent[6] == -2

        SeaPearl.setUnion!(disjointSet, 4, 6)
        @test disjointSet.parent[3] == 4
        @test disjointSet.parent[4] == 6
        @test disjointSet.parent[5] == 6
        @test disjointSet.parent[6] == -3

        representative = SeaPearl.findRepresentative!(disjointSet, 3)
        @test representative == 6
        @test disjointSet.parent[3] == 6

        representative = SeaPearl.findRepresentative!(disjointSet, 1)
        @test representative == 1
    end

    @testset "greatest!(disjointSet::DisjointSet, element::Int)::Int" begin
        disjointSet = SeaPearl.DisjointSet(6)
        SeaPearl.setUnion!(disjointSet, 3, 4)
        SeaPearl.setUnion!(disjointSet, 5, 6)
        SeaPearl.setUnion!(disjointSet, 4, 6)

        max = SeaPearl.greatest!(disjointSet, 1)
        @test max == 1

        max = SeaPearl.greatest!(disjointSet, 3)
        @test max == 6
        @test disjointSet.parent[3] == 6
    end

    @testset "smallest!(disjointSet::DisjointSet, element::Int)::Int" begin
        disjointSet = SeaPearl.DisjointSet(6)
        SeaPearl.setUnion!(disjointSet, 3, 4)
        SeaPearl.setUnion!(disjointSet, 5, 6)
        SeaPearl.setUnion!(disjointSet, 4, 6)

        min = SeaPearl.smallest!(disjointSet, 1)
        @test min == 1

        min = SeaPearl.smallest!(disjointSet, 3)
        @test min == 3
        @test disjointSet.parent[3] == 6
    end
end


text:
using Conda

ENV["PYTHON"] = ""  # Configure PyCall.jl to use Conda.jl's Python

# We need install a specific version of llvmlite which supports the same minor version of
# llvm used with Julia. https://github.com/numba/llvmlite#compatibility
const LLVM_VERSION = Base.libllvm_version
const LLVMLITE_VERSION = if LLVM_VERSION > v"7.0.0"
    v"0.29.0"
elseif LLVM_VERSION >= v"6.0.0"
    v"0.26.0"
elseif LLVM_VERSION >= v"5.0.0"
    v"0.22.0"
elseif LLVM_VERSION >= v"4.0.0"
    v"0.20.0"
elseif LLVM_VERSION >= v"3.9.0"
    v"0.17.0"
elseif LLVM_VERSION >= v"3.8.0"
    v"0.15.0"
elseif LLVM_VERSION >= v"3.7.0"
    v"0.12.1"
elseif LLVM_VERSION >= v"3.6.0"
    v"0.8.0"
elseif LLVM_VERSION >= v"3.5.0"
    v"0.5.1"
else
    error("Can't find a version of llvmlite that supports $LLVM_VERSION")
end

Conda.update()
Conda.add_channel("conda-forge")
Conda.add("numpy==1.17.0")
Conda.add_channel("numba")  # channel containing llvmlite
Conda.add("llvmlite==$LLVMLITE_VERSION")  # https://github.com/numba/llvmlite
Conda.add_channel("invenia")
Conda.add("pandapower==2.1.0")
@info "Verifying pandapower install..."
python_bin = joinpath(Conda.PYTHONDIR, "python")
run(`$python_bin -c "import pandapower, pandapower.networks, pandapower.topology, pandapower.plotting, pandapower.converter, pandapower.estimation"`)
run(`$python_bin -c "from pandapower.networks import case9; case9()"`)


text:
@testset "Induced measure constructor" begin
	pts = rand(3,1000)
	ϵF, ϵs = [4, 4, 4], [ [3, 5, 2], [5, 5, 5]]
	iv1 = inducedrectangularinvariantmeasure(pts, ϵF, ϵs[1])
	iv2 = inducedrectangularinvariantmeasure(pts, ϵF, ϵs[2])
	@test iv1 isa InducedRectangularInvariantMeasure
	@test iv2 isa InducedRectangularInvariantMeasure
	@test sum(iv1.measure_induced.dist) ≈ 1
	@test sum(iv2.measure_induced.dist) ≈ 1
end


text:
using Test
using LogicCircuits

@testset "MNIST Loader tests" begin
    data = LogicCircuits.IO.mnist();

    @test typeof(data) == XYDataset{Float32, UInt8}
    @test size(feature_matrix(train(data))) == (60000, 784)
    @test isnothing(valid(data))
    @test size(feature_matrix(test(data))) == (10000, 784)

    @test typeof(feature_matrix(train(data))) == Array{Float32,2};
    @test typeof(labels(train(data))) == Array{UInt8,1};

    data2 = LogicCircuits.IO.sampled_mnist();
    @test typeof(data2) == XDataset{Bool}

    @test size(feature_matrix(train(data2))) == (50000, 784)
    @test size(feature_matrix(test(data2))) == (10000, 784)
end


text:
function fixargs(args::Vector{Any}, stub::Any)
    n = length(args)
    data = Array{Any}(uninitialized, n)
    na = BitArray(uninitialized, n)
    for i in 1:n
        if args[i] == :missing || args[i] == :NA
            data[i] = stub
            na[i] = true
        else
            data[i] = args[i]
            na[i] = false
        end
    end
    return data, na
end

# We assume that data has at least one "value" that isn't missing
function findstub_vector(ex::Expr)
    if ex.args[1] != :missing && ex.args[1] != :NA
        return ex.args[1]
    end
    for i in 2:length(ex.args)
        if ex.args[i] != :missing && ex.args[i] != :NA
            return ex.args[i]
        end
    end
    return missing
end

# We assume that data has at least one "value" that isn't missing
function findstub_matrix(ex::Expr)
    if ex.args[1].args[1] != :missing && ex.args[1].args[1] != :NA
        return ex.args[1].args[1]
    end
    nrows = length(ex.args)
    for row in 1:nrows
        subex = ex.args[row]
        for i in 1:length(subex.args)
            if subex.args[i] != :missing && subex.args[i] != :NA
                return subex.args[i]
            end
        end
    end
    return missing
end

function parsevector(ex::Expr)
    if ex.head in (:ref, :typed_hcat, :typed_vcat)
        data, na = fixargs(ex.args[2:end], :(zero($(ex.args[1]))))
        return Expr(ex.head, ex.args[1], data...), ex.head == :typed_hcat ? reshape(na, 1, length(na)) : na
    else
        stub = findstub_vector(ex)
        data, na = fixargs(ex.args, stub)
        if ex.head == :hcat
            na = reshape(na, 1, length(na))
        end

        if isequal(stub, missing)
            return Expr(ex.head == :hcat ? (:typed_hcat) : (:ref), Any, data...), na
        else
            return Expr(ex.head, data...), na
        end
    end
end

function parsematrix(ex::Expr)
    if ex.head == :typed_vcat
        stub = :(zero($(ex.args[1])))
        rows = 2:length(ex.args)
    else
        stub = findstub_matrix(ex)
        rows = 1:length(ex.args)
    end

    nrows = length(rows)
    datarows = Array{Expr}(uninitialized, nrows)
    narows   = Array{Expr}(uninitialized, nrows)
    for irow in 1:nrows
        data, na = fixargs(ex.args[rows[irow]].args, stub)
        datarows[irow] = Expr(:row, data...)
        narows[irow] = Expr(:row, na...)
    end
    if ex.head == :typed_vcat
        return Expr(:typed_vcat, ex.args[1], datarows...), Expr(:vcat, narows...)
    elseif isequal(stub, missing)
        return Expr(:typed_vcat, Any, datarows...), Expr(:vcat, narows...)
    else
        return Expr(:vcat, datarows...), Expr(:vcat, narows...)
    end
end

function parsedata(ex::Expr)
    if length(ex.args) == 0
        return :([]), Expr(:call, :(Array{Bool}), 0)
    end
    if ex.head == :typed_vcat || (isa(ex.args[1], Expr) && ex.args[1].head == :row)
        return parsematrix(ex)
    else
        return parsevector(ex)
    end
end

"""
    @data expr

Create a [`DataArray`](@ref) based on the given expression.

# Examples

```jldoctest
julia> @data [1, missing, 3]
3-element DataArrays.DataArray{Int64,1}:
 1
  missing
 3

julia> @data hcat(1:3, 4:6)
3×2 DataArrays.DataArray{Int64,2}:
 1  4
 2  5
 3  6
```
"""
macro data(ex)
    if !(ex.head in (:vect, :vcat, :hcat, :ref, :typed_vcat, :typed_hcat))
        return quote
            tmp = $(esc(ex))
            DataArray(tmp, broadcast(x->isequal(x, missing), tmp))
        end
    end
    dataexpr, naexpr = parsedata(ex)
    return Expr(:call, :DataArray, esc(dataexpr), esc(naexpr))
end

"""
    @pdata expr

Create a [`PooledDataArray`](@ref) based on the given expression.

# Examples

```jldoctest
julia> @pdata ["Hello", missing, "World"]
3-element DataArrays.PooledDataArray{String,UInt32,1}:
 "Hello"
 missing
 "World"
```
"""
macro pdata(ex)
    if !(ex.head in (:vect, :vcat, :hcat, :ref, :typed_vcat, :typed_hcat))
        return quote
            tmp = $(esc(ex))
            PooledDataArray(tmp, broadcast(x->isequal(x, missing), tmp))
        end
    end
    dataexpr, naexpr = parsedata(ex)
    return Expr(:call, :PooledDataArray, esc(dataexpr), esc(naexpr))
end


text: @inline executor_type(x::Referenceables.Referenceable) = executor_type(parent(x))


text:
# Copyright 2020 Google LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

@testset "test_saddle_point" begin
  # Define problems that we will use for tests
  lp2 = deepcopy(example_lp())
  lp2.objective_vector = [0.0, 0.0, 0.0, 0.0]

  lp3 = deepcopy(example_lp())
  lp3.right_hand_side = [0.0, 0.0, 0.0]

  lp4 = deepcopy(lp3)
  lp4.variable_upper_bound = [Inf, Inf, Inf, Inf]

  lp5 = deepcopy(lp4)
  lp5.variable_lower_bound = [0.0, 0.0, 0.0, 1.0]

  qp = deepcopy(lp2)
  qp.objective_matrix[1, 1] = 1.0

  @testset "test_select_initial_primal_weight" begin
    verbosity = 0
    lp1 = example_lp()
    primal_importance = 1.3
    primal_weight = FirstOrderLp.select_initial_primal_weight(
      lp1,
      [1.0, 1.0, 1.0, 1.0],
      [1.0, 1.0, 1.0],
      primal_importance,
      verbosity,
    )
    @test primal_weight ==
          primal_importance * norm([5.0, 2.0, 1.0, 1.0], 2) /
          norm([12.0, 7.0, 1.0], 2)

    primal_weight = FirstOrderLp.select_initial_primal_weight(
      lp2,
      [1.0, 1.0, 1.0, 1.0],
      [1.0, 1.0, 1.0],
      primal_importance,
      verbosity,
    )
    @test primal_weight == primal_importance

    primal_weight = FirstOrderLp.select_initial_primal_weight(
      lp3,
      [1.0, 1.0, 1.0, 1.0],
      [1.0, 1.0, 1.0],
      primal_importance,
      verbosity,
    )
    @test primal_weight == primal_importance
  end

  @testset "test_compute_langragian_value" begin
    lp = example_lp()
    @test FirstOrderLp.compute_lagrangian_value(lp, zeros(4), zeros(3)) == -14.0
    qp = example_qp()
    @test FirstOrderLp.compute_lagrangian_value(qp, [1.0, 1.0], [0.0]) == 0.5
    @test FirstOrderLp.compute_lagrangian_value(qp, [1.0, 1.0], [1.0]) == 1.5
    @test FirstOrderLp.compute_lagrangian_value(qp, [0.25, 0.0], [0.0]) ==
          -0.125
  end
end


text:
# Benchmark multiple integration schemes against ground truth for comparison.
function benchmark_integration(integration_schemes::Array{AbstractIntegrator,1}, noise_levels::AbstractArray, shape = SynthSphere(); trials::Int = 5, debug_mode = false)
    # setup ground truth
    Zₜ = ground_truth(shape)
    Zₜ = Zₜ./maximum(Zₜ)

    # setup plotting
    scene = Scene()
    runs = first(size(noise_levels))
    labels = string.(typeof.(integration_schemes))
    colours = Array{RGB{Float64}, 1}(undef, 0)
    # Benchmark each scheme
    for i in eachindex(integration_schemes)

        # Setup error_list
        RMSE = zeros(Float64, 0)
        error_bars = Array{Pair{Point{2,Float32},Point{2,Float32}}, 1}(undef, 0)

        # For each noise level bechmark the scheme
        for j in eachindex(noise_levels)
            error_list = Array{Float64}(undef, 0)
            println("Part $i: $j/$runs")

            # Average across multiple trials
            for k = 1:trials
                # Setup and run test
                p, q = synthetic_gradient(shape, noise = noise_levels[j])
                Z = integration_schemes[i](p,q)
                Z = Z./maximum(Z)

                # Calculate error
                E,a = benchmark_shape(Z, Zₜ)
                push!(error_list, E)
            end

            # Average error and calculate error bars.
            temp_error = mean(error_list)
            push!(RMSE, temp_error)
            push!(error_bars, Point2f0(noise_levels[j], temp_error - std(error_list./(maximum(Zₜ) - minimum(Zₜ)))) => Point2f0(noise_levels[j], temp_error + std(error_list./(maximum(Zₜ) - minimum(Zₜ)))))
        end

        #plot result
        c = RGB{Float64}(rand(),rand(),rand())
        lines!(scene, noise_levels, RMSE, color = c)
        # linesegments!(scene, error_bars, color = c)
        push!(colours , c)
        a = [scene[2]]
        if i > 2
            for k = 3:(i+1)
                push!(a, scene[k])
            end
        end
        lgd = legend(a, labels[1:i]; camera=campixel!, raw=true)
        s = Makie.vbox(scene, lgd)
        if !debug_mode
            display(s)
        end
    end
    a = [scene[2]]
    for i = 3:(length(labels)+1)
        push!(a, scene[i])
    end
    lgd = legend(a, labels; camera=campixel!, raw=true)
    scene = Makie.vbox(scene, lgd)
    if !debug_mode
        display(s)
    end
end

# Benchmark reconsrtucted shape against ground truth and return error map
function benchmark_shape(Z::AbstractArray, ground::AbstractArray)
    # Setup test
    z = copy(Z)
    g = copy(ground)
    z = z./maximum(z)
    g = g./maximum(g)
    a = zeros(size(z))

    # Calculate error
    E = 0.0
    for i in CartesianIndices(Z)
        E = E + (g[i] - z[i])^2
        a[i] = sqrt((g[i] - z[i])^2)
    end
    error = sqrt(E / (first(size(Z)) * last(size(Z))))
    return error, a
end

function benchmark_image(img::AbstractArray, ground::AbstractArray)
    # Setup test
    R = Float64.(img)
    g = Float64.(ground)

    # Calculate error
    E = 0.0
    for i in CartesianIndices(R)
        E = E + (g[i] - R[i])^2
    end
    error = sqrt(E / (first(size(R)) * last(size(R))))
    return error
end

function benchmark_normals(P::AbstractArray, Q::AbstractArray, Pground::AbstractArray, Qground::AbstractArray)
    # Setup test
    p = copy(P)
    q = copy(Q)
    pground = copy(Pground)
    qground = copy(Qground)
    # Calculate error
    E = 0.0
    for i in CartesianIndices(p)
        E = E + (p[i] - pground[i])^2 + (q[i] - qground[i])^2
    end
    error = sqrt(E / (2 * first(size(p)) * last(size(q))))
    return error
end


text:
# order 2
h = [1//2, 1, 1, 1, 1, 1, 1]
e = [1//1, 0, 0, 0, 0, 0, 0]
d1 = [-3//2, 2, -1//2, 0, 0, 0, 0]
d2 = [1//1, -2, 1, 0, 0, 0, 0]
d3 = [-1//1, 3, -3, 1, 0, 0, 0]
H = Diagonal(h)

q1 = [0 1//2 0 0 0 0 0]
q2 = [-1//2 0 1//2 0 0 0 0]
q3 = [0 -1//2 0 1//2 0 0 0]
q4 = [0 0 -1//2 0 1//2 0 0]
q5 = [0 0 0 -1//2 0 1//2 0]
q6 = [0 0 0 0 -1//2 0 1//1]
q7 = [0 0 0 0 0 -1//2 0]
Q = vcat(q1, q2, q3, q4, q5, q6, q7)

m1 = [1//1 -1 0 0 0 0 0]
m2 = [m1[2] 2 -1 0 0 0 0]
m3 = [0 -1 2 -1 0 0 0]
m4 = [0 0 -1 2 -1 0 0]
m5 = [0 0 0 -1 2 -1 0]
m6 = [0 0 0 0 -1 2 -1]
m7 = [0 0 0 0 0 -1 2 ]
M = vcat(m1, m2, m3, m4, m5, m6, m7)

r1 = [0 -13//16 7//8 -1//16 0 0 0]
r2 = [-r1[2] 0 -23//16 5//8 0 0 0]
r3 = [-r1[3] -r2[3] 0 -17//16 1//2 0 0]
r4 = [0 -1//2 2//2 0 -2//2 1//2 0]
r5 = [0 0 -1//2 2//2 0 -2//2 1//2]
r6 = [0 0 0 -1//2 2//2 0 -2//2]
r7 = [0 0 0 0 -1//2 2//2 0]
R = vcat(r1, r2, r3, r4, r5, r6, r7)

n1 = [13//10 -12//5 9//10 1//5 0 0 0]
n2 = [n1[2] 26//5 -16//5 2//5 0 0 0]
n3 = [n1[3] n2[3] 47//10 -17//5 1 0 0]
n4 = [n1[4] n2[4] n3[4] -29//5 -4 1 0]
n5 = [0 0 1//1 -4 6 -4 1]
n6 = [0 0 0 1//1 -4 6 -4]
n7 = [0 0 0 0 1//1 -4 6]
N = vcat(n1, n2, n3, n4, n5, n6, n7)

D1 = H \ (Q - 1//2*e*e'); display(D1)
D2 = H \ (-M - e*d1'); display(D2)
D3 = H \ (R + 1//2*d1*d1' - e*d2'); display(D3)
D4 = H \ (N - e*d3' + d1*d2'); display(D4)

x=1//1:length(h)
ks = 0:5
for k in ks
    u=x.^k
    println("k = ", k)
    println(D1*u - k*x.^(k-1))
end
for k in ks
    u=x.^k
    println("k = ", k)
    println(D2*u - k*(k-1)*x.^(k-2))
end
for k in ks
    u=x.^k
    println("k = ", k)
    println(D3*u - k*(k-1)*(k-2)*x.^(k-3))
end
for k in ks
    u=x.^k
    println("k = ", k)
    println(D4*u - k*(k-1)*(k-2)*(k-3)*x.^(k-4))
end

# order 4
h = [35809//100800,
    13297//11200,
    5701//5600,
    45109//50400,
    35191//33600,
    33503//33600,
    1, 1, 1, 1]
e = [1//1, 0, 0, 0, 0, 0, 0, 0, 0, 0]
d1 = [-11//6, 18//6, -9//6, 2//6, 0, 0, 0, 0, 0, 0]
d2 = [2//1, -5, 4, -1, 0, 0, 0, 0, 0, 0]
d3 = [-1//1, 3, -3, 1, 0, 0, 0, 0, 0, 0]
H = Diagonal(h)

q1 = [0 526249//907200 -10819//777600 -50767//907200 -631//28800 91//7776 0 0 0 0]
q2 = [-q1[2] 0 1421209//2721600 16657//201600 -8467//453600 -33059//5443200 0 0 0 0]
q3 = [-q1[3] -q2[3] 0 631187//1360800 400139//5443200 -8789//302400 0 0 0 0]
q4 = [-q1[4] -q2[4] -q3[4] 0 496403//907200 -308533//5443200 0 0 0 0]
q5 = [-q1[5] -q2[5] -q3[5] -q4[5] 0 1805647//2721600 -1//12 0 0 0]
q6 = [-q1[6] -q2[6] -q3[6] -q4[6] -q5[6] 0 8//12 -1//12 0 0]
q7 = [0 0 0 0 1//12 -8//12 0 8//12 -1//12 0]
q8 = [0 0 0 0 0 1//12 -8//12 0 8//12 -1//12]
q9 = [0 0 0 0 0 0 1//12 -8//12 0 8//12]
q10 = [0 0 0 0 0 0 0 1//12 -8//12 0]
Q = vcat(q1, q2, q3, q4, q5, q6, q7, q8, q9, q10)

m1 = [2386127//2177280 -515449//453600 -10781//777600 61567//1360800 6817//403200 -1069//136080 0 0 0 0]
m2 = [m1[2] 4756039//2177280 -1270009//1360800 -3751//28800 3067//680400 119459//10886400 0 0 0 0]
m3 = [m1[3] m2[3] 111623//60480 -555587//680400 -551339//5443200 8789//453600 0 0 0 0]
m4 = [m1[4] m2[4] m3[4] 1025327//544320 -464003//453600 222133//5443200 0 0 0 0]
m5 = [m1[5] m2[5] m3[5] m4[5] 5074159//2177280 -1784047//1360800 1//12 0 0 0]
m6 = [m1[6] m2[6] m3[6] m4[6] m5[6] 1812749//725760 -16//12 1//12 0 0]
m7 = [0 0 0 0 1//12 -16//12 30//12 -16//12 1//12 0]
m8 = [0 0 0 0 0 1//12 -16//12 30//12 -16//12 1//12]
m9 = [0 0 0 0 0 0 1//12 -16//12 30//12 -16//12]
m10 = [0 0 0 0 0 0 0 1//12 -16//12 30//12]
M = vcat(m1, m2, m3, m4, m5, m6, m7, m8, m9, m10)

r1 = [0 -88471//67200 58139//33600 -1167//2800 -89//11200 7//640 0 0 0 0]
r2 = [-r1[2] 0 -43723//16800 46783//33600 -191//3200 -1567//33600 0 0 0 0]
r3 = [-r1[3] -r2[3] 0 -4049//2400 29083//33600 -71//1400 0 0 0 0]
r4 = [-r1[4] -r2[4] -r3[4] 0 -8591//5600 10613//11200 -1//8 0 0 0]
r5 = [-r1[5] -r2[5] -r3[5] -r4[5] 0 -108271//67200 8//8 -1//8 0 0]
r6 = [-r1[6] -r2[6] -r3[6] -r4[6] -r5[6] 0 -13//8 8//8 -1//8 0]
r7 = [0 0 0 1//8 -8//8 13//8 0 -13//8 8//8 -1//8]
r8 = [0 0 0 0 1//8 -8//8 13//8 0 -13//8 8//8]
r9 = [0 0 0 0 0 1//8 -8//8 13//8 0 -13//8]
r10 = [0 0 0 0 0 0 1//8 -8//8 13//8 0]
R = vcat(r1, r2, r3, r4, r5, r6, r7, r8, r9, r10)

n1 = [4596181//1814400 -10307743//1814400 160961//43200 -535019//907200 109057//1814400 -29273//604800 0 0 0 0]
n2 = [n1[2] 8368543//604800 -9558943//907200 2177057//907200 -11351//86400 204257//1814400 0 0 0 0]
n3 = [n1[3] n2[3] 4938581//453600 -786473//151200 1141057//907200 -120619//907200 0 0 0 0]
n4 = [n1[4] n2[4] n3[4] 3146581//453600 -4614143//907200 24587//14400 -1//6 0 0 0]
n5 = [n1[5] n2[5] n3[5] n4[5] 185709//22400 -11293343//1814400 12//6 -1//6 0 0]
n6 = [n1[6] n2[6] n3[6] n4[6] n5[6] 16787381//1814400 -39//6 12//6 -1//6 0]
n7 = [0 0 0 -1//6 12//6 -39//6 56//6 -39//6 12//6 -1//6]
n8 = [0 0 0 0 -1//6 12//6 -39//6 56//6 -39//6 12//6]
n9 = [0 0 0 0 0 -1//6 12//6 -39//6 56//6 -39//6]
n10 = [0 0 0 0 0 0 -1//6 12//6 -39//6 56//6]
N = vcat(n1, n2, n3, n4, n5, n6, n7, n8, n9, n10)

D1 = H \ (Q - 1//2*e*e'); display(D1)
D2 = H \ (-M - e*d1'); display(D2)
D3 = H \ (R + 1//2*d1*d1' - e*d2'); display(D3)
D4 = H \ (N - e*d3' + d1*d2'); display(D4)

x=1//1:length(h)
ks = 0:7
for k in ks
    u=x.^k
    println("k = ", k)
    println(Float64.(D1*u - k*x.^(k-1)))
end
for k in ks
    u=x.^k
    println("k = ", k)
    println(Float64.(D2*u - k*(k-1)*x.^(k-2)))
end
for k in ks
    u=x.^k
    println("k = ", k)
    println(Float64.(D3*u - k*(k-1)*(k-2)*x.^(k-3)))
end
for k in ks
    u=x.^k
    println("k = ", k)
    println(Float64.(D4*u - k*(k-1)*(k-2)*(k-3)*x.^(k-4)))
end

# order 6
h = [318365//1016064,
    145979//103680,
    139177//241920,
    964969//725760,
    593477//725760,
    52009//48384,
    141893//145152,
    1019713//1016064,
    1, 1, 1, 1]
e = [1//1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
d1 = [-25//12, 48//12, -36//12, 16//12, -3//12, 0, 0, 0, 0, 0, 0, 0]
d2 = [35//12, -104//12, 114//12, -56//12, 11//12, 0, 0, 0, 0, 0, 0, 0]
d3 = [-5//2, 18//2, -24//2, 14//2, -3//2, 0, 0, 0, 0, 0, 0, 0]
H = Diagonal(h)

q1 = reshape(
    [0, 1547358409//2421619200, -422423//11211200, -1002751721//8717829120,
    -15605263//484323840, 1023865//24216192, 291943739//21794572800, -24659//2534400,
    0, 0, 0, 0], 1, length(h))
q2 = reshape(
    [-q1[2], 0, 23031829//62899200, 10784027//34594560,
    2859215//31135104, -45982103//345945600, -26681//1182720, 538846039//21794572800,
    0, 0, 0, 0], 1, length(h))
q3 = reshape(
    [-q1[3], -q2[3], 0, 28368209//69189120,
    -9693137//69189120, 1289363//17740800, -39181//5491200, -168647//24216192,
    0, 0, 0, 0], 1, length(h))
q4 = reshape(
    [-q1[4], -q2[4], -q3[4], 0,
    5833151//10644480, 4353179//69189120, 2462459//155675520, -215471//10762752,
    0, 0, 0, 0], 1, length(h))
q5 = reshape(
    [-q1[5], -q2[5], -q3[5], -q4[5], 0, 7521509//13837824, -1013231//11531520, 103152839//8717829120,
    0, 0, 0, 0], 1, length(h))
q6 = reshape(
    [-q1[6], -q2[6], -q3[6], -q4[6], -q5[6], 0, 67795697//98841600, -17263733//151351200,
    1//60, 0, 0, 0], 1, length(h))
q7 = reshape(
    [-q1[7], -q2[7], -q3[7], -q4[7], -q5[7], -q6[7], 0, 1769933569//2421619200,
    -9//60, 1//60, 0, 0], 1, length(h))
q8 = reshape(
    [-q1[8], -q2[8], -q3[8], -q4[8], -q5[8], -q6[8], -q7[8], 0,
    45//60, -9//60, 1//60, 0], 1, length(h))
q9 = [0 0 0 0 0 -1//60 9//60 -45//60 0 45//60 -9//60 1//60]
q10 = [0 0 0 0 0 0 -1//60 9//60 -45//60 0 45//60 -9//60]
q11 = [0 0 0 0 0 0 0 -1//60 9//60 -45//60 0 45//60]
q12 = [0 0 0 0 0 0 0 0 -1//60 9//60 -45//60 0]
Q = vcat(q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12)

m1 = reshape(
    [4347276223//3736212480, -1534657609//1210809600, 68879//3057600, 1092927401//13076743680,
    18145423//968647680, -1143817//60540480, -355447739//65383718400, 56081//16473600,
    0, 0, 0, 0], 1, length(h))
m2 = reshape(
    [m1[2], 42416226217//18681062400, -228654119//345945600, -12245627//34594560,
    -2995295//46702656, 52836503//691891200, 119351//12812800, -634102039//65383718400,
    0, 0, 0, 0], 1, length(h))
m3 = reshape(
    [m1[3], m2[3], 5399287//4193280, -24739409//34594560,
    7878737//69189120, -1917829//31449600, 39727//3660800, 10259//4656960,
    0, 0, 0, 0], 1, length(h))
m4 = reshape(
    [m1[4], m2[4], m3[4], 7780367599//3736212480,
    -70085363//69189120, -500209//6289920, -311543//17962560, 278191//21525504,
    0, 0, 0, 0], 1, length(h))
m5 = reshape(
    [m1[5], m2[5], m3[5], m4[5],
    7116321131//3736212480, -545081//532224, 811631//11531520, -84101639//13076743680,
    0, 0, 0, 0], 1, length(h))
m6 = reshape(
    [m1[6], m2[6], m3[6], m4[6],
    m5[6], 324760747//138378240, -65995697//49420800, 1469203//13759200,
    -2//180, 0, 0, 0], 1, length(h))
m7 = reshape(
    [m1[7], m2[7], m3[7], m4[7],  m5[7], m6[7], 48284442317//18681062400, -1762877569//1210809600,
    27//180, -2//180, 0, 0], 1, length(h))
m8 = reshape(
    [m1[8], m2[8], m3[8], m4[8],  m5[8], m6[8], m7[8], 10117212851//3736212480,
    -270//180, 27//180, -2//180, 0], 1, length(h))
m9 = [0 0 0 0 0 -2//180 27//180 -270//180 490//180 -270//180 27//180 -2//180]
m10 = [0 0 0 0 0 0 -2//180 27//180 -270//180 490//180 -270//180 27//180]
m11 = [0 0 0 0 0 0 0 -2//180 27//180 -270//180 490//180 -270//180]
m12 = [0 0 0 0 0 0 0 0 -2//180 27//180 -270//180 490//180]
M = vcat(m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11, m12)

r1 = reshape(
    [0, -10882810591//5811886080, 398713069//132088320, -1746657571//1162377216,
    56050639//145297152, -11473393//1162377216, -38062741//1452971520, 30473//4392960,
    0, 0, 0, 0], 1, length(h))
r2 = reshape(
    [-r1[2], 0, -3720544343//830269440, 767707019//207567360,
    -1047978301//830269440, 1240729//14826240, 6807397//55351296, -50022767//1452971520,
    0, 0, 0, 0], 1, length(h))
r3 = reshape(
    [-r1[3], -r2[3], 0, -2870078009//830269440,
    74962049//29652480, -12944857//30750720, -17846623//103783680, 68707591//1162377216,
    0, 0, 0, 0], 1, length(h))
r4 = reshape(
    [-r1[4], -r2[4], -r3[4], 0,
    -727867087//276756480, 327603877//207567360, -175223717//830269440, 1353613//726485760,
    0, 0, 0, 0], 1, length(h))
r5 = reshape(
    [-r1[5], -r2[5], -r3[5], -r4[5],
    0, -1804641793//830269440, 311038417//207567360, -1932566239//5811886080,
    7//240, 0, 0, 0], 1, length(h))
r6 = reshape(
    [-r1[6], -r2[6], -r3[6], -r4[6], -r5[6], 0, -1760949511//830269440, 2105883973//1452971520,
    -72//240, 7//240, 0, 0], 1, length(h))
r7 = reshape(
    [-r1[7], -r2[7], -r3[7], -r4[7], -r5[7], -r6[7], 0, -1081094773//528353280,
    338//240, -72//240, 7//240, 0], 1, length(h))
r8 = reshape(
    [-r1[8], -r2[8], -r3[8], -r4[8], -r5[8], -r6[8], -r7[8], 0,
    -488//240, 338//240, -72//240, 7//240], 1, length(h))
r9 = [0 0 0 0 -7//240 72//240 -338//240 488//240 0 -488//240 338//240 -72//240]
r10 = [0 0 0 0 0 -7//240 72//240 -338//240 488//240 0 -488//240 338//240]
r11 = [0 0 0 0 0 0 -7//240 72//240 -338//240 488//240 0 -488//240]
r12 = [0 0 0 0 0 0 0 -7//240 72//240 -338//240 488//240 0]
R = vcat(r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11, r12)

n1 = reshape(
    [40833734273//10761070320, -162181998421//16397821440, 4696168417//521748864, -245714671483//68870850048,
    2185939219//2869618752, -15248255797//114784750080, 345156907//12298366080, 6388381//1093188096,
    0, 0, 0, 0], 1, length(h))
n2 = reshape(
    [n1[2], 147281127041//5380535160, -3072614435609//114784750080, 320122985851//28696187520,
    -768046031383//344354250240, 7861605187//14348093760, -803762437//4251287040, 167394281//86088562560,
    0, 0, 0, 0], 1, length(h))
n3 = reshape(
    [n1[3], n2[3], 139712483333//4782697920, -1634124842747//114784750080,
    90855193447//28696187520, -26412188989//38261583360, 668741173//1793511720, -132673781//2342545920,
    0, 0, 0, 0], 1, length(h))
n4 = reshape(
    [n1[4], n2[4], n3[4], 437353997177//43044281280,
    -172873969321//38261583360, 34759553483//28696187520, -98928859751//344354250240, 295000207//3587023440,
    0, 0, 0, 0], 1, length(h))
n5 = reshape(
    [n1[5], n2[5], n3[5], n4[5],
    126711914423//21522140640, -520477408939//114784750080, 49581230003//28696187520, -99640101991//344354250240,
    7//240, 0, 0, 0], 1, length(h))
n6 = reshape(
    [n1[6], n2[6], n3[6], n4[6],
    n5[6], 19422074929//2391348960, -772894368601//114784750080, 10579712849//4099455360,
    -96//240, 7//240, 0, 0], 1, length(h))
n7 = reshape(
    [n1[7], n2[7], n3[7], n4[7],
    n5[7], n6[7], 456715296239//43044281280, -915425403107//114784750080,
    676//240, -96//240, 7//240, 0], 1, length(h))
n8 = reshape(
    [n1[8], n2[8], n3[8], n4[8], n5[8], n6[8], n7[8], 488029542379//43044281280,
    -1952//240, 676//240, -96//240, 7//240], 1, length(h))
n9 = [0 0 0 0 7//240 -96//240 676//240 -1952//240 2730//240 -1952//240 676//240 -96//240]
n10 = [0 0 0 0 0 7//240 -96//240 676//240 -1952//240 2730//240 -1952//240 676//240]
n11 = [0 0 0 0 0 0 7//240 -96//240 676//240 -1952//240 2730//240 -1952//240]
n12 = [0 0 0 0 0 0 0 7//240 -96//240 676//240 -1952//240 2730//240]
N = vcat(n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11, n12)

D1 = H \ (Q - 1//2*e*e'); display(D1)
D2 = H \ (-M - e*d1'); display(D2)
D3 = H \ (R + 1//2*d1*d1' - e*d2'); display(D3)
D4 = H \ (N - e*d3' + d1*d2'); display(D4)


text:
# GCM-specific Sources
# This file contains helpers and lists currently available options

# Current options for GCM-specific sources:

"""
    HeldSuarezForcing{PV <: Union{Momentum,Energy}} <: TendencyDef{Source, PV}

Defines a forcing that parametrises radiative and frictional effects using
Newtonian relaxation and Rayleigh friction, following Held and Suarez (1994)
"""
struct HeldSuarezForcing{PV <: Union{Momentum, Energy}} <:
       TendencyDef{Source, PV} end

HeldSuarezForcing() =
    (HeldSuarezForcing{Momentum}(), HeldSuarezForcing{Energy}())

import ClimateMachine.BalanceLaws: source
import ClimateMachine.Atmos: filter_source, atmos_source!

filter_source(pv::PV, m, s::HeldSuarezForcing{PV}) where {PV} = s
atmos_source!(::HeldSuarezForcing, args...) = nothing

function held_suarez_forcing_coefficients(bl, args)
    @unpack state, aux = args
    @unpack ts = args.precomputed
    FT = eltype(state)

    # Parameters
    T_ref = FT(255)

    _R_d = FT(R_d(bl.param_set))
    _day = FT(day(bl.param_set))
    _grav = FT(grav(bl.param_set))
    _cp_d = FT(cp_d(bl.param_set))
    _p0 = FT(MSLP(bl.param_set))

    # Held-Suarez parameters
    k_a = FT(1 / (40 * _day))
    k_f = FT(1 / _day)
    k_s = FT(1 / (4 * _day))
    ΔT_y = FT(60)
    Δθ_z = FT(10)
    T_equator = FT(315)
    T_min = FT(200)
    σ_b = FT(7 / 10)

    # Held-Suarez forcing
    φ = latitude(bl, aux)
    p = air_pressure(ts)

    #TODO: replace _p0 with dynamic surface pressure in Δσ calculations to account
    #for topography, but leave unchanged for calculations of σ involved in T_equil
    σ = p / _p0
    exner_p = σ^(_R_d / _cp_d)
    Δσ = (σ - σ_b) / (1 - σ_b)
    height_factor = max(0, Δσ)
    T_equil = (T_equator - ΔT_y * sin(φ)^2 - Δθ_z * log(σ) * cos(φ)^2) * exner_p
    T_equil = max(T_min, T_equil)
    k_T = k_a + (k_s - k_a) * height_factor * cos(φ)^4
    k_v = k_f * height_factor
    return (k_v = k_v, k_T = k_T, T_equil = T_equil)
end

function source(s::HeldSuarezForcing{Energy}, m, args)
    @unpack state = args
    @unpack ts = args.precomputed
    nt = held_suarez_forcing_coefficients(m, args)
    FT = eltype(state)
    _cv_d = FT(cv_d(m.param_set))
    @unpack k_T, T_equil = nt
    T = air_temperature(ts)
    return -k_T * state.ρ * _cv_d * (T - T_equil)
end

function source(s::HeldSuarezForcing{Momentum}, m, args)
    nt = held_suarez_forcing_coefficients(m, args)
    return -nt.k_v * projection_tangential(m, args.aux, args.state.ρu)
end

max_forks_repo_name: jm-c/CLIMA


text:
@testset "remove characters" begin
    s = "Hello_World"
    new_s = PSI.remove_chars(s, "_")
    @test new_s == "HelloWorld"
end

@testset "test key with value" begin
    d = Dict("foo" => "bar")
    @test_throws ErrorException PSI.find_key_with_value(d, "fake")
end

@testset "dense axis to dataframe" begin
    one = JuMP.Containers.DenseAxisArray{Float64}(undef, 1:2)
    fill!(one, 1.0)
    one_df = PSI.axis_array_to_dataframe(one, [:name])
    test_df = DataFrames.DataFrame(:name => [1.0, 1.0])
    @test one_df == test_df
    three = JuMP.Containers.DenseAxisArray{Float64}(undef, [:a], 1:2, 1:3)
    fill!(three, 1.0)
    three_df = PSI.axis_array_to_dataframe(three)
    test_df = DataFrames.DataFrame(
        :S1 => [1.0, 1.0, 1.0, 2.0, 2.0, 2.0],
        :a => [1.0, 1.0, 1.0, 1.0, 1.0, 1.0],
    )
    @test three_df == test_df
    four = JuMP.Containers.DenseAxisArray{Float64}(undef, [:a], 1:2, 1:3, 1:5)
    @test_throws ErrorException PSI.axis_array_to_dataframe(four)

    three_int = JuMP.Containers.DenseAxisArray{Int}(undef, [:a], 1:2, 1:3)
    fill!(three_int, 2)
    three_df = PSI.axis_array_to_dataframe(three_int)
    test_df = DataFrames.DataFrame(
        :S1 => [1.0, 1.0, 1.0, 2.0, 2.0, 2.0],
        :a => [2.0, 2.0, 2.0, 2.0, 2.0, 2.0],
    )
    @test three_df == test_df
    four = JuMP.Containers.DenseAxisArray{Int}(undef, [:a], 1:2, 1:3, 1:5)
    fill!(four, 2)
    @test isempty(PSI.axis_array_to_dataframe(four))
end


text:
"""
    funnel(names, values)

Creates an `EChart` funnel.

## Methods
```julia
funnel(names::AbstractVector, values::AbstractVector{<:Union{Missing, Real}})
```

## Arguments
* `legend::Bool` : display legend?
* `kwargs` : varargs to set any field of resulting `EChart` struct

## Notes

Given the concept of a funnel (largest at top, narrowing), if you pass in the values out of order, the echarts.js library automatically sorts the numbers highest to lowest.

## Examples
```julia
v = [100,60,80,70,50]
n = ["A", "B", "C", "D", "E"]
fn = funnel(n, v)
```
"""
function funnel(names::AbstractVector, values::AbstractVector{<:Union{Missing, Real}};
 legend::Bool = false,
 sort::String = "descending",
 kwargs...)

    ec = newplot(kwargs, ec_charttype = "funnel")

    ec.series = [FunnelSeries(name = "Series 1",
                            data =  arrayofdicts(value = values, name = names),
                            sort = sort
                            )]

    #Add legend if requested
    legend ? legend!(ec) : nothing

    return ec

end


text:
function LineSegments(t::Tube{T})::Vector{AbstractLine{T,3,:cartesian}} where {T <: SSDFloat}
    ls = AbstractLine{T, 3, :cartesian}[]
    translate::CartesianVector{T} = ismissing(t.translate) ? CartesianVector{T}([0, 0, 0]) : t.translate
    for r in (t.r_interval.left == 0 ? [t.r_interval.right] : [t.r_interval.left, t.r_interval.right])
        for z in [t.z_interval.left, t.z_interval.right]
            push!(ls, PartialCircle(r, t.φ_interval.left, t.φ_interval.right, translate + CartesianVector{T}([0, 0, z])))
        end
    end
    for r in [t.r_interval.left, t.r_interval.right]
        if r != 0
            for φ in ((t.φ_interval.right - t.φ_interval.left ≈ 2π) ? [t.φ_interval.left] : [t.φ_interval.left, t.φ_interval.right])
                push!(ls, LineSegment(
                    CartesianPoint{T}(r * cos(φ), r * sin(φ), t.z_interval.left) + translate,
                    CartesianPoint{T}(r * cos(φ), r * sin(φ), t.z_interval.right) + translate))
            end
        end
    end
    for φ in ((t.φ_interval.right - t.φ_interval.left ≈ 2π) ? [] : [t.φ_interval.left, t.φ_interval.right])
        for z in [t.z_interval.left, t.z_interval.right]
            push!(ls, LineSegment(
                CartesianPoint{T}(t.r_interval.left * cos(φ), t.r_interval.left * sin(φ), z) + translate,
                CartesianPoint{T}(t.r_interval.right * cos(φ), t.r_interval.right * sin(φ), z) + translate))
        end
    end
    return ls
end

@recipe function f(t::Tube{T}; n = 30, seriescolor = :green) where {T}
    linewidth --> 2
    n --> n
    @series begin
        seriescolor --> seriescolor
        label --> "Tube"
        []
    end
    label := ""
    seriescolor := seriescolor
    LineSegments(t)
end

function LineSegments(c::Cone{T})::Vector{AbstractLine{T, 3, :cartesian}} where {T <: SSDFloat}
    ls = AbstractLine{T, 3, :cartesian}[]
    translate::CartesianVector{T} = ismissing(c.translate) ? CartesianVector{T}([0, 0, 0]) : c.translate
    for r in [c.rStart1, c.rStop1]
        push!(ls, PartialCircle(r, c.φStart, c.φStop, translate + CartesianVector{T}([0, 0, c.zStart])))
    end
    for r in [c.rStart2, c.rStop2]
        push!(ls, PartialCircle(r, c.φStart, c.φStop, translate + CartesianVector{T}([0, 0, c.zStop])))

    end
    for φ in ((c.φStop - c.φStart ≈ 2π) ? [c.φStart] : [c.φStart, c.φStop])
        push!(ls, LineSegment(
            CartesianPoint{T}(c.rStart1 * cos(φ), c.rStart1 * sin(φ), c.zStart) + translate,
            CartesianPoint{T}(c.rStart2 * cos(φ), c.rStart2 * sin(φ), c.zStop) + translate))
        push!(ls, LineSegment(
            CartesianPoint{T}(c.rStop1 * cos(φ), c.rStop1 * sin(φ), c.zStart) + translate,
            CartesianPoint{T}(c.rStop2 * cos(φ), c.rStop2 * sin(φ), c.zStop) + translate))
    end
    for φ in ((c.φStop - c.φStart ≈ 2π) ? [c.φStart] : [c.φStart, c.φStop])
        push!(ls, LineSegment(
            CartesianPoint{T}(c.rStart1 * cos(φ),  c.rStart1 * sin(φ),  c.zStart) + translate,
            CartesianPoint{T}(c.rStop1 * cos(φ),   c.rStop1 * sin(φ), c.zStart) + translate))
        push!(ls, LineSegment(
            CartesianPoint{T}(c.rStart2 * cos(φ),  c.rStart2 * sin(φ),  c.zStop) + translate,
            CartesianPoint{T}(c.rStop2 * cos(φ),   c.rStop2 * sin(φ), c.zStop) + translate))
    end
    return ls
end

@recipe function f(c::Cone{T}; n = 30, seriescolor = :orange) where {T}
    linewidth --> 2
    n --> n
    @series begin
        seriescolor --> seriescolor
        label --> "Cone"
        []
    end
    seriescolor := seriescolor
    label := ""
    LineSegments(c)
end


text:
include("../oracles/shortest_path_oracle.jl")

using Test

@testset "Shortest Path Oracle" begin
    sources = [1; 1; 2; 3]
    destinations = [2; 3; 4; 4]
    start_node = 1
    end_node = 4
    oracle = sp_flow_jump_setup(sources, destinations, start_node, end_node)

    weights = [4.0; 1; 4; 6]
    (z, w) = oracle(weights)
    @test z ≈ 7.0
    @test w ≈ [0, 1.0, 0, 1.0]

    weights2 = [2.5; 2.0; 1.1; 3.2]
    (z2, w2) = oracle(weights2)
    @test z2 ≈ 3.6
    @test w2 ≈ [1.0, 0, 1.0, 0]

    # Try a 5 edge case
    sources = [1; 1; 2; 3; 3]
    destinations = [2; 3; 4; 4; 2]
    start_node = 1
    end_node = 4
    oracle = sp_flow_jump_setup(sources, destinations, start_node, end_node)

    weights = [4.0; 1; 4; 6; 3]
    (z, w) = oracle(weights)
    @test z ≈ 7.0
    @test w ≈ [0, 1.0, 0, 1.0, 0]

    weights2 = [4.0; 1; 4; 6; 1]
    (z2, w2) = oracle(weights2)
    @test z2 ≈ 6.0
    @test w2 ≈ [0, 1.0, 1.0, 0, 1.0]

    # Big grid test
    sources, destinations = convert_grid_to_list(20,20)
    start_node, end_node = 1, 400
    oracle = sp_flow_jump_setup(sources, destinations, start_node, end_node)
    weights = ones(1520)
    (z, w) = oracle(weights)
    @test z ≈ 38.0
end

max_stars_repo_name: paulgrigas/SmartPredictThenOptimize

max_issues_repo_name: paulgrigas/SmartPredictThenOptimize

max_forks_repo_name: paulgrigas/SmartPredictThenOptimize


text:
using Distributions, Plots; pyplot()

mu, sig = 2, 3
eta = sqrt(3)*sig/pi
n, N = 15, 10^7
dNormal   = Normal(mu, sig)
dLogistic = Logistic(mu, eta)
xGrid = -8:0.1:12

sNormal   = [var(rand(dNormal,n)) for _ in 1:N]
sLogistic = [var(rand(dLogistic,n)) for _ in 1:N]

p1 = plot(xGrid, pdf.(dNormal,xGrid), c=:blue, label="Normal")
p1 = plot!(xGrid, pdf.(dLogistic,xGrid), c=:red, label="Logistic",
	xlabel="x",ylabel="Density", xlims= (-8,12), ylims=(0,0.16))

p2 = stephist(sNormal, bins=200, c=:blue, normed=true, label="Normal")
p2 = stephist!(sLogistic, bins=200, c=:red, normed=true, label="Logistic",
	xlabel="Sample Variance", ylabel="Density", xlims=(0,30), ylims=(0,0.14))

plot(p1, p2, size=(800, 400))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using Base.Test

using POMDPs
import POMDPs: transition, reward, initial_state_distribution
import POMDPs: generate_o
type A <: POMDP{Int,Bool,Bool} end

println("Warning expected:")
@test_throws MethodError initial_state(A(), Base.GLOBAL_RNG)
println("Warning expected:")
@test_throws MethodError generate_s(A(), 1, true, Base.GLOBAL_RNG)

type B <: POMDP{Int, Bool, Bool} end

transition(b::B, s::Int, a::Bool) = Int[s+a]
@test implemented(generate_s, Tuple{B, Int, Bool, MersenneTwister})
@test generate_s(B(), 1, false, Base.GLOBAL_RNG) == 1

@test !@implemented generate_sor(::B, ::Int, ::Bool, ::MersenneTwister)
# @test_throws MethodError generate_sor(B(), 1, false, Base.GLOBAL_RNG)

reward(b::B, s::Int, a::Bool, sp::Int) = -1.0
@test generate_sr(B(), 1, false, Base.GLOBAL_RNG) == (1, -1.0)

generate_o(b::B, s::Int, a::Bool, sp::Int, rng::AbstractRNG) = sp
@test @implemented generate_o(b::B, s::Int, a::Bool, sp::Int, rng::AbstractRNG)
@test @implemented generate_so(b::B, s::Int, a::Bool, rng::MersenneTwister)
@test @implemented generate_sor(b::B, s::Int, a::Bool, rng::MersenneTwister)
@test generate_sor(B(), 1, true, Base.GLOBAL_RNG) == (2, 2, -1.0) # should throw sor error

initial_state_distribution(b::B) = Int[1,2,3]
@test initial_state(B(), Base.GLOBAL_RNG) in initial_state_distribution(B())


text:
typeof_model_elements(model::StateSpaceModel) = eltype(model.system.y)
num_states(model::StateSpaceModel) = num_states(system(model))
num_series(model::StateSpaceModel) = num_series(system(model))
system(model::StateSpaceModel) = model.system
isunivariate(model::StateSpaceModel) = isa(model.system.y, Vector)
model_name(model::StateSpaceModel) = "$(typeof(model))"
length_observations(model::StateSpaceModel) = length(model.system.y)
observations(model::StateSpaceModel) = model.system.y

function lagmat(y::Vector{Fl}, k::Int) where Fl
    X = Matrix{Fl}(undef, length(y) - k, k)
    for i in 1:k
        X[:, i] = lag(y, i)[k + 1:end]
    end
    return X
end


text:
const plot_functions = Dict(
    :plot => [plot, (args...; kwargs...) -> plot!(shared.plt, args...; kwargs...)],
    :scatter => [scatter, (args...; kwargs...) -> scatter!(shared.plt, args...; kwargs...)],
    :groupedbar => [groupedbar, (args...; kwargs...) -> groupedbar!(shared.plt, args...; kwargs...)],
    :boxplot => [boxplot, (args...; kwargs...) -> boxplot!(shared.plt, args...; kwargs...)],
    :violin => [violin, (args...; kwargs...) -> violin!(shared.plt, args...; kwargs...)],
    :marginalhist => [marginalhist, marginalhist]
)

function get_plotvalues(df)
    xvalues = ComboBoxType("X AXIS", ComboBoxEntry.(string.(colnames(df))), true)
    ylist = string.(union([:hazard, :density, :cumulative],colnames(df)))
    yvalues = ComboBoxType("Y AXIS", ComboBoxEntry.(ylist), true)
    plotlist = [
        "plot",
        "scatter",
        "groupedbar",
        "boxplot",
        "violin",
        "marginalhist"
    ]
    plot_type = ComboBoxType("PLOT TYPE", ComboBoxEntry.(plotlist), false)
    axislist = ["pointbypoint", "continuous", "binned", "discrete"]
    axis_type = ComboBoxType("AXIS TYPE", ComboBoxEntry.(axislist) , false)
    errorlist = union(["none", "bootstrap", "across"], "across " .* string.(colnames(df)))
    compute_error = ComboBoxType("COMPUTE ERROR",  ComboBoxEntry.(errorlist), false)
    namelist = union(["pointbypoint"], "point=" .* string.(colnames(df)))
    dataperpoint = ComboBoxType("DATA PER POINT",  ComboBoxEntry.(namelist), false)
    return [xvalues, yvalues, plot_type, axis_type, compute_error, dataperpoint]
end

get_kwargs(s) = s == "" ? [] : [(x.args[1], eval(x.args[2])) for x in parse("("*s*",)").args]

get_func(s) = (s == "") ? mean : eval(parse(s))

function get_plot!(shared, in_place)
    df, selectlist, plotvalues = shared.df, shared.selectlist, shared.plotvalues
    selectdata = choose_data(shared)
    xval, yval, line, axis_type, compute_error, dataperpoint =
        getfield.(plotvalues, :chosen_value)
    extra_kwargs = get_kwargs(shared.plotkwargs.value)
    x_info, y_info = plotvalues[1].text_info, plotvalues[2].text_info
    xfunc = get_func(x_info)
    yfunc = get_func(y_info)
    isrecipe = !(Symbol(line) in [:plot, :scatter, :groupedbar])
    isgroupapply = !isrecipe && !(Symbol(axis_type) == :pointbypoint)
    xlabel, ylabel = xval, yval
    plot_func = (args...; kwargs...) -> plot_functions[Symbol(line)][in_place ? 2 : 1](args...;
                xlabel = xlabel, ylabel = ylabel, kwargs..., extra_kwargs...)
    group_cols = [Symbol(col.name) for col in selectlist if col.split]

    maybe_nbins = Symbol(axis_type) == :binned ? (round(Int64, 101-shared.smoother.value),) : ()

    s = GroupedErrors.ColumnSelector(selectdata)
    s = GroupedErrors._splitby(s, group_cols)
    if (Symbol(axis_type) == :pointbypoint) || isrecipe
        if !isrecipe && ('=' in dataperpoint)
            s = GroupedErrors._across(s, Symbol(split(dataperpoint, '=')[2]))
            if check_equality() && (':' in shared.splitting_var.chosen_value)
                splitting_var = Symbol(split(shared.splitting_var.chosen_value, ':')[2])
                s = GroupedErrors._compare(s, splitting_var)
            end
            s = GroupedErrors._x(s, Symbol(xval), xfunc)
            s = GroupedErrors._y(s, Symbol(yval), yfunc)
        else isrecipe
            s = GroupedErrors._x(s, Symbol(xval))
            s = GroupedErrors._y(s, Symbol(yval))
        end
    else
        compute_error = convert_error_type(compute_error)
        if compute_error[1] == :across
            s = GroupedErrors._across(s, compute_error[2])
        elseif compute_error[1] == :bootstrap
            s = GroupedErrors._bootstrap(s, compute_error[2])
        end
        s = GroupedErrors._x(s, Symbol(xval), Symbol(axis_type), maybe_nbins...)
        if Symbol(axis_type) == :continuous
            if Symbol(yval) in colnames(selectdata)
                s = GroupedErrors._y(s, :locreg, Symbol(yval), span = (shared.smoother.value+1.0)/100)
            elseif Symbol(yval) in [:density, :hazard]
                s = GroupedErrors._y(s, Symbol(yval), bandwidth = (shared.smoother.value+1.0)*std(column(selectdata, Symbol(xval)))/200)
            else
                s = GroupedErrors._y(s, Symbol(yval))
            end
        else
            s = GroupedErrors._x(s, Symbol(xval), xfunc)
            if Symbol(yval) in colnames(selectdata)
                s = GroupedErrors._y(s, :locreg, Symbol(yval), estimator = yfunc)
            else
                s = GroupedErrors._y(s, Symbol(yval))
            end
        end
    end
    shared.plt = Symbol(line) == :plot ? @plot(s, plot_func(), :ribbon) : @plot(s, plot_func())
    get(s.kw, :compare, false) && Plots.abline!(shared.plt, 1, 0, legend = false)
end

function convert_error_type(s::AbstractString)
    n = search(s, ' ')
    if n == 0
        if Symbol(s) == :bootstrap
            return (:bootstrap, 1000)
        elseif Symbol(s) == :across
            return (:across, :all)
        else
            return (Symbol(s),)
        end
    else
        return (Symbol(s[1:n-1]), Symbol(s[n+1:end]))
    end
end


text:
function hvcat!(C, Tiles)
    libgb.GxB_Matrix_concat(C, Tiles, size(Tiles,1), size(Tiles,2), C_NULL)
    return C
end

function Base.hvcat(Tiles)
    ncols = sum(size.(Tiles[1,:], 2))
    nrows = sum(size.(Tiles[:, 1], 1))
    types = eltype.(Tiles)
    t = types[1]
    for type ∈ types[2:end]
        t = promote_type(t, type)
    end
    if Tiles isa AbstractArray{<:GBVector}
        Tiles = GBMatrix.(Tiles)
    end
    C = GBMatrix{t}(nrows,ncols)
    return hvcat!(C, Tiles)
end

vcat!(C, A::GBArray...) = hvcat!(C, collect(A))
Base.vcat(A::GBArray...) = hvcat(collect(A))

hcat!(C, A::GBArray...) = hvcat!(C, reshape(collect(A), 1, :))
Base.hcat(A::GBArray...) = hvcat(reshape(collect(A), 1, :))


text:
function TXTLDictionary(time_start,time_stop,time_step)

#Blocked Reactions
degradation_rxn = [100:1:116;];

AA_synthesis_rxn = [81;87;88;89;90;91;92;93;94;96;97;98;99;];

# setup promoter models -
Promoter = Float64[]
push!(Promoter,10) #K1

# Setup the mRNA elongation rate, and global translation
RNAP_concentration_nM = 75;                          # 60-75nM (ACS SynBio Garamella 2016)
RNAP_elongation_rate = 25;                          # >5 NT/s (ACS SynBio Garamella 2016)
RIBOSOME_concentration = 0.0016;                      # 0.0016mM with 72% MaxActive (Underwood, Swartz, Puglisi 2005 Biotech Bioeng) & <0.0023mM (ACS SynBio Garamella 2016)
RIBOSOME_elongation_rate = 2;                      # >1 (ACS SynBio Garamella 2016) & 1.5 AA/sec (Underwood, Swartz, Puglisi 2005 Biotech Bioeng)
#number_of_cells = 2e8;                                # 1e9 cells/ml

# Protein of Interest Sequences
mRNA_length = 660;
protein_length = 219;
gene_copies = 3.125e10;
mRNA_degradation_rate = 10#5.2;
polysome_gain = 5#10;
volume = 15e-6
inducer = 35;
plasmid_saturation_coefficient = 3.5;
mRNA_saturation_coefficient = 0.045;

# =============================== DO NOT EDIT BELOW THIS LINE ============================== #
TXTL_dictionary = Dict{AbstractString,Any}()
TXTL_dictionary["AA_degradation_rxn"] = degradation_rxn
TXTL_dictionary["AA_synthesis_rxn"] = AA_synthesis_rxn
TXTL_dictionary["Promoter"] = Promoter
TXTL_dictionary["RNAP_concentration_nM"] = RNAP_concentration_nM
TXTL_dictionary["RNAP_elongation_rate"] = RNAP_elongation_rate
TXTL_dictionary["RIBOSOME_concentration"] = RIBOSOME_concentration
TXTL_dictionary["RIBOSOME_elongation_rate"] = RIBOSOME_elongation_rate
TXTL_dictionary["mRNA_length"] = mRNA_length
TXTL_dictionary["protein_length"] = protein_length
TXTL_dictionary["gene_copies"] = gene_copies
TXTL_dictionary["mRNA_degradation_rate"] = mRNA_degradation_rate
TXTL_dictionary["polysome_gain"] = polysome_gain
TXTL_dictionary["volume"] = volume
TXTL_dictionary["inducer"] = inducer
TXTL_dictionary["plasmid_saturation_coefficient"] = plasmid_saturation_coefficient
TXTL_dictionary["mRNA_saturation_coefficient"] = mRNA_saturation_coefficient
# =============================== DO NOT EDIT ABOVE THIS LINE ============================== #
return TXTL_dictionary
end


text:
function fib(n)
	if n <= 1
		return n
	end
	return fib(n - 2) + fib(n - 1)
end

for i in 0:24
	fib(i)
end

max_stars_repo_name: bichanna/slap

max_issues_repo_name: adiyenuubarii/slap

max_forks_repo_name: adiyenuubarii/slap


text:
module TracedRandom

using Random

export here, rand!, randn!, randexp, randexp!,
       bitrand, randstring, randsubseq, randsubseq!,
       shuffle, shuffle!, randperm, randperm!, randcycle, randcycle!

const Address = Union{Symbol,Pair{Symbol}}

Base.rand(::Address, args...) = rand(args...)
Base.rand(rng::AbstractRNG, ::Address, args...) = rand(rng, args...)

Base.randn(::Address, args...) = randn(args...)
Base.randn(rng::AbstractRNG, ::Address, args...) = randn(rng, args...)

for fn in (:rand!, :randn!, :randexp, :randexp!, :bitrand, :randstring,
           :randsubseq, :randsubseq!, :shuffle, :shuffle!,
           :randperm, :randperm!, :randcycle, :randcycle!)
    fn = GlobalRef(Random, fn)
    @eval $fn(::Address, args...) = $fn(args...)
    @eval $fn(rng::AbstractRNG, ::Address, args...) = $fn(rng, args...)
end

for fn in (:rand, :rand!, :randn!, :randexp, :randexp!, :bitrand, :randstring,
           :randsubseq, :randsubseq!, :shuffle, :shuffle!,
           :randperm, :randperm!, :randcycle, :randcycle!)
    @eval begin
    @doc """
        $($fn)([rng=GLOBAL_RNG,] addr::Address, args...)

    Traced execution of `$($fn)`, where `addr` specifies the name / address of
    the random choice. By default, the address `addr` is ignored, but it can
    be intercepted to support inference in a probabilistic programming system.
    An `Address` is either a `Symbol`, or a `Pair` that begins with a `Symbol`.
    """
    $fn
    end
end

"Singleton type for the special [`here`](@ref) address."
struct Here end

"""
    here

Special address that refers to the address namespace of the calling context.
Supported when [`rand`](@ref) is called on a non-primitive stochastic function,
causing all named random variables sampled by that function to be spliced into
the address namespace of the calling context.
"""
const here = Here()

"""
    rand([rng=GLOBAL_RNG,] addr::Union{Here,Address}, fn::Function, args...)

Traced execution of a stochastic function `fn`, where `addr` specifies the
address namespace for any primitive random choices made within the call to `fn`.
By default, the address `addr` is ignored, but it can be intercepted to support
inference in a probabilistic programming system.

An `Address` is either a `Symbol`, or a `Pair` that begins with a `Symbol`.
Alternatively, if the special address [`here`](@ref) is specified, then all
traced random choices made by `fn` will be spliced into the namespace of the
calling context.
"""
Base.rand(::Here, fn::Function, args...) = fn(args...)
Base.rand(::AbstractRNG, ::Here, fn::Function, args...) = fn(args...)

Base.rand(::Address, fn::Function, args...) = fn(args...)
Base.rand(::AbstractRNG, ::Address, fn::Function, args...) = fn(args...)

end


text:
struct RuntimeArrays{T}
    # Iterparticle Arrays
    indij       :: Vector{SVector{2, Int}}
    rij         :: Vector{SVector{3, <:T}}
    inside      :: Vector{<:T}
    wj          :: Vector{<:T}
    rcutij      :: Vector{<:T}
    element     :: Vector{Int}

    # U Arrays: Real, Imaginary Parts
    ulisttot_r  :: Vector{<:T}
    ulisttot_i  :: Vector{<:T}
    dulist_r    :: Matrix{<:T}
    dulist_i    :: Matrix{<:T}

    ulist_r_ij  :: Matrix{<:T}
    ulist_i_ij  :: Matrix{<:T}

    # Z Arrays: Real, Imaginary Parts
    zlist_r     :: Vector{<:T}
    zlist_i     :: Vector{<:T}

    # Y Arrays: Real, Imaginary Parts
    ylist_r     :: Vector{<:T}
    ylist_i     :: Vector{<:T}

    # Bispectrum Arrays
    blist       :: Vector{<:T}
    dblist      :: Matrix{<:T}
end

function initialize_runtime_arrays(i::Int, nnlist::NeighborList, A::AbstractSystem, snap :: SNAPParams)
    # Initialize Arrays
    ind_ij  = SVector{2, Int}[]  # Contains (i,j) pairs
    rcut_ij = AbstractFloat[]   # Contains rcut between (i,j)
    inside  = AbstractFloat[]   # Unused
    wj      = AbstractFloat[]   # Weight of particle (j) w.r.t particle (i)
    rij     = SVector{3, AbstractFloat}[] # vector from particle j to particle i
    element = Int[]            # Element of particle (j)
    elements = snap.elements

    i_list = nnlist.i[i]
    j = nnlist.j[i]
    r = nnlist.r[i]
    i_element = findall(x->x==Symbol(A[i].atomic_symbol), elements)[1]
    for (jj, rj) in zip(j, r)
        push!(ind_ij, SVector{2}([i, jj]))

        j_element = findall(x->x==Symbol(A[jj].atomic_symbol), elements)[1]
        push!(element, j_element)
        push!(rcut_ij, (snap.radii[i_element] + snap.radii[j_element])*snap.rcutfac)
        push!(wj, snap.weight[j_element])
        push!(inside, 1)
        push!(rij, rj)
    end

    n_interactions = length(rij)
    num_elements = snap.prebuilt_arrays.num_elements
    n_doubles = num_elements * num_elements
    n_triples = n_doubles * num_elements
    return RuntimeArrays{AbstractFloat}(
        ind_ij,                                                 # indij
        rij,                                                    # rij
        inside,                                                 # inside
        wj,                                                     # wj
        rcut_ij,                                                 # rcutij
        element,                                                # element
        zeros(snap.prebuilt_arrays.idxu_max*num_elements),   # ulisttot_r
        zeros(snap.prebuilt_arrays.idxu_max*num_elements),   # ulisttot_i
        zeros(snap.prebuilt_arrays.idxu_max, 3),                # dulist_r
        zeros(snap.prebuilt_arrays.idxu_max, 3),                # dulist_i
        zeros(n_interactions+1, snap.prebuilt_arrays.idxu_max),                       # ulist_r_ij
        zeros(n_interactions+1, snap.prebuilt_arrays.idxu_max),                       # ulist_i_ij
        zeros(snap.prebuilt_arrays.idxz_max*n_doubles),         # zlist_r
        zeros(snap.prebuilt_arrays.idxz_max*n_doubles),         # zlist_i
        zeros(snap.prebuilt_arrays.idxu_max*num_elements),   # ylist_r
        zeros(snap.prebuilt_arrays.idxu_max*num_elements),   # ylist_i
        zeros(snap.prebuilt_arrays.idxb_max*n_triples),         # blist
        zeros(snap.prebuilt_arrays.idxb_max*n_triples, 3)       # dblist
    )
end


text:
include("apply_noise.jl")
include("apply_dyson.jl")

    """
        (p, t, pmax, tmax) = spatialsearch(psi0, Adjacency; kwargs...)

    Evaluate a (noisy) spatial search run on the graph described by the adjacency
    matrix Adjacency, starting from an initial state `psi0`. Return a `NamedTuple`
    consisting of
    - `p`: the probability of finding the walker in the target node as a
    f unction of `t`
    - `t`: the time instants at which probArray is evaluated
    - `pmax`: the maximum probability
    - `tmax`: the optimal time at which `pmax` is reached for the first time

    This function implements the spatial search as devised by Childs & Goldstone,
    Phys. Rev. A 70, 22314 (2004), with the addition of RTN noise as described in
    the manuscript

    # Arguments
    ## Required
    - `psi0`: a complex vector representing the initial state
    - `Adjacency`: a function that returns a sparse adjacency matrix for the desired
        graph topology

    ## Keyword arguments
    - `maxdt::Real = 0.05`: maximum allowed dt for the output
    - `gamma::Real = 1.`: Coupling γ in the Hamiltonian
    - `posW::Integer = 1`: Position of the target node
    - `time::Real = 1.`: Final time
    - `mu::Real = 1.`: Switching rate μ of the noise
    - `noiseStrength::Real = 0.`: Coupling ν to the noise (noise stregth)
    - `noiseRealizations::Integer = 1`: Number of Montecarlo noise realizations
    - `dysonOrder::Integer = 4`: Order of the Dyson series expansion

    # Examples

    In this example we run the noiseless spatial search on a complete graph with
    `N = 4` and we check that `p_max = 1` and `t_opt = 1` (up to error due to
    `dt`)

    ```jldoctest
    julia> (probArray, t, p_max, t_opt) = spatialsearch(superposition_state(4),
    complete_graph_Ad; time=5, gamma = 1. / 4, maxdt = .001)

    julia> isapprox(p_max, 1, atol=1e-6)
    true

    julia> isapprox(t_opt, pi, atol=1e-3)
    true
    ```
    """
    function spatialsearch(psi0, Adjacency;
                        maxdt::Real=0.05, # maximum allowed time for dt
                        gamma::Real=1.,
                        posW::Integer=1,
                        time::Real=1.,
                        mu::Real=1.,
                        coupling::Real=1.,
                        noiseStrength::Real=0.,
                        noiseRealizations::Integer=1,
                        dysonOrder::Integer=4)

        # maxdt = maximum value of dt of the evolution.
        # posW = position of the state we want to find in the basis of the nodes
        # of the graph (position in the state array).

        N = length(psi0)   # N (dimension of the lattice)

        # Required number of output timesteps
        output_timesteps::Int = Int(cld(time, maxdt) + 1)

        # We evaluate the number of internal steps required so that the noise is
        # sampled correctly. kout is the number of output timesteps t
        (Mtot, kout) = internal_M(time, mu, output_timesteps)

        # every internal step has a length = dt
        dtev = internal_dt(time, mu, output_timesteps)

        Ad = Adjacency(N)  # We define the whole adjacency matrix

        # We obtain the Hamiltonian (Laplacian matrix) from the adjacency matrix
        H = discLapl(Ad)

        # The number of links is the number of nonzero elements in the upper
        # triangular adjacency matrix.
        link_number = Int(length(Ad.nzval) / 2)
        tvec = range(0., stop=time, length=Mtot)

        # Contains the output
        probArray = zeros(Mtot)

        probArray = @distributed (+) for n = 1 : noiseRealizations

            # array containing the probability of getting the required state |posW>
            # if we measure the evoluted state, at each internal timestep t
            tmp = zeros(Float64, Mtot)
            tmp[1] = abs2(psi0[posW]) / noiseRealizations

            psi = copy(psi0)  #psi0 is the initial state (array of complex floats)

            # This is the noise array generated for each indipedent link.
            # It's a Mtot x NumberOfIndipendentLinks matrix.
            temparray = - gamma * (coupling .+
                        noiseStrength .* generateRTN(tvec, mu, link_number))

            for t = 2 : Mtot    # evolution over time
                apply_noise!(H, Ad, view(temparray, t-1, :))

                # we subtract the oracle Hamiltonian from the noisy Laplacian
                H[posW, posW] -= 1.

                apply_dyson!(H, psi, dtev, dysonOrder)

                psi ./= norm(psi)     # state psi at time=t

                # we calculate the probability
                tmp[t] = abs2(psi[posW]) ./ noiseRealizations

            end  # end of loop over time

            tmp # The parallel-for reduction sums this value to probArray
        end  # end of loop over noise realizations

        # we find the maximum and its position in the probability array
        (max, position) = findmax(probArray)

        return  (p=probArray, t=tvec, pmax=max, tmax=(position - 1) * dtev)
    end

max_stars_repo_name: matteoacrossi/QuantumSpatialSearch

max_issues_repo_name: matteoacrossi/QuantumSpatialSearch

max_forks_repo_name: matteoacrossi/QuantumSpatialSearch


text:
# Abstract type for MIDI events
abstract TrackEvent

#=
All track events begin with a variable length time value (see readvariablelength()) and must have a field named dT.
MIDI events start with a MIDI channel message defined in constants.jl. They're followed by 1 or 2 bytes, depending on the
channel message (see EVENTTYPETOLENGTH). If no valid channel message is identified, the previous seen channel message is used.
Meta events and sysex events both begin with a specific byte (see constants.jl)
=#

function dt(e::TrackEvent)
    e.dT
end


text:
struct Pca
    T::Array{Float64}
    P::Array{Float64}
    sv::Vector{Float64}
    xmeans::Vector{Float64}
    weights::Vector{Float64}
    ## For consistency with PCA Nipals
    niter::Union{Int64, Nothing}
    conv::Union{Bool, Nothing}
end

"""
    pcasvd(X, weights = ones(size(X, 1)); nlv)
PCA by SVD factorization.
* `X` : X-data.
* `weights` : Weights of the observations.
* `nlv` : Nb. principal components (PCs).

Noting D a (n, n) diagonal matrix of weights for the observations (rows of X),
the function makes a SVD factorization of D^(1/2) * X, using LinearAlgebra.svd.

`X` is internally centered.
"""
function pcasvd(X, weights = ones(size(X, 1)); nlv)
    pcasvd!(copy(X), weights; nlv = nlv)
end

function pcasvd!(X, weights = ones(size(X, 1)); nlv)
    X = ensure_mat(X)
    n, p = size(X)
    nlv = min(nlv, n, p)
    weights = mweights(weights)
    sqrtw = sqrt.(weights)
    xmeans = colmeans(X, weights)
    center!(X, xmeans)
    res = LinearAlgebra.svd!(Diagonal(sqrtw) * X)
    P = res.V[:, 1:nlv]
    ## by default in LinearAlgebra.svd
    ## "full = false" ==> [1:min(n, p)]
    sv = res.S
    sv[sv .< 0] .= 0
    ## eig = sv.^2
    ## = Eigenvalues of X'DX = Cov(X) in metric D
    ## where D = Diagonal(weights)
    ## = variances of scores T in metric D
    ## = TT = colSums(weights * T .* T)
    ## = norms^2 of the scores T in metric D
    ## = colnorms(T, weights = weights)^2
    ## -------- Strictly, T should be generated by transform(object, X)
    ## But T is returned below to keep consistency with PLS,
    ## where T is a direct output of the algorithm
    ## T = D^(-1/2) * U * S
    ## Faster but requires weights > 0 (or mananage NA vales for weights = 0)
    ## Alternative: T = X * P
    T = Diagonal(1 ./ sqrtw) * vcol(res.U, 1:nlv) * (Diagonal(vrow(sv, 1:nlv)))
    Pca(T, P, sv, xmeans, weights, nothing, nothing)
end

"""
    summary(object::Pca, X)
Summarize the fitted model.
* `object` : The fitted model.
* `X` : The X-data that was used to fit the model.
"""
function Base.summary(object::Pca, X)
    nlv = size(object.T, 2)
    X = center(X, object.xmeans)
    sstot = sum(object.weights' * (X.^2))
    D = Diagonal(object.weights)
    TT = D * object.T.^2
    tt = vec(sum(TT, dims = 1))
    pvar = tt / sstot
    cumpvar = cumsum(pvar)
    explvar = DataFrame(pc = 1:nlv, var = tt, pvar = pvar, cumpvar = cumpvar)
    contr_ind = DataFrame(scale(TT, tt), :auto)
    xvars = colvars(X, object.weights)
    zX = scale(X, sqrt.(xvars))
    zT = D * scale(object.T, sqrt.(tt))
    cor_circle = DataFrame(zX' * zT, :auto)
    z = X' * zT
    coord_var = DataFrame(z, :auto)
    z = z .* z
    contr_var = DataFrame(scale(z, sum(z, dims = 1)), :auto)
    nam = map(string, repeat(["pc"], nlv), 1:nlv)
    rename!(contr_ind, nam)
    rename!(cor_circle, nam)
    rename!(coord_var, nam)
    rename!(contr_var, nam)
    (explvar = explvar, contr_ind = contr_ind,
        contr_var = contr_var, coord_var = coord_var, cor_circle = cor_circle)
end

"""
    transform(object::Pca, X; nlv = nothing)
Compute components (scores matrix "T") from a fitted model and X-data.
* `object` : The maximal fitted model.
* `X` : X-data for which PCs are computed.
* `nlv` : Nb. components to compute. If nothing, it is the maximum nb. PCs.
"""
function transform(object::Union{Pca, Fda}, X; nlv = nothing)
    a = size(object.T, 2)
    isnothing(nlv) ? nlv = a : nlv = min(nlv, a)
    center(X, object.xmeans) * vcol(object.P, 1:nlv)
end


text:
using HPAT
using MPI
using DocOpt
ParallelAccelerator.ParallelIR.PIRHoistParfors(1)
ParallelAccelerator.CGen.setRecreateLoops(true)
ParallelAccelerator.CGen.setRecreateConds(true)
HPAT.DistributedPass.DistOptimize(true)
#HPAT.DistributedPass.set_debug_level(3)

@acc hpat function linear_regression(iterations, N)
    D = 10  # Number of features
    p = 4
    labels = rand(p,N)
    points = rand(D,N)
    w = zeros(p,D)
    alphaN = 0.01/N

    for i in 1:iterations
       w -= alphaN*((w*points)-labels)*points'
    end
    w
end

function main()
    doc = """linear regression statistical method.

Usage:
  linear_regression.jl -h | --help
  linear_regression.jl [--iterations=<iterations>] [--instances=<instances>]

Options:
  -h --help                  Show this screen.
  --iterations=<iterations>  Specify number of iterations; defaults to 20.
  --instances=<instances>    Specify number of instances; defaults to 10^7.
"""
    arguments = docopt(doc)
    iterations = 20
    instances = 10^7

    if (arguments["--iterations"] != nothing)
        iterations = parse(Int, arguments["--iterations"])
    end

    if (arguments["--instances"] != nothing)
        instances = parse(Int, arguments["--instances"])
    end

    W = linear_regression(iterations, instances)

    if MPI.Comm_rank(MPI.COMM_WORLD)==0 println("result = ", W) end
end

main()

max_stars_repo_name: ehsantn/data-analytics-benchmarks

max_issues_repo_name: ehsantn/data-analytics-benchmarks

max_forks_repo_name: ehsantn/data-analytics-benchmarks