text:
###########################################################################################################
###########################################################################################################
##################### OVERLAP LIKE G VİA SERRIES REPRESENTATION OF BETA FUNCTOINS #########################
###########################################################################################################
###########################################################################################################
export AuxiliaryG
###########################################################################################################
function AuxiliaryG(n1::arb, q::Int, n2::arb, n3::arb, p1::arb, p2::arb, p3::arb, lim::Int)

    zero = parent(n1)(0)
    one = parent(n1)(1)
    two = parent(n1)(2)
    three = parent(n1)(3)
    four = parent(n1)(4)

    lowq = lowlim = lowlimq = low2limq = low2lim2q = low3lim2q = low3lim3q = low4lim3q = low4lim4q =
    low5lim4q = low5lim5q = 0

    upq = q
    uplim = lim
    uplimq = lim + q
    up2limq = 2*lim + q
    up2lim2q = 2*lim + 2*q
    up3lim2q = 3*lim + 2*q
    up3lim3q = 3*lim + 3*q
    up4lim3q = 4*lim + 3*q
    up4lim4q = 4*lim + 4*q
    up5lim4q = 5*lim + 4*q
    up5lim5q = 5*lim + 5*q
    ###################################################
    ###################################################
    m1c1 = zeros(RF, uplim - lowlim + 1, up5lim4q - low5lim4q + 1)
    m1c2 = zeros(RF, uplim - lowlim + 1, up5lim4q - low5lim4q + 1)
    m2c1 = zeros(RF, uplim - lowlim + 1, up5lim4q - low5lim4q + 1)
    m2c2 = zeros(RF, uplim - lowlim + 1, up5lim4q - low5lim4q + 1)

    m1r1 = zeros(RF, uplim - lowlim + 1, 2)
    m1r2 = zeros(RF, uplim - lowlim + 1, 2)
    m2r1 = zeros(RF, uplim - lowlim + 1, 2)
    m2r2 = zeros(RF, uplim - lowlim + 1, 2)

    m1 = zeros(RF, uplim - lowlim + 1, up5lim4q - low5lim4q + 1)
    m2 = zeros(RF, uplim - lowlim + 1, up5lim4q - low5lim4q + 1)
    #################
    #################
    p13rc = zeros(RF, uplim - lowlim + 1, 2)
    p23rc = zeros(RF, uplim - lowlim + 1, 2)

    p13cc = zeros(RF, uplim - lowlim + 1, up5lim4q - low5lim4q + 1)
    p23cc = zeros(RF, uplim - lowlim + 1, up5lim4q - low5lim4q + 1)

    p13 = zeros(RF, uplim - lowlim + 1, up5lim4q - low5lim4q + 1)
    p23 = zeros(RF, uplim - lowlim + 1, up5lim4q - low5lim4q + 1)
    #################
    #################
    s4c1 = zeros(RF, uplim - lowlim + 1, up5lim4q - low5lim4q + 1)
    s4c2 = zeros(RF, uplim - lowlim + 1, up5lim4q - low5lim4q + 1)
    #################
    #################
    s4f1 = s4f2 = zeros(RF, up5lim4q - low5lim4q + 1)
    #fs2 is t he the gamma function Gamma(S_2 + 1)#
    fs2 = zeros(RF, uplim - lowlim + 1)
    #################
    #################
    p12 = zeros(RF, up2lim2q - low2lim2q + 1, up3lim2q - low3lim2q + 1)
    p22 = zeros(RF, up2lim2q - low2lim2q + 1, up3lim2q - low3lim2q + 1)
    #################
    #################
    p11 = zeros(RF, uplim - lowlim + 1, uplimq - lowlimq + 1)
    p21 = zeros(RF, uplim - lowlim + 1, uplimq - lowlimq + 1)
    #################
    #################
    b1rc = zeros(RF, uplim - lowlim + 1, 2)
    b2rc = zeros(RF, uplim - lowlim + 1, 2)

    b1cc = zeros(RF, uplim - lowlim + 1, uplimq - lowlimq + 1)
    b2cc = zeros(RF, uplim - lowlim + 1, uplimq - lowlimq + 1)

    B = zeros(RF, uplim - lowlim + 1, uplimq - lowlimq + 1)
    B = zeros(RF, uplim - lowlim + 1, uplimq - lowlimq + 1)
    #################
    #################
    Bins32 = zeros(RF,  uplim - lowlim + 1, uplim - lowlim + 1)
    #################
    #################
    E = zeros(RF, uplim - lowlim + 1)
    ###################################################
    ###################################################
    m1[1,1] = real(AuxiliaryGm(n2 + 2*q, n3 + one, p2))
    m1[2,1] = real(AuxiliaryGm(n2 + 2*q + one, n3 + one, p2))
    m1[1,2] = real(AuxiliaryGm(n2 + 2*q - one, n3 + two, p2))
    m1[2,2] = real(AuxiliaryGm(n2 + 2*q, n3 + two, p2))

    m2[1,1] = real(AuxiliaryGm(n3 + 2*q, n2 + one, p2))
    m2[2,1] = real(AuxiliaryGm(n3 + 2*q + one, n2 + one, p2))
    m2[1,2] = real(AuxiliaryGm(n3 + 2*q - one, n2 + two, p2))
    m2[2,2] = real(AuxiliaryGm(n3 + 2*q, n2 + two, p2))
    #################
    #################
    p13[1,1] = Pochhammer(-(n2 + 2*q), 0)
    p13[2,1] = Pochhammer(-(n2 + 2*q + two), one)
    p13[1,2] = Pochhammer(-(n2 + 2*q) + one, 0)
    p13[2,2] = Pochhammer(-(n2 + 2*q + two) + one, one)

    p23[1,1] = Pochhammer(-(n3 + 2*q), 0)
    p23[2,1] = Pochhammer(-(n3 + 2*q + two), one)
    p23[1,2] = Pochhammer(-(n3 + 2*q) + one, 0)
    p23[2,2] = Pochhammer(-(n3 + 2*q + two) + one, one)
    #################
    #################
    s4c1[1,1] = s4c1[2,1] = (n3 + one)
    s4c1[1,2] = s4c1[2,2] = (n3 + two)

    s4c2[1,1] = s4c2[2,1] = (n2 + one)
    s4c2[1,2] = s4c2[2,2] = (n2 + two)
    #################
    #################
    s4f1[1] = s4f2[1] = fs2[1] = RF(1)
    s4f1[2] = s4f2[2] = fs2[2] = RF(1)
    #################
    #################
    z1 = n2 + 2*q + one
    z2 = n3 + one
    B[1,1] = Beta(z1, z2)
    B[1,2] = Beta(z1 - two, z2 + two)

    B[2,1] = Beta(z1 + two, z2)
    B[2,2] = Beta(z1, z2 + two)
    #################
    #################
    p11[1,1] = Pochhammer(-(n2 + 2*q), 0)
    p11[2,1] = Pochhammer(-(n2 + 2*q + two), one)
    p11[1,2] = Pochhammer(-(n2 + 2*q - two), 0)
    p11[2,2] = Pochhammer(-(n2 + 2*q), one)

    p21[1,1] = Pochhammer(-(n3 + 2*q), 0)
    p21[2,1] = Pochhammer(-(n3 + 2*q + two), one)
    p21[1,2] = Pochhammer(-(n3 + 2*q - two), 0)
    p21[2,2] = Pochhammer(-(n3 + 2*q), one)
    #################
    #################
    Bins32[1,1] = RF(1)
    Bins32[2,1] = RF(1)
    Bins32[1,2] = zero
    Bins32[2,2] = RF(1)
    #################
    #################
    α = (n2 + n3 + 2*q + one)
    E[1] =  ExpIntegralE(-α, p2)
    E[2] =  ExpIntegralE(-(α + one), p2)
    ###################################################
    ###################################################
    for s1 in lowlim : uplim - 2
        s11r = s1 - lowlim + 3

        m1r1[s11r,1] = (-(n2 + s11r - 1) - 2*q + one) // (p2)
        m1r2[s11r,1] = ((n2 + s11r - 1) + n3 + 2*q + p2 + one) // (p2)
        m1r1[s11r,2] = (one - (n2 + s11r - 1) - 2*q + one) // (p2)
        m1r2[s11r,2] = ((n2 + s11r - 1) + n3 + 2*q + p2 + one) // (p2)

        m2r1[s11r,1] = (-(n3 + s11r - 1) - 2*q + one) // (p2)
        m2r2[s11r,1] = ((n3 + s11r - 1) + n2 + 2*q + p2 + one) // (p2)
        m2r1[s11r,2] = (one - (n3 + s11r - 1) - 2*q + one) // (p2)
        m2r2[s11r,2] = ((n3 + s11r - 1) + n2 + 2*q + p2 + one) // (p2)

        m1[s11r,1] = (40//10) * m1r1[s11r,1] * m1[s11r - 2,1] + (20//10) * m1r2[s11r,1] * m1[s11r - 1,1]
        m1[s11r,2] = (40//10) * m1r1[s11r,2] * m1[s11r - 2,2] + (20//10) * m1r2[s11r,2] * m1[s11r - 1,2]

        m2[s11r,1] = (40//10) * m2r1[s11r,1] * m2[s11r - 2,1] + (20//10) * m2r2[s11r,1] * m2[s11r - 1,1]
        m2[s11r,2] = (40//10) * m2r1[s11r,2] * m2[s11r - 2,2] + (20//10) * m2r2[s11r,2] * m2[s11r - 1,2]

        p13rc[s11r,1] = (
            (
                (-(n2 + 2*q) - 2*(s11r - 1))
                //
                (-(n2 + 2*q) - (s11r - 1))
            ) * (-(n2 + 2*q) - 2*(s11r - 1) + 1)
        )
        p23rc[s11r,1] = (
            (
                (-(n3 + 2*q) - 2*(s11r - 1))
                //
                (-(n3 + 2*q) - (s11r - 1))
            ) * (-(n3 + 2*q) - 2*(s11r - 1) + 1)
        )
        p13rc[s11r,2] = (
            (
                (-(n2 - one + 2*q) - 2*(s11r - 1))
                //
                (-(n2 - one + 2*q) - (s11r - 1))
            ) * (-(n2 - one + 2*q) - 2*(s11r - 1) + 1)
        )
        p23rc[s11r,2] = (
            (
                (-(n3 - one + 2*q) - 2*(s11r - 1))
                //
                (-(n3 - one + 2*q) - (s11r - 1))
            ) * (-(n3 - one + 2*q) - 2*(s11r - 1) + 1)
        )

        p13[s11r,1] = p13rc[s11r,1] * p13[s11r - 1,1]
        p13[s11r,2] = p13rc[s11r,2] * p13[s11r - 1,2]

        p23[s11r,1] = p23rc[s11r,1] * p23[s11r - 1,1]
        p23[s11r,2] = p23rc[s11r,2] * p23[s11r - 1,2]

        s4c1[s11r,1] = s4c1[1,1]
        s4c1[s11r,2] = s4c1[1,2]

        s4c2[s11r,1] = s4c2[1,1]
        s4c2[s11r,2] = s4c2[1,2]

        fs2[s11r] = (s11r - 1) * fs2[s11r - 1]

        b1rc[s11r,1] = (
            ((z1 + 2*(s11r - 1) - 1) * (z1 + 2*(s11r - 1) - 2))
            //
            ((z1 + z2 + 2*(s11r - 1) - 1) * (z1 + z2 + 2*(s11r - 1) - 2))
        )

        b2rc[s11r,1] = (
            ((z1 + 2*(s11r - 1) - 3) * (z1 + 2*(s11r - 1) - 4))
            //
            ((z1 + z2 + 2*(s11r - 1) - 1) * (z1 + z2 + 2*(s11r - 1) - 2))
        )

        B[s11r,1] =  b1rc[s11r,1] *  B[s11r - 1,1]
        B[s11r,2] =  b2rc[s11r,1] *  B[s11r - 1,2]

        Bins32[s11r,1] = RF(1)
        Bins32[s11r,2] = RF(s11r - 1)

        E[s11r] = (one // p2) * (Exp(-p2) +  (α + s11r - one) * E[s11r - 1])
    end
    #################
    #################
    for s1 in lowlim : uplim
        s11c = s1 - lowlim + 1

        for s2 in low5lim4q : lim + 4*q + 4*s1 - 2
            s21 = s2 - low5lim4q + 3

            m1c1[s11c,s21] = (n3+(s21-1))//((n2+s1)+2*q+one-(s21-1))
            m1c2[s11c,s21] = ((n2+s1)-n3-2*(s21-1)+2*q+one-p2)//((n2+s1)+2*q+one-(s21-1))

            m2c1[s11c,s21] = (n2+(s21-1))//((n3+s1)+2*q+one-(s21-1))
            m2c2[s11c,s21] = ((n3+s1)-n2-2*(s21-1)+2*q+one-p2)//((n3+s1)+2*q+one-(s21-1))

            m1[s11c,s21] = (10//40)*m1c1[s11c,s21]*m1[s11c,s21-2]+(10//20)*m1c2[s11c,s21]*m1[s11c,s21-1]
            m2[s11c,s21] = (10//40)*m2c1[s11c,s21]*m2[s11c,s21-2]+(10//20)*m2c2[s11c,s21]*m2[s11c,s21-1]

            p13cc[s11c,s21] = ((-(n2 + 2*q) - s1 + (s21 - 1) - one) // (-(n2 + 2*q)-2*s1 + (s21 - 1) - 1))
            p23cc[s11c,s21] = ((-(n3 + 2*q) - s1 + (s21 - 1) - one) // (-(n3 + 2*q)-2*s1 + (s21 - 1) - 1))

            p13[s11c,s21] =  p13cc[s11c,s21] * p13[s11c,s21 - 1]
            p23[s11c,s21] =  p23cc[s11c,s21] * p23[s11c,s21 - 1]

            s4c1[s11c,s21] = (n3 + s2 + two + one)
            s4c2[s11c,s21] = (n2 + s2 + two + one)

            if s21 <= lim + 2*q + 2*s11c - 1
                s4f1[s21] = s4f2[s21] = (s21 - 1) * s4f1[s21 - 1]
            end
            #################
            if s11c > 2 && s21 <= s11c
                Bins32[s11c,s21] = Bins32[s11c - 1,s21 - 1] + Bins32[s11c - 1,s21]
            end
        end
    end
    #################
    #################
    for s1 in low2lim2q : up2lim2q
        s11 = s1 - low2lim2q + 1

        p12[s11,1] = one
        p12[s11,2] = -(n2 + 2*q + lim - s11 + 1)

        p22[s11,1] = one
        p22[s11,2] = -(n3 + 2*q + lim - s11 + 1)

        for s2 in low3lim2q : up3lim2q - 2
            s21 = s2 - low3lim2q + 3

            p12[s11,s21] = (-(n2 + 2*q + lim - s11 + 1) + s21 - 1 - 1) * p12[s11,s21 - 1]

            p22[s11,s21] = (-(n3 + 2*q + lim - s11 + 1) + s21 - 1 - 1) * p22[s11,s21 - 1]

        end
    end
    #################
    #################
    for s1 in lowlim : uplim
        s11 = s1 - lowlim + 1

        p11[s11,1] = p13[s11,1]
        p21[s11,1] = p23[s11,1]

        p11[s11,2] = p13[s11,3]
        p21[s11,2] = p23[s11,3]

        for s2 in lowq : upq + s1 - 2
            s21 = s2 - lowq + 3

            p11[s11,s21] = p13[s11,2*s21-1]
            p21[s11,s21] = p23[s11,2*s21-1]

            B[s11,s21] = (
                ((z2 + 2*s21 - 3) * (z2 + 2*s21 - 4))
                //
                ((z1 + 2*s1 - 2*(s21 - 1)) * (z1 + 2*s1 - 2*(s21 - 1) + 1))
            ) *  B[s11,s21 - 1]
        end
    end
    #################
    #################
    LES1 = zeros(RF, upq - lowq + 1, uplim - lowlim + 1, uplim - lowlim + 1, up3lim2q - low3lim2q + 1)
    LES2 = zeros(RF, upq - lowq + 1, uplim - lowlim + 1, uplim - lowlim + 1)
    BES = zeros(RF, upq - lowq + 1, uplim - lowlim + 1, uplim - lowlim + 1)
    JES1 = zeros(RF, upq - lowq + 1, uplim - lowlim + 1, uplim - lowlim + 1)
    JES2 = zeros(RF, upq - lowq + 1, uplim - lowlim + 1)
    GES1 = zeros(RF, upq - lowq + 1)
    GES2 = RF(0)
    GES3 = RF(0)
    #################
    #################
    for s1 in lowq : upq
        for s2 in lowlim : uplim
            for s3 in lowlim : s2

                BES[s1 + 1,s2 + 1,s3 + 1] =
                ((p3^s2)  // fs2[s2 + 1]) *
                (
                B[s2 + 1,s1 + s3 + 1] *
                E[s2 + 1]
                )

                for s4 in low3lim2q : lim + 2*q + 2*s2

                    LES1[s1 + 1,s2 + 1,s3 + 1,s4 + 1] =
                    ((p3^s2)  // fs2[s2 + 1]) *
                    (
                    (   p11[s2 + 1, s1 + s3 + 1] *
                        m1[s2 + 1,2*s1 + 2*s3 + s4 + 1]
                        //
                        (
                            p13[s2 + 1,2*s1 + 2*s3 + s4 + 1]
                            *
                            s4c1[s2 + 1,2*s1 + 2*s3 + s4 + 1]
                            *
                            s4f1[s4 + 1]
                        ) *
                        p12[lim + 2*s1 - s2 + 2*s3 + 1,s4 + 1]
                    ) + #-
                    (   p21[s2 + 1,(q - s1) + (s2 - s3) + 1]  *
                        m2[s2 + 1,2*(q - s1) + 2*(s2 - s3) + s4 + 1]
                        //
                       (
                            p23[s2 + 1,2*(q - s1) + 2*(s2 - s3) + s4 + 1]
                            *
                            s4c2[s2 + 1,2*(q - s1) + 2*(s2 - s3) + s4 + 1]
                            *
                            s4f2[s4 + 1]
                       ) *
                       p22[lim + 2*(q - s1) - s2 + 2*(s2 - s3) + 1,s4 + 1]
                    )
                    )

                    LES2[s1 + 1,s2 + 1,s3 + 1] += LES1[s1 + 1,s2 + 1,s3 + 1,s4 + 1]

                end

                JES1[s1 + 1,s2 + 1,s3 + 1] =
                Power((20//10), (s2)) * (Power((20//10),(n2 + n3 + 2*q + s2 + 1)) * BES[s1 + 1,s2 + 1,s3 + 1]) -
                Power((20//10),(s2)) * LES2[s1 + 1,s2 + 1,s3 + 1]

                JES2[s1 + 1,s2 + 1] +=
                ((10//10) // Power(20//10, 2*s2)) * Power((-10//10),(s3)) * Bins32[s2 + 1, s3 + 1] *
                JES1[s1 + 1,s2 + 1,s3 + 1]

            end

            GES1[s1 + 1] += Power((-10//10),(s1)) * Bins32[q + 1,s1 + 1] *
            Power(-10//10, s2) * JES2[s1 + 1,s2 + 1]

        end
        GES2 += GES1[s1 + 1]
    end

    GES3 = ((10//10) // ((20//10)^(2*q))) * GES2
    Res = ((p1^(n1))//Gamma(n1+1)) * GES3

    return Res
end


text:
include("plot_utils.jl")

θ = 180.0
ϕ = 0.0

#title = "boost_regular_12"

fname = "../data/voronoi_ul7n12_2e6.h5"
title = split(fname, "/")[end]
title = split(title, ".")[1]

plot_convergence(fname, string(title))
write_convergence(fname, string(title))

atmos, line, S_λ, α_tot = plotter(read_irregular(fname)..., θ, ϕ, fname)
# atmos, line, S_λ, α_tot = plotter(read_quantities(fname, periodic=true)..., θ, ϕ, fname)
write_top_intensity(atmos, line, S_λ, α_tot, θ, ϕ, "$(title)_disk_centre_1")

# plot_top_line(atmos, line, S_λ, α_tot, θ, ϕ, title)

# for idλ in line.λidx[1]+1 : line.λidx[2]
#     plot_top_intensity(atmos, line, S_λ, α_tot, θ, ϕ, idλ, title*string(idλ))
# end

# plot_convergence("../data/voronoi_ul2n3_2.h5", "Irregular grid convergence")
# plot_convergence("../data/regular_ul2n3_zero_radiation_1.h5", "Regular grid convergence")

ext: jl

max_stars_repo_name: meudnaes/VoronoiRT

max_issues_repo_name: meudnaes/VoronoiRT

max_forks_repo_name: meudnaes/VoronoiRT


text:
module ITensorsGateEvolution

using ITensors

export ProductOps,
       ops,
       apply,
       findsiteinds,
       firstsiteinds,
       movesite,
       movesites,
       get_warn_itensor_order,
       set_warn_itensor_order!,
       reset_warn_itensor_order!,
       disable_warn_itensor_order!

# Extensions to ITensors
include("itensors/global_variables.jl")
include("itensors/tupletools.jl")
include("itensors/index.jl")
include("itensors/indexset.jl")
include("itensors/itensor.jl")
include("itensors/mps.jl")
include("itensors/sitetype.jl")

# ProductOps type
include("productops.jl")

# Qubit site type
include("qubit.jl")

# Apply function
include("apply.jl")

end


text:
include("../src/ZZDiffusionBridge.jl")
using Distributions
using Plots

struct SinSDE <: AbstractModel
    α::Float64
    # optional precompiled factors
end

"""

`ss` is the stepsize, `n` batchsize for Monte Carlo integration, `nstep' sampler
"""
struct NewSystem
    ξ::Array{Float64, 1}
    ϕ::Vector{Fs}
    L::Int64
    T::Float64
    drift::Array{Float64, 1}
    function NewSystem(L, T)
        new(fill(0.0, 2<<L-1), generate(L, T), L, T, fill(0.0, 2<<L-1))
    end
end

function stoch_langevin(S::NewSystem, X::AbstractModel, ss, nstep, skip, u, v, n)
    d = 2<<S.L - 1
    if ss >= 1
        error("stepsize has to be strictly less the 1")
    end
    Ξ = []
    for i in 1:nstep
        stoch_drift!(S, X, u, v, n)
        S.ξ .=  S.ξ - S.drift*0.5*ss .+ sqrt(ss)*randn(d)  #otherwise: setfield! immutable struct of type NewSystem cannot be changed

        if i%skip == 0
            push!(Ξ, copy(S.ξ))
        end
    end
    Ξ
end

function stoch_drift!(S::NewSystem, X::AbstractModel, u, v, nn)
    for n in 1:2<<S.L-1
        ave = 0.0
        w = 1/nn
        for i in 1:nn
            t = MCintegration(S.ϕ[n])
            XX = fs_expansion(S, t, u, v)
            λ = Λ(t, S.ϕ[n].i , S.ϕ[n].j, S.T)
            ave += w*(0.5*S.ϕ[n].range*λ*(X.α*X.α*sin(2.0*(XX)) - X.α*sin(XX)) + S.ξ[n])
        end
        S.drift[n] = ave
    end
end

function fs_expansion(S::NewSystem, t::Float64, u::Float64, v::Float64, n = i -> 2^-(1 + i/2))
        dt = 0:S.T/(2<<S.L):S.T
        k = (searchsortedfirst(dt, t) - 1)
        j0 =  Int(ceil(k/2))-1
        n0 = Faber(S.L, j0)
        if k % 2 != 0
                return interpolate([dt[k], dt[k + 1]], fs_expansion(S.ϕ[n0], S.ξ, u, v, S.L, S.T)[1:2], t)
        else
                return interpolate([dt[k], dt[k + 1]], fs_expansion(S.ϕ[n0], S.ξ, u, v,  S.L, S.T)[2:3], t)
        end
end

function runall(SHORT = false)
    L = 6
    T = 100.0
    S = NewSystem(L, T)
    α = 0.7
    X = SinSDE(α)
    step_size = 0.002
    nstep = 2*10^5
    skip = 4*10^3
    u = -3.0
    v = 3.0
    n = 1
    Y = stoch_langevin(S, X, step_size, nstep, skip, u, v, n)
    if SHORT == false
        burning  = 10
        p = plot(leg = false)
        dt = 0:T/2<<L:T
        for i in burning:length(Y)
            plot!(p, dt, fs_expansion(Y[i], u, v, L, T),  color= RGB(i/length(Y), 0.4, 1-i/length(Y)), linewidth=0.3, alpha = 0.4)
        end
        hline!(p, [n*π for n in -7:2:7])
        display(p)
    end
    return Y, p
end

Y, p = runall()

error("STOP HERE")
png("output/langevin_sampler_sin_07.png")

max_stars_repo_name: SebaGraz/ZZDiffusionBridge

max_issues_repo_name: SebaGraz/ZZDiffusionBridge

max_forks_repo_name: SebaGraz/ZZDiffusionBridge


text:
using Plots
using LaTeXStrings
using Formatting
using ComplexMixtures, PDBTools
const CM = ComplexMixtures

# Default plot options (too look good with LaTeXStrings)
plot_font = "Computer Modern"
default(fontfamily=plot_font,
        linewidth=2, framestyle=:box, label=nothing)

# PDB file of the system simulated
pdb = readPDB("../Data/system.pdb")

# Load results of a ComplexMixtures run
R = CM.load("../Data/results_glyc50.json")

# Inform which is the solute
protein = select(pdb,"protein")
solute = CM.Selection(protein,nmols=1)

# Inform which is the solvent
glycerol = select(pdb,"resname GLYC")
solvent = CM.Selection(glycerol,natomspermol=14)

# Plot a 2D map showing the contributions of some residues
residues = collect(eachresidue(protein))
irange=70:110
rescontrib = zeros(length(R.mddf),length(residues))
for (ires,residue) in pairs(residues)
  rescontrib[:,ires] .= CM.contrib(solute,R.solute_atom,residue)
end

# Plot only for distances within 1.5 and 3.5:
idmin = findfirst( d -> d > 1.5, R.d)
idmax = findfirst( d -> d > 3.5, R.d)

# Obtain pretty labels for the residues in the x-axis
labels = PDBTools.oneletter.(resname.(residues)).*format.(resnum.(residues))

contourf(irange,R.d[idmin:idmax],rescontrib[idmin:idmax,irange],
         color=cgrad(:tempo),linewidth=0.1,linecolor=:black,
         colorbar=:none,levels=5,
         xlabel="Residue",ylabel=L"r/\AA",
         xticks=(irange,labels[irange]),xrotation=60,
         xtickfont=font(6,plot_font),
         size=(500,280))

savefig("./density2D.pdf")


text:
path = @get_scratch!("rocket")

# rocket
nx = rocket.nq
nu = rocket.nu
nw = rocket.nw

nz = nx
nθ = nx + nu + 1

@variables z[1:nz] θ[1:nθ] κ[1:1]

# Residual
function residual(rocket, z, θ, κ)
    y = z[1:nx]
    x = θ[1:nx]
    u = θ[nx .+ (1:nu)]
    h = θ[nx + nu .+ (1:1)]
    w = zeros(rocket.nw)

    return y - (x + h[1] * f(rocket, 0.5 * (x + y), u, w))
end

r_rocket = residual(rocket, z, θ, κ)
r_rocket = Symbolics.simplify.(r_rocket)
rz_rocket = Symbolics.jacobian(r_rocket, z, simplify=true)
rθ_rocket = Symbolics.jacobian(r_rocket, θ, simplify=true)

# Build function
r_rocket_func = build_function(r_rocket, z, θ, κ)[2]
rz_rocket_func = build_function(rz_rocket, z, θ)[2]
rθ_rocket_func = build_function(rθ_rocket, z, θ)[2]

rz_rocket_array = similar(rz_rocket, Float64)
rθ_rocket_array = similar(rθ_rocket, Float64)

@save joinpath(path, "residual.jld2") r_rocket_func rz_rocket_func rθ_rocket_func rz_rocket_array rθ_rocket_array
@load joinpath(path, "residual.jld2") r_rocket_func rz_rocket_func rθ_rocket_func rz_rocket_array rθ_rocket_array

# projection
nz = 3 + 1 + 1 + 1 + 1 + 3
nθ = 3 + 1

function residual_projection(z, θ, κ)
    u = z[1:3]
    p = z[4:4]
    s = z[5:5]
    w = z[6:6]
    y = z[7:7]
    v = z[8:10]

    ū = θ[1:3]
    uu = θ[3 .+ (1:1)]
    idx = [3; 1; 2]
    [
     u - ū - v - [0.0; 0.0; y[1] + p[1]];
     uu[1] - u[3] - s[1];
     -y[1] - w[1];
     w[1] * s[1] - κ[1];
     p[1] * u[3] - κ[1];
     cone_product(u[idx], v[idx]) - [κ[1]; 0.0; 0.0]
    ]
end

@variables z[1:nz], θ[1:nθ], κ[1:1]

r_proj = residual_projection(z, θ, κ)
r_proj .= simplify(r_proj)
r_proj_func = Symbolics.build_function(r_proj, z, θ, κ)[2]

rz_proj = Symbolics.jacobian(r_proj, z)
rz_proj = simplify.(rz_proj)
rz_proj_func = Symbolics.build_function(rz_proj, z, θ)[2]
rθ_proj = Symbolics.jacobian(r_proj, θ)
rθ_proj = simplify.(rθ_proj)
rθ_proj_func = Symbolics.build_function(rθ_proj, z, θ)[2]

rz_proj_array = similar(rz_proj, Float64)
rθ_proj_array = similar(rθ_proj, Float64)

@save joinpath(path, "projection.jld2") r_proj_func rz_proj_func rθ_proj_func rz_proj_array rθ_proj_array
@load joinpath(path, "projection.jld2") r_proj_func rz_proj_func rθ_proj_func rz_proj_array rθ_proj_array


text:
driver = "470.57.02"

sources_linux_x86 = [
    FileSource("https://us.download.nvidia.com/tesla/$driver/NVIDIA-Linux-x86_64-$driver.run",
               "55d7ae104827faa79e975321fe2b60f9dd42fbff65642053443c0e56fdb4c47d", "installer.run")
]
sources_linux_ppc64le = [
    FileSource("https://us.download.nvidia.com/tesla/$driver/NVIDIA-Linux-ppc64le-$driver.run",
               "f430c6f898e0e2e1530e75f3770a29fd6b4e5628fea0b305fdd501a1a951b067", "installer.run")
]
sources_linux_aarch64 = [
    FileSource("https://us.download.nvidia.com/tesla/$driver/NVIDIA-Linux-aarch64-$driver.run",
               "4b3e51963b88e8de304a6bca68c64954844ed490156a28483e9272eb79a00451", "installer.run")
]

if should_build_platform("x86_64-linux-gnu-cuda+11.4")
    build_tarballs(non_reg_ARGS, name, version, sources_linux_x86, script,
                   [Platform("x86_64", "linux"; cuda="11.4")], products, dependencies;
                   lazy_artifacts=true, skip_audit=true)
end

if should_build_platform("powerpc64le-linux-gnu-cuda+11.4")
    build_tarballs(non_reg_ARGS, name, version, sources_linux_ppc64le, script,
                   [Platform("powerpc64le", "linux"; cuda="11.4")], products, dependencies;
                   lazy_artifacts=true, skip_audit=true)
end

if should_build_platform("aarch64-linux-gnu-cuda+11.4")
    build_tarballs(ARGS, name, version, sources_linux_aarch64, script,
                   [Platform("aarch64", "linux"; cuda="11.4")], products, dependencies;
                   lazy_artifacts=true, skip_audit=true)
end

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
function Base.showarg(io::IO, A::BSplineInterpolation{T,N,TW,ST}, toplevel) where {T,N,TW,ST}
    print(io, "interpolate(")
    Base.showarg(io, A.coefs, false)
    print(io, ", ")
    show(io, itpflag(A))
    if toplevel
        print(io, ") with element type ",T)
    else
        print(io, ')')
    end
end

function Base.showarg(io::IO, A::GriddedInterpolation{T,N,TC,ST,K}, toplevel) where {T,N,TC,ST,K}
    print(io, "interpolate(")
    _showknots(io, A.knots)
    print(io, ", ")
    Base.showarg(io, A.coefs, false)
    print(io, ", ")
    show(io, itpflag(A))
    if toplevel
        print(io, ") with element type ",T)
    else
        print(io, ')')
    end
end

_showknots(io, A) = Base.showarg(io, A, false)
function _showknots(io, tup::NTuple{N,Any}) where N
    print(io, '(')
    for (i, A) in enumerate(tup)
        Base.showarg(io, A, false)
        i < N && print(io, ',')
    end
    N == 1 && print(io, ',')
    print(io, ')')
end

function Base.showarg(io::IO, A::ScaledInterpolation{T}, toplevel) where {T}
    print(io, "scale(")
    Base.showarg(io, A.itp, false)
    print(io, ", ", A.ranges, ')')
    if toplevel
        print(io, " with element type ",T)
    end
end

function Base.showarg(io::IO, A::MonotonicInterpolation{T, TCoeffs, Tel, Type, K}, toplevel) where {T, TCoeffs, Tel, Type, K}
    print(io, "interpolate(")
    _showknots(io, A.knots)
    print(io, ", ")
    Base.showarg(io, A.A, false)
    print(io, ", ")
    show(io, A.it)
    if toplevel
        print(io, ") with element type ",T)
    else
        print(io, ')')
    end
end

function Base.showarg(io::IO, A::Extrapolation{T,N,TI,IT,ET}, toplevel) where {T,N,TI,IT,ET}
    print(io, "extrapolate(")
    Base.showarg(io, A.itp, false)
    print(io, ", ")
    show(io, etpflag(A))
    print(io, ')')
    if toplevel
        print(io, " with element type ",T)
    end
end

function Base.showarg(io::IO, A::FilledExtrapolation{T,N,TI,IT}, toplevel) where {T,N,TI,IT}
    print(io, "extrapolate(")
    Base.showarg(io, A.itp, false)
    print(io, ", ", A.fillvalue, ')')
    if toplevel
        print(io, " with element type ",T)
    end
end


text:
function create_grid(mesh::RawMesh, format, user_elements)
    dim = getdim(mesh)
    cells = create_cells(getelementsdicts(mesh), user_elements, format)
    nodes = create_nodes(getnodes(mesh), Val(dim))
    cellsets = create_cellsets(getelementsdicts(mesh), getelementsets(mesh))
    nodesets = create_nodesets(getnodesets(mesh))
    return Grid(cells, nodes; cellsets=cellsets, nodesets=nodesets)
end

function create_nodes(rawnodes::RawNodes, ::Val{dim}) where{dim}
    num = getnumnodes(rawnodes)
    nodes=Array{Node{dim, Float64}}(undef, num)
    for (index, node_number) in enumerate(getnumbers(rawnodes))
        x=Vec{dim}(getcoordinate(rawnodes, index))
        nodes[node_number] = Node(x)
    end
    return nodes
end

function create_cells(rawelementsdict::Dict{String,RawElements}, user_elements::Dict, format)
    builtin_elements = get_element_type_dict(format)
    num_elements = sum(getnumelements.(values(rawelementsdict)))
    cells_generic = Array{Ferrite.AbstractCell}(undef, num_elements)
    for (key, rawelements) in rawelementsdict
        if haskey(user_elements, key)   # user_elements are prioritized over builtin
            addcells!(cells_generic, user_elements[key], rawelements, format)
        elseif haskey(builtin_elements, key)
            addcells!(cells_generic, builtin_elements[key], rawelements, format)
        else
            throw(UnsupportedElementType(key))
        end
    end
    cell_type = promote_type(unique(typeof.(cells_generic))...)
    cells = Array{cell_type}([cell for cell in cells_generic])  # Ensures that an as concrete type as possible returned.
    return cells
end

function addcells!(cells, elementtype, rawelements::RawElements, format)
    for (i, element_number) in enumerate(getnumbers(rawelements))
        node_numbers = gettopology(rawelements)[:,i]
        cells[element_number] = create_cell(elementtype, node_numbers, format)
    end
end

# Creating sets
function create_cellsets(rawelementsdict, rawelementsets)
    cellsets = Dict(key => Set(getnumbers(rawelements)) for (key, rawelements) in rawelementsdict)
    merge!(cellsets, Dict(key => Set(nums) for (key, nums) in rawelementsets))
    return cellsets
end

function create_nodesets(rawnodesets)
    return Dict(key => Set(nums) for (key, nums) in rawnodesets)
end

"""
    function generate_facesets!(grid::Ferrite.Grid)

Based on all nodesets in `grid`, generate facesets for those sets.
"""
function generate_facesets!(grid::Ferrite.Grid)
    cellset_or_nothing = Dict(key => (haskey(getcellsets(grid), key) ? getcellsets(grid)[key] : nothing)
                                for (key, set) in getnodesets(grid))
    merge!(getfacesets(grid), Dict(key => create_faceset(grid, set, cellset_or_nothing[key])
                                for (key, set) in getnodesets(grid)))
end


text:
type SharedArray{T,N} <: DenseArray{T,N}
    dims::NTuple{N,Int}
    pids::Vector{Int}
    refs::Array{RemoteRef}

    # The segname is currently used only in the test scripts to ensure that
    # the shmem segment has been unlinked.
    segname::String

    # Fields below are not to be serialized
    # Local shmem map.
    s::Array{T,N}

    # idx of current workers pid into the pids vector, 0 if this shared array is not mapped locally.
    pidx::Int

    # the local partition into the array when viewed as a single dimensional array.
    # this can be removed when @parallel or its equivalent supports looping on
    # a subset of workers.
    loc_subarr_1d

    SharedArray(d,p,r,sn) = new(d,p,r,sn)
end

function SharedArray(T::Type, dims::NTuple; init=false, pids=Int[])
    N = length(dims)

    isbits(T) || error("Type of Shared Array elements must be bits types")
    @windows_only error(" SharedArray is not supported on Windows yet.")

    if isempty(pids)
        # only use workers on the current host
        pids = procs(myid())
        onlocalhost = true
    else
        onlocalhost = assert_same_host(pids)
    end

    local shm_seg_name = ""
    local s
    local S = nothing
    local shmmem_create_pid
    try
        # On OSX, the shm_seg_name length must be < 32 characters
        shm_seg_name = string("/jl", getpid(), int64(time() * 10^9))
        if onlocalhost
            shmmem_create_pid = myid()
            s = shm_mmap_array(T, dims, shm_seg_name, JL_O_CREAT | JL_O_RDWR)
        else
            # The shared array is created on a remote machine....
            shmmem_create_pid = pids[1]
            remotecall_fetch(pids[1], () -> begin shm_mmap_array(T, dims, shm_seg_name, JL_O_CREAT | JL_O_RDWR); nothing end)
        end

        func_mapshmem = () -> shm_mmap_array(T, dims, shm_seg_name, JL_O_RDWR)

        refs = Array(RemoteRef, length(pids))
        for (i, p) in enumerate(pids)
            refs[i] = remotecall(p, func_mapshmem)
        end

        # Wait till all the workers have mapped the segment
        for i in 1:length(refs)
            wait(refs[i])
        end

        # All good, immediately unlink the segment.
        if onlocalhost
            shm_unlink(shm_seg_name)
        else
            remotecall_fetch(shmmem_create_pid, shm_unlink, shm_seg_name)
        end
        S = SharedArray{T,N}(dims, pids, refs, shm_seg_name)
        shm_seg_name = ""

        if onlocalhost
            init_loc_flds(S)

            # In the event that myid() is not part of pids, s will not be set
            # in the init function above, hence setting it here if available.
            S.s = s
        else
            S.pidx = 0
        end

        # if present init function is called on each of the parts
        if isa(init, Function)
            @sync begin
                for p in pids
                    @async remotecall_wait(p, init, S)
                end
            end
        end

    finally
        if shm_seg_name != ""
            remotecall_fetch(shmmem_create_pid, shm_unlink, shm_seg_name)
        end
    end
    S
end

SharedArray(T, I::Int...; kwargs...) = SharedArray(T, I; kwargs...)

length(S::SharedArray) = prod(S.dims)
size(S::SharedArray) = S.dims

procs(S::SharedArray) = S.pids
indexpids(S::SharedArray) = S.pidx

sdata(S::SharedArray) = S.s
sdata(A::AbstractArray) = A

localindexes(S::SharedArray) = S.pidx > 0 ? range_1dim(S, S.pidx) : error("SharedArray is not mapped to this process")

convert{T}(::Type{Ptr{T}}, S::SharedArray) = convert(Ptr{T}, sdata(S))

convert(::Type{SharedArray}, A::Array) = (S = SharedArray(eltype(A), size(A)); copy!(S, A))
convert{T}(::Type{SharedArray{T}}, A::Array) = (S = SharedArray(T, size(A)); copy!(S, A))
convert{TS,TA,N}(::Type{SharedArray{TS,N}}, A::Array{TA,N}) = (S = SharedArray(TS, size(A)); copy!(S, A))

function deepcopy_internal(S::SharedArray, stackdict::ObjectIdDict)
    haskey(stackdict, S) && return stackdict[S]
    # Note: copy can be used here because SharedArrays are restricted to isbits types
    R = copy(S)
    stackdict[S] = R
    return R
end

function range_1dim(S::SharedArray, pidx)
    l = length(S)
    nw = length(S.pids)
    partlen = div(l, nw)

    if l < nw
        if pidx <= l
            return pidx:pidx
        else
            return 1:0
        end
    elseif pidx == nw
        return (((pidx-1) * partlen) + 1):l
    else
        return (((pidx-1) * partlen) + 1):(pidx*partlen)
    end
end

sub_1dim(S::SharedArray, pidx) = sub(S.s, range_1dim(S, pidx))

function init_loc_flds(S)
    if myid() in S.pids
        S.pidx = findfirst(S.pids, myid())
        S.s = fetch(S.refs[S.pidx])
        S.loc_subarr_1d = sub_1dim(S, S.pidx)
    else
        S.pidx = 0
    end
end

# Don't serialize s (it is the complete array) and
# pidx, which is relevant to the current process only
function serialize(s, S::SharedArray)
    serialize_type(s, typeof(S))
    serialize(s, length(SharedArray.names))
    for n in SharedArray.names
        if n in [:s, :pidx, :loc_subarr_1d]
            writetag(s, UndefRefTag)
        else
            serialize(s, getfield(S, n))
        end
    end
end

function deserialize{T,N}(s, t::Type{SharedArray{T,N}})
    S = invoke(deserialize, (Any, DataType), s, t)
    init_loc_flds(S)
    if (S.pidx == 0)
        error("SharedArray cannot be used on a non-participating process")
    end
    S
end

convert(::Type{Array}, S::SharedArray) = S.s

# # pass through getindex and setindex! - they always work on the complete array unlike DArrays
getindex(S::SharedArray) = getindex(S.s)
getindex(S::SharedArray, I::Real) = getindex(S.s, I)
getindex(S::SharedArray, I::AbstractArray) = getindex(S.s, I)
@nsplat N 1:5 getindex(S::SharedArray, I::NTuple{N,Any}...) = getindex(S.s, I...)

setindex!(S::SharedArray, x) = (setindex!(S.s, x); S)
setindex!(S::SharedArray, x, I::Real) = (setindex!(S.s, x, I); S)
setindex!(S::SharedArray, x, I::AbstractArray) = (setindex!(S.s, x, I); S)
@nsplat N 1:5 setindex!(S::SharedArray, x, I::NTuple{N,Any}...) = (setindex!(S.s, x, I...); S)

# convenience constructors
function shmem_fill(v, dims; kwargs...)
    SharedArray(typeof(v), dims; init = S->fill!(S.loc_subarr_1d, v), kwargs...)
end
shmem_fill(v, I::Int...; kwargs...) = shmem_fill(v, I; kwargs...)

# rand variant with range
function shmem_rand(TR::Union(DataType, UnitRange), dims; kwargs...)
    if isa(TR, UnitRange)
        SharedArray(Int, dims; init = S -> map!((x)->rand(TR), S.loc_subarr_1d), kwargs...)
    else
        SharedArray(TR, dims; init = S -> map!((x)->rand(TR), S.loc_subarr_1d), kwargs...)
    end
end
shmem_rand(TR::Union(DataType, UnitRange), i::Int; kwargs...) = shmem_rand(TR, (i,); kwargs...)
shmem_rand(TR::Union(DataType, UnitRange), I::Int...; kwargs...) = shmem_rand(TR, I; kwargs...)

shmem_rand(dims; kwargs...) = shmem_rand(Float64, dims; kwargs...)
shmem_rand(I::Int...; kwargs...) = shmem_rand(I; kwargs...)

function shmem_randn(dims; kwargs...)
    SharedArray(Float64, dims; init = S-> map!((x)->randn(), S.loc_subarr_1d), kwargs...)
end
shmem_randn(I::Int...; kwargs...) = shmem_randn(I; kwargs...)

similar(S::SharedArray, T, dims::Dims) = SharedArray(T, dims; pids=procs(S))
similar(S::SharedArray, T) = similar(S, T, size(S))
similar(S::SharedArray, dims::Dims) = similar(S, eltype(S), dims)
similar(S::SharedArray) = similar(S, eltype(S), size(S))

map(f::Callable, S::SharedArray) = (S2 = similar(S); S2[:] = S[:]; map!(f, S2); S2)

reduce(f::Function, S::SharedArray) =
    mapreduce(fetch, f,
              { @spawnat p reduce(f, S.loc_subarr_1d) for p in procs(S) })

function map!(f::Callable, S::SharedArray)
    @sync for p in procs(S)
        @spawnat p begin
            for idx in localindexes(S)
                S.s[idx] = f(S.s[idx])
            end
        end
    end
    return S
end

copy!(S::SharedArray, A::Array) = copy!(S.s, A)

function copy!(S::SharedArray, R::SharedArray)
    length(S) == length(R) || throw(BoundsError())
    ps = intersect(procs(S), procs(R))
    isempty(ps) && error("source and destination arrays don't share any process")
    l = length(S)
    length(ps) > l && (ps = ps[1:l])
    nw = length(ps)
    partlen = div(l, nw)

    @sync for i = 1:nw
        p = ps[i]
        idx = i < nw ?  ((i-1)*partlen+1:i*partlen) : ((i-1)*partlen+1:l)
        @spawnat p begin
            S.s[idx] = R.s[idx]
        end
    end

    return S
end

function print_shmem_limits(slen)
    try
        @linux_only pfx = "kernel"
        @osx_only pfx = "kern.sysv"

        shmmax_MB = div(int(split(readall(readsfrom(`sysctl $(pfx).shmmax`)[1]))[end]), 1024*1024)
        page_size = int(split(readall(readsfrom(`getconf PAGE_SIZE`)[1]))[end])
        shmall_MB = div(int(split(readall(readsfrom(`sysctl $(pfx).shmall`)[1]))[end]) * page_size, 1024*1024)

        println("System max size of single shmem segment(MB) : ", shmmax_MB,
            "\nSystem max size of all shmem segments(MB) : ", shmall_MB,
            "\nRequested size(MB) : ", div(slen, 1024*1024),
            "\nPlease ensure requested size is within system limits.",
            "\nIf not, increase system limits and try again."
        )
    catch e
        nothing # Ignore any errors in this...
    end
end

# utilities
function shm_mmap_array(T, dims, shm_seg_name, mode)
    local s = nothing
    local A = nothing
    try
        fd_mem = shm_open(shm_seg_name, mode, S_IRUSR | S_IWUSR)
        systemerror("shm_open() failed for " * shm_seg_name, fd_mem <= 0)

        s = fdio(fd_mem, true)

        # On OSX, ftruncate must to used to set size of segment, just lseek does not work.
        # and only at creation time
        if (mode & JL_O_CREAT) == JL_O_CREAT
            rc = ccall(:ftruncate, Int, (Int, Int), fd_mem, prod(dims)*sizeof(T))
            systemerror("ftruncate() failed for shm segment " * shm_seg_name, rc != 0)
        end

        A = mmap_array(T, dims, s, zero(FileOffset), grow=false)
    catch e
        print_shmem_limits(prod(dims)*sizeof(T))
        rethrow(e)

    finally
        if s != nothing
            close(s)
        end
    end
    A
end

@unix_only begin
function shm_unlink(shm_seg_name)
    rc = ccall(:shm_unlink, Cint, (Ptr{Uint8},), shm_seg_name)
    systemerror("Error unlinking shmem segment " * shm_seg_name, rc != 0)
    rc
end
end

@unix_only shm_open(shm_seg_name, oflags, permissions) = ccall(:shm_open, Int, (Ptr{Uint8}, Int, Int), shm_seg_name, oflags, permissions)

function assert_same_host(procs)
    first_privip = getprivipaddr(procs[1])
    if !all(x -> getprivipaddr(x) == first_privip, procs)
        error("SharedArray requires all requested processes to be on the same machine.")
    end

    return myid() in procs
end

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
# Code derived from Robin Deits' RegionTrees and AdaptiveDistanceFields
# The RegionTrees.jl package is licensed under the MIT "Expat" License
# The AdaptiveDistanceFields.jl package is licensed under the MIT "Expat" License

# this is designed to avoid allocation of a full octtree and generate a mesh
# while sampling

function vertices(h::GB.HyperRectangle, ::Type{SV}) where SV
    z = zero(eltype(SV))
    @inbounds begin
        o = SV(h.origin[1],h.origin[2],h.origin[3])
        w = SV(h.widths[1],h.widths[2],h.widths[3])
        (o,
        o.+SV(w[1],z,z),
        o.+SV(w[1],w[2],z),
        o.+SV(z,w[2],z),
        o.+SV(z,z,w[3]),
        o.+SV(w[1],z,w[3]),
        o.+w,
        o.+SV(z,w[2],w[3]))
    end
end

function vertices_mt(h::GB.HyperRectangle, ::Type{SV}) where SV
    z = zero(eltype(SV))
    @inbounds begin
        o = SV(h.origin[1],h.origin[2],h.origin[3])
        w = SV(h.widths[1],h.widths[2],h.widths[3])
        (o,
        o.+SV(   z,w[2],   z),
        o.+SV(w[1],w[2],   z),
        o.+SV(w[1],   z,   z),
        o.+SV(   z,   z,w[3]),
        o.+SV(   z,w[2],w[3]),
        o.+w,
        o.+SV(w[1],   z,w[3]))
    end
end

function interpolate_mt(c)
    @inbounds begin
        (sum(c)*0.125, # center
        (c[1]+c[2]+c[3]+c[4])*0.25,
        (c[1]+c[4]+c[5]+c[8])*0.25,
        (c[1]+c[2]+c[5]+c[6])*0.25,
        (c[5]+c[6]+c[7]+c[8])*0.25,
        (c[2]+c[3]+c[6]+c[7])*0.25,
        (c[3]+c[4]+c[7]+c[8])*0.25)
    end
end

function face_center_vertices(h::GB.HyperRectangle{N,T}) where{N,T}
    SV = SVector{3,T}
    z = zero(T)
    @inbounds begin
        o = SV(h.origin[1],h.origin[2],h.origin[3])
        w = SV(h.widths[1],h.widths[2],h.widths[3])
        hw = w.*0.5
        c = center(h)
        (SV(c[1],c[2],c[3]),
        o.+SV(hw[1],hw[2],    z),
        o.+SV(hw[1],    z,hw[3]),
        o.+SV(z    ,hw[2],hw[3]),
        o.+SV(hw[1],hw[2], w[3]),
        o.+SV(hw[1], w[2],hw[3]),
        o.+SV( w[1],hw[2],hw[3]))
    end
end

center(rect::GB.HyperRectangle) = rect.origin + 0.5 * rect.widths

function octsplit(h::GB.HyperRectangle)
    ET = eltype(h.origin)
    VT = Vec{3,ET}
    z = zero(ET)
    hw = h.widths
    nw = h.widths.*ET(0.5)
    no = h.origin
    @inbounds begin
        (GB.HyperRectangle(no, nw),
        GB.HyperRectangle(no.+nw, nw),
        GB.HyperRectangle(no.+VT(nw[1],z,z), nw),
        GB.HyperRectangle(no.+VT(z,nw[2],z), nw),
        GB.HyperRectangle(no.+VT(z,z,nw[3]), nw),
        GB.HyperRectangle(no.+VT(nw[1],nw[2],z), nw),
        GB.HyperRectangle(no.+VT(nw[1],z,nw[3]), nw),
        GB.HyperRectangle(no.+VT(z,nw[2],nw[3]), nw))
    end
end

function get_iso_vals(f, val_store, points::NTuple{N,T}) where {N, T}
    @inbounds begin
        return ntuple(N) do i
            get_iso_vals(f,val_store,points[i])
        end
    end
end

function get_iso_vals(f, val_store, pt)
    # tuples are faster to hash than SVec
    key_tup = pt.data # tuple in SVectors/GeometryTypes, probably shoudl make more generic
    if haskey(val_store,key_tup)
        val_store[key_tup]
    else
        val_store[key_tup] = f(pt)::valtype(val_store)
    end
end

@inline function _get_interpindex(iso_vals, iso)
    cubeindex = iso_vals[1] < iso ? 0x01 : 0x00
    iso_vals[2] < iso && (cubeindex |= 0x02)
    iso_vals[3] < iso && (cubeindex |= 0x04)
    iso_vals[4] < iso && (cubeindex |= 0x08)
    iso_vals[5] < iso && (cubeindex |= 0x10)
    iso_vals[6] < iso && (cubeindex |= 0x20)
    iso_vals[7] < iso && (cubeindex |= 0x40)
    cubeindex
end

"""
    _get_cubeindex_pos(iso_vals, iso)

given `iso_vals` and iso, return an 8 bit value corresponding
to each corner of a cube. In each bit position,
0 indicates in the isosurface and 1 indicates outside the surface,
where the sign convention indicates positive inside the surface
"""
@inline function _get_interpindex_pos(iso_vals, iso)
    cubeindex = iso_vals[1] > iso ? 0x01 : 0x00
    iso_vals[2] > iso && (cubeindex |= 0x02)
    iso_vals[3] > iso && (cubeindex |= 0x04)
    iso_vals[4] > iso && (cubeindex |= 0x08)
    iso_vals[5] > iso && (cubeindex |= 0x10)
    iso_vals[6] > iso && (cubeindex |= 0x20)
    iso_vals[7] > iso && (cubeindex |= 0x40)
    cubeindex
end

function isosurface(f::Function, method::AdaptiveMarchingCubes, ::Type{VertType}=SVector{3,Float64}, ::Type{FaceType}=SVector{3, Int};
                    origin=VertType(-1,-1,-1), widths=VertType(2,2,2)) where {VertType, FaceType}

    ET = eltype(VertType)

    # arrays for vertices and faces
    vts = VertType[]
    fcs = FaceType[]

    # refinement queue
    refinement_queue = GB.HyperRectangle{3,ET}[GB.HyperRectangle{3,ET}(origin,widths)]

    val_store = Dict{NTuple{3,ET},ET}();

    @inbounds while true

        if isempty(refinement_queue)
            break
        end

        cell = pop!(refinement_queue)
        points = vertices(cell, VertType)

        iso_vals = get_iso_vals(f,val_store,points)

        #Determine the index into the edge table which
        #tells us which vertices are inside of the surface
        cubeindex = method.insidepositive ? _get_cubeindex_pos(iso_vals, method.iso) : _get_cubeindex(iso_vals, method.iso)
        #interpindex = method.insidepositive ? _get_interpindex_pos(iso_vals, method.iso) : _get_interpindex(iso_vals, method.iso)

        value_interp = sum(iso_vals)*0.125
        value_true = get_iso_vals(f, val_store, center(cell))

        if (cubeindex == 0xff && value_true < 0) || (iszero(cubeindex) && value_true > 0)
            continue
        elseif minimum(cell.widths) > method.atol && !isapprox(value_interp, value_true, rtol=method.rtol, atol=method.atol)
            append!(refinement_queue, octsplit(cell))
        else
            # Find the vertices where the surface intersects the cube
            # The underlying space is non-linear so there will be error otherwise
            vertlist = find_vertices_interp(points, iso_vals, cubeindex, method.iso, method.eps)

            # Create the triangle
            method.reduceverts && _mc_unique_triangles!(vts, fcs, vertlist, cubeindex, FaceType)
            !method.reduceverts && _mc_create_triangles!(vts, fcs, vertlist, cubeindex, FaceType)
        end
    end
    vts,fcs
end

@inline function vertPos(f, e, width, origin, ::Type{VertType}) where {VertType}
    T = eltype(eltype(VertType))

    ixs     = voxEdgeCrnrs[e]
    c1 = voxCrnrPos(VertType)[ixs[1]].*width .+ origin
    c2 = voxCrnrPos(VertType)[ixs[2]].*width .+ origin
    a = brent(f, c1, c2.-c1) # find root using brents method

    c1 + a.*(c2.-c1)
end

@inline function getVertId(f, e, width, vals, iso::Real, origin, vtsAry::Vector, vertex_store, eps::Real)

    VertType = eltype(vtsAry)

    # calculate vert position
    v = vertPos(f, e, width, origin, VertType)
    vt_key = v.data
    if haskey(vertex_store, vt_key)
        return vertex_store[vt_key]
    else
        push!(vtsAry, v)
        return vertex_store[vt_key] = length(vtsAry)
    end
end

slivers(x,y,z) = x == y || y == z || x == z

"""
    procVox(vals, iso::Real, x, y, z, nx, ny,
                    vts::Dict, vtsAry::Vector, fcs::Vector,
                    eps::Real)

Processes a voxel, adding any new vertices and faces to the given
containers as necessary.
"""
function procVox(f, vals, iso::Real, width, origin, vtsAry::Vector, vertex_store, fcs::Vector,
                 eps::Real, cubeindex)
    VertType = eltype(vtsAry)
    FaceType = eltype(fcs)
    # check each sub-tetrahedron in the voxel
    @inbounds for i = 1:6
        tIx = tetIx(i, cubeindex)
        (tIx == 0x00 || tIx == 0x0f) && continue

        e = tetTri[tIx]
        v1 = getVertId(f, voxEdgeId(i, e[1]), width, vals, iso, origin, vtsAry, vertex_store, eps)
        v2 = getVertId(f, voxEdgeId(i, e[2]), width, vals, iso, origin, vtsAry, vertex_store, eps)
        v3 = getVertId(f, voxEdgeId(i, e[3]), width, vals, iso, origin, vtsAry, vertex_store, eps)

        # add the face to the list
        !slivers(v1,v2,v3) && push!(fcs, FaceType(v1,v2,v3))

        # bail if there are no more faces
        iszero(e[4]) && continue
        v4 = getVertId(f, voxEdgeId(i, e[4]), width, vals, iso, origin, vtsAry, vertex_store, eps)
        v5 = getVertId(f, voxEdgeId(i, e[5]), width, vals, iso, origin, vtsAry, vertex_store, eps)
        v6 = getVertId(f, voxEdgeId(i, e[6]), width, vals, iso, origin, vtsAry, vertex_store, eps)
        !slivers(v4,v5,v6) && push!(fcs, FaceType(v4,v5,v6))

    end
end

function isosurface(f::Function, method::AdaptiveMarchingTetrahedra, ::Type{VertType}=SVector{3,Float64}, ::Type{FaceType}=SVector{3, Int};
                    origin=VertType(-1,-1,-1), widths=VertType(2,2,2)) where {VertType, FaceType}

    ET = eltype(VertType)

    # arrays for vertices and faces
    vts = VertType[]
    fcs = FaceType[]

    # refinement queue
    # initialize with depth two
    refinement_queue = GB.HyperRectangle{3,ET}[]
    for h in octsplit(GB.HyperRectangle{3,ET}(origin,widths))
        append!(refinement_queue,octsplit(h))
    end

    val_store = Dict{NTuple{3,ET},ET}();
    vertex_store = Dict{NTuple{3,ET},ET}();
    sizehint!(val_store, trunc(Int, maximum(widths)*12/max(method.atol,method.rtol)))
    sizehint!(vertex_store, trunc(Int, maximum(widths)*12/max(method.atol,method.rtol)))

    @inbounds while true

        if isempty(refinement_queue)
            break
        end

        cell = pop!(refinement_queue)
        points = vertices_mt(cell, VertType)
        interp_points = face_center_vertices(cell)

        iso_vals = get_iso_vals(f,val_store,points)
        true_vals = get_iso_vals(f,val_store,interp_points)

        #Determine the index into the edge table which
        #tells us which vertices are inside of the surface
        cubeindex = method.insidepositive ? _get_cubeindex_pos(iso_vals, method.iso) : _get_cubeindex(iso_vals, method.iso)
        interpindex = method.insidepositive ? _get_interpindex_pos(true_vals, method.iso) : _get_interpindex(true_vals, method.iso)

        if (cubeindex == 0xff && interpindex == 0x7f) || (iszero(cubeindex) && iszero(interpindex))
            continue
        else
            value_interp = interpolate_mt(iso_vals)
            accurate = true
            for i = 1:7
                if !isapprox(true_vals[i], value_interp[i], rtol=method.rtol, atol=method.atol)
                    append!(refinement_queue, octsplit(cell))
                    accurate = false
                    break
                end
            end
            if accurate
                procVox(f, iso_vals, method.iso, cell.widths, cell.origin, vts, vertex_store, fcs, method.eps, cubeindex)
            end
        end
    end
    vts,fcs
end


text:
# More convenient constructors for types and conversions to types

RealQuantity(value::Real) = RealQuantity(value=value)
IntegerQuantity(value::Number) = IntegerQuantity(value=value)

# Allow filling of fields which are types with only a single field
# e.g., EventParameters(public_id="example string")
# rather than EventParameters(public_id=ResourceReference("example string"))
for T in Base.uniontypes(ValueTypes)
    @eval Base.convert(::Type{$T}, s::AbstractString) = $T(s)
end

# Allow contruction of types which use the following more easily.
# e.g.,
#     QuakeML.Origin(time=DateTime("2012-01-01T00:00:00"),
#         longitude=1, latitude=2, public_id="smi:a.com/b")
# rather than
#     QuakeML.Origin(time=TimeQuantity(value=DateTime("2012-01-01T00:00:00")),
#         longitude=RealQuantity(value=1), etc...)
Base.convert(::Type{RealQuantity}, value::Real) = RealQuantity(value=value)
Base.convert(::Type{IntegerQuantity}, value::Integer) = IntegerQuantity(value=value)
Base.convert(::Type{TimeQuantity}, value::DateTime) = TimeQuantity(value=value)

max_stars_repo_name: anowacki/QuakeML


text:
module mmmmmmmmmmmm
using FinEtools
using Compat.Test
function test()
    println("""
    pnpConcreteColumnsymb_conv.
    Two-dimensional heat transfer with convection.
    """
    )

    a=2.5; dy=a/2*sin(15/180*pi); dx=a/2*cos(15/180*pi); Q=4.5; k=1.8; Dz=1.0;
    h= 5.;

    m = MatHeatDiff([k 0; 0 k])

    modeldata = nothing

    fens = FENodeSet([0 0; dx -dy; dx dy; 2*dx -2*dy; 2*dx 2*dy])
    fes = FESetT3([1 2 3; 2 4 5; 2 5 3]);

    geom = NodalField(fens.xyz)
    Temp = NodalField(zeros(size(fens.xyz,1),1))

    applyebc!(Temp)
    numberdofs!(Temp)
    Temp.dofnums[1] = 3
    Temp.dofnums[2] = 1
    Temp.dofnums[3] = 2
    Temp.dofnums[4] = 5
    Temp.dofnums[5] = 4

    bfes = FESetL2([4 5]);

    cfemm = FEMMHeatDiffSurf(IntegData(bfes, GaussRule(1, 2), Dz), h)
    femm = FEMMHeatDiff(IntegData(fes, TriRule(1), Dz), m)
    fi = ForceIntensity(FFlt[Q]);
    F1 = distribloads(femm, geom, Temp, fi, 3);
    K = conductivity(femm, geom, Temp)
    F2 = nzebcloadsconductivity(femm, geom, Temp);
    H = surfacetransfer(cfemm, geom, Temp);
    show(H)
    # F3 = surfacetransferloads(femm, geom, temp, amb);
    Factor = cholesky(K+H)
    U = Factor\(F1+F2)
    scattersysvec!(Temp, U[:])
    display(Temp)

end
end
using .mmmmmmmmmmmm
mmmmmmmmmmmm.test()


text:
# Example: how to read mds format output on poseidon.whoi.edu.
using MITgcmTools, MeshArrays

path = "/poseidon/ecco/ECCOv4r4/MITgcm/exps/iter129_bulkformula/run/"
θfile = "T.0000000001.data" # it's actually potential temperature snapshot
@time θmds     = read_mdsio(path,θfile); # output is a 3D array

# To use a gcmarray variable type, need to have grid information in γ
@time θstart   = γ.read(path*θfile,MeshArray(γ,Float32));

# Do same thing with final snapshot.
θfile = "T.0000227904.data"
θend=γ.read(path*θfile,MeshArray(γ,Float32))

# operations proceed normally
Δθ = θend - θstart
maximum(ΔΘ) # etc.

# compare to reading monthly average fields
path = "/poseidon/ecco/ECCOv4r4/MITgcm/exps/iter129_fluxforced/run/diags/"
θfile = "state_3d_set1.0000000732.data"
@time xmds     = read_mdsio(path,θfile);

# To use a gcmarray variable type, need to have grid information in γ
@time xstart   = γ.read(path*θfile,MeshArray(γ,Float32));


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Compute the sun position.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] The Astronomical Almanac for the year 2000 (p. C24).
#
#   [2] http://aa.usno.navy.mil/faq/docs/SunApprox.php
#
#   [3] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorne, CA.
#
#   [4] The Astronomical Almanac for the year 2006.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export sun_position_i

"""
    sun_position_i(JD::Number)

Compute the Sun position represented in the Mean Equinox of Date (MOD) at the
Julian Day `JD`. The algorithm was adapted from [3, p. 277-279].

"""
function sun_position_i(JD::Number)
    # Constants.
    deg2rad = π/180

    # Number of Julian centuries from J2000 epoch.
    T_UT1 = (JD-JD_J2000)/36525.0

    # Here, we will assume that T_TBD ≈ T_UT1. This can be done because we are
    # using a low-precision algorithm [3].
    T_TBD = T_UT1

    # Mean anomaly of the Sun [deg].
    Ms = 357.529_109_2 + 35_999.050_34T_TBD

    # Convert Ms to [rad] and limit to the interval [0,2π].
    Ms = mod2pi(Ms*deg2rad)

    # Compute auxiliary variables.
    sinMs,  cosMs  = sincos(Ms)
    sin2Ms, cos2Ms = sincos(2Ms)

    # Mean longitude of the Sun [deg].
    λ_m = 280.460 + 36_000.771T_UT1

    # Ecliptic latitude of the Sun [deg].
    λ_e = λ_m + 1.914_666_471sinMs + 0.019_994_643sin2Ms

    # Obliquity of the ecliptic [deg].
    ϵ = 23.439_291 - 0.013_004_2T_TBD

    # Convert λ_e and ϵ to [rad] and limit to the interval [0,2π].
    λ_e = mod2pi(λ_e*deg2rad)
    ϵ   = mod2pi(  ϵ*deg2rad)

    # Auxiliary variables.
    sinϵ   , cosϵ   = sincos(ϵ)
    sinλ_e , cosλ_e = sincos(λ_e)

    # Distance of the Sun from Earth [m].
    r = ( 1.000_140_612 - 0.016_708_617cosMs - 0.000_139_589cos2Ms )*au2m

    # Compute the Sun vector represented in the Mean Equinox of Date (MOD).
    S_MOD = SVector{3,Float64}(r*cosλ_e, r*cosϵ*sinλ_e, r*sinϵ*sinλ_e)

    # TODO: This vector must be transformed to the inertial reference frame that
    # will be used in this toolbox (J2000, True of Date, etc.).
end


text:
function add_entry!(names, values, entry; default)
    i = findfirst(==(entry), names)
    if isnothing(i)
        push!(names, entry)
        push!(values, default)
        i = lastindex(names)
    end
    return values[i]
end

struct LegendSection
    title::String
    names::Vector{String}
    plots::Vector{Vector{AbstractPlot}}
end
LegendSection(title::String="") = LegendSection(title, String[], Vector{AbstractPlot}[])

# Add an empty trace list with name `entry` to the legend section
function add_entry!(legendsection::LegendSection, entry::String)
    names, plots = legendsection.names, legendsection.plots
    return add_entry!(names, plots, entry; default=AbstractPlot[])
end

struct Legend
    names::Vector{String}
    sections::Vector{LegendSection}
end
Legend() = Legend(String[], LegendSection[])

# Add an empty section with name `entry` and title `title` to the legend
function add_entry!(legend::Legend, entry::String; title::String="")
    names, sections = legend.names, legend.sections
    return add_entry!(names, sections, entry; default=LegendSection(title))
end

function create_legend(scene, legend::Legend)
    legend = remove_duplicates(legend)
    sections = legend.sections
    MakieLayout.LLegend(
        scene,
        getproperty.(sections, :plots),
        getproperty.(sections, :names),
        # LLegend needs `nothing` to remove the space for a missing title
        [t == " " ? nothing : t for t in getproperty.(sections, :title)]
    )
end

function remove_duplicates(legend)
    sections = legend.sections
    titles = getproperty.(sections, :title)
    # check if there are duplicate titles
    unique_inds = unique_indices(titles; keep = " ")
    has_duplicates = length(unique_inds) < length(titles)
    # if so: remove duplicates, generate new names
    if has_duplicates
        sections_new = sections[unique_inds]
        names_new = legend.names[unique_inds]
        return Legend(names_new, sections_new)
    else
        return legend
    end
end

function unique_indices(x; keep)
    first_inds = indexin(x, x)
    inds_keep = findall(==(keep), x)
    sort!(union!(inds_keep, first_inds))
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

mutable struct Set{T} <: AbstractSet{T}
    dict::Dict{T,Void}

    Set{T}() where {T} = new(Dict{T,Void}())
    Set{T}(itr) where {T} = union!(new(Dict{T,Void}()), itr)
end
Set() = Set{Any}()

"""
    Set([itr])

Construct a [`Set`](@ref) of the values generated by the given iterable object, or an
empty set. Should be used instead of [`IntSet`](@ref) for sparse integer sets, or
for sets of arbitrary objects.
"""
Set(itr) = Set{eltype(itr)}(itr)
function Set(g::Generator)
    T = _default_eltype(typeof(g))
    (_isleaftype(T) || T === Union{}) || return grow_to!(Set{T}(), g)
    return Set{T}(g)
end

eltype(::Type{Set{T}}) where {T} = T
similar(s::Set{T}) where {T} = Set{T}()
similar(s::Set, T::Type) = Set{T}()

function show(io::IO, s::Set)
    print(io, "Set")
    if isempty(s)
        print(io, "{", eltype(s), "}()")
        return
    end
    print(io, "(")
    show_vector(io, s, "[", "]")
    print(io, ")")
end

isempty(s::Set) = isempty(s.dict)
length(s::Set)  = length(s.dict)
in(x, s::Set) = haskey(s.dict, x)
push!(s::Set, x) = (s.dict[x] = nothing; s)
pop!(s::Set, x) = (pop!(s.dict, x); x)
pop!(s::Set, x, deflt) = x in s ? pop!(s, x) : deflt

function pop!(s::Set)
    isempty(s) && throw(ArgumentError("set must be non-empty"))
    idx = start(s.dict)
    val = s.dict.keys[idx]
    _delete!(s.dict, idx)
    val
end

delete!(s::Set, x) = (delete!(s.dict, x); s)

copy(s::Set) = union!(similar(s), s)

sizehint!(s::Set, newsz) = (sizehint!(s.dict, newsz); s)
empty!(s::Set) = (empty!(s.dict); s)
rehash!(s::Set) = (rehash!(s.dict); s)

start(s::Set)       = start(s.dict)
done(s::Set, state) = done(s.dict, state)
# NOTE: manually optimized to take advantage of Dict representation
next(s::Set, i)     = (s.dict.keys[i], skip_deleted(s.dict, i+1))

"""
    union(s1,s2...)
    ∪(s1,s2...)

Construct the union of two or more sets. Maintains order with arrays.

# Examples
```jldoctest
julia> union([1, 2], [3, 4])
4-element Array{Int64,1}:
 1
 2
 3
 4

julia> union([1, 2], [2, 4])
3-element Array{Int64,1}:
 1
 2
 4

julia> union([4, 2], [1, 2])
3-element Array{Int64,1}:
 4
 2
 1
```
"""
function union end

union(s::Set) = copy(s)
function union(s::Set, sets...)
    u = Set{join_eltype(s, sets...)}()
    union!(u, s)
    for t in sets
        union!(u, t)
    end
    return u
end
const ∪ = union

"""
    union!(s, iterable)

Union each element of `iterable` into set `s` in-place.

# Examples
```jldoctest
julia> a = Set([1, 3, 4, 5]);

julia> union!(a, 1:2:8);

julia> a
Set([7, 4, 3, 5, 1])
```
"""
function union!(s::Set{T}, xs) where T
    haslength(xs) && sizehint!(s, length(xs))
    for x=xs
        push!(s, x)
        length(s) == max_values(T) && break
    end
    s
end

join_eltype() = Bottom
join_eltype(v1, vs...) = typejoin(eltype(v1), join_eltype(vs...))

"""
    intersect(s1,s2...)
    ∩(s1,s2)

Construct the intersection of two or more sets.
Maintains order and multiplicity of the first argument for arrays and ranges.

# Examples
```jldoctest
julia> intersect([1, 2, 3], [3, 4, 5])
1-element Array{Int64,1}:
 3

julia> intersect([1, 4, 4, 5, 6], [4, 6, 6, 7, 8])
3-element Array{Int64,1}:
 4
 4
 6
```
"""
function intersect end

intersect(s::Set) = copy(s)
function intersect(s::Set, sets...)
    i = similar(s)
    for x in s
        inall = true
        for t in sets
            if !in(x, t)
                inall = false
                break
            end
        end
        inall && push!(i, x)
    end
    return i
end
const ∩ = intersect

function setdiff(a::Set, b)
    d = similar(a)
    for x in a
        if !(x in b)
            push!(d, x)
        end
    end
    d
end

"""
    setdiff!(s, iterable)

Remove each element of `iterable` from set `s` in-place.

# Examples
```jldoctest
julia> a = Set([1, 3, 4, 5]);

julia> setdiff!(a, 1:2:6);

julia> a
Set([4])
```
"""
setdiff!(s::Set, xs) = (for x=xs; delete!(s, x); end; s)

==(l::Set, r::Set) = (length(l) == length(r)) && (l <= r)
<( l::Set, r::Set) = (length(l) < length(r)) && (l <= r)
<=(l::Set, r::Set) = issubset(l, r)

"""
    issubset(a, b)
    ⊆(a,b) -> Bool
    ⊈(a,b) -> Bool
    ⊊(a,b) -> Bool

Determine whether every element of `a` is also in `b`, using [`in`](@ref).

# Examples
```jldoctest
julia> issubset([1, 2], [1, 2, 3])
true

julia> issubset([1, 2, 3], [1, 2])
false
```
"""
function issubset(l, r)
    for elt in l
        if !in(elt, r)
            return false
        end
    end
    return true
end
const ⊆ = issubset
⊊(l::Set, r::Set) = <(l, r)
⊈(l::Set, r::Set) = !⊆(l, r)
⊇(l, r) = issubset(r, l)
⊉(l::Set, r::Set) = !⊇(l, r)
⊋(l::Set, r::Set) = <(r, l)

"""
    unique(itr)

Return an array containing only the unique elements of collection `itr`,
as determined by [`isequal`](@ref), in the order that the first of each
set of equivalent elements originally appears. The element type of the
input is preserved.

# Examples
```jldoctest
julia> unique([1, 2, 6, 2])
3-element Array{Int64,1}:
 1
 2
 6

julia> unique(Real[1, 1.0, 2])
2-element Array{Real,1}:
 1
 2
```
"""
function unique(itr)
    T = _default_eltype(typeof(itr))
    out = Vector{T}()
    seen = Set{T}()
    i = start(itr)
    if done(itr, i)
        return out
    end
    x, i = next(itr, i)
    if !_isleaftype(T) && iteratoreltype(itr) == EltypeUnknown()
        S = typeof(x)
        return _unique_from(itr, S[x], Set{S}((x,)), i)
    end
    push!(seen, x)
    push!(out, x)
    return unique_from(itr, out, seen, i)
end

_unique_from(itr, out, seen, i) = unique_from(itr, out, seen, i)
@inline function unique_from(itr, out::Vector{T}, seen, i) where T
    while !done(itr, i)
        x, i = next(itr, i)
        S = typeof(x)
        if !(S === T || S <: T)
            R = typejoin(S, T)
            seenR = convert(Set{R}, seen)
            outR = convert(Vector{R}, out)
            if !in(x, seenR)
                push!(seenR, x)
                push!(outR, x)
            end
            return _unique_from(itr, outR, seenR, i)
        end
        if !in(x, seen)
            push!(seen, x)
            push!(out, x)
        end
    end
    return out
end

"""
    unique(f, itr)

Returns an array containing one value from `itr` for each unique value produced by `f`
applied to elements of `itr`.

# Examples
```jldoctest
julia> unique(x -> x^2, [1, -1, 3, -3, 4])
3-element Array{Int64,1}:
 1
 3
 4
```
"""
function unique(f::Callable, C)
    out = Vector{eltype(C)}()
    seen = Set()
    for x in C
        y = f(x)
        if !in(y, seen)
            push!(seen, y)
            push!(out, x)
        end
    end
    out
end

# If A is not grouped, then we will need to keep track of all of the elements that we have
# seen so far.
function _unique!(A::AbstractVector)
    seen = Set{eltype(A)}()
    idxs = eachindex(A)
    i = state = start(idxs)
    for x in A
        if x ∉ seen
            push!(seen, x)
            i, state = next(idxs, state)
            A[i] = x
        end
    end
    resize!(A, i - first(idxs) + 1)
end

# If A is grouped, so that each unique element is in a contiguous group, then we only
# need to keep track of one element at a time. We replace the elements of A with the
# unique elements that we see in the order that we see them. Once we have iterated
# through A, we resize A based on the number of unique elements that we see.
function _groupedunique!(A::AbstractVector)
    isempty(A) && return A
    idxs = eachindex(A)
    y = first(A)
    state = start(idxs)
    i, state = next(idxs, state)
    for x in A
        if !isequal(x, y)
            i, state = next(idxs, state)
            y = A[i] = x
        end
    end
    resize!(A, i - first(idxs) + 1)
end

"""
    unique!(A::AbstractVector)

Remove duplicate items as determined by [`isequal`](@ref), then return the modified `A`.
`unique!` will return the elements of `A` in the order that they occur. If you do not care
about the order of the returned data, then calling `(sort!(A); unique!(A))` will be much
more efficient as long as the elements of `A` can be sorted.

# Examples
```jldoctest
julia> unique!([1, 1, 1])
1-element Array{Int64,1}:
 1

julia> A = [7, 3, 2, 3, 7, 5];

julia> unique!(A)
4-element Array{Int64,1}:
 7
 3
 2
 5

julia> B = [7, 6, 42, 6, 7, 42];

julia> sort!(B);  # unique! is able to process sorted data much more efficiently.

julia> unique!(B)
3-element Array{Int64,1}:
  6
  7
 42
```
"""
function unique!(A::Union{AbstractVector{<:Real}, AbstractVector{<:AbstractString},
                          AbstractVector{<:Symbol}})
    if isempty(A)
        return A
    elseif issorted(A) || issorted(A, rev=true)
        return _groupedunique!(A)
    else
        return _unique!(A)
    end
end
# issorted fails for some element types, so the method above has to be restricted to
# elements with isless/< defined.
function unique!(A)
    if isempty(A)
        return A
    else
        return _unique!(A)
    end
end

"""
    allunique(itr) -> Bool

Return `true` if all values from `itr` are distinct when compared with [`isequal`](@ref).

# Examples
```jldoctest
julia> a = [1; 2; 3]
3-element Array{Int64,1}:
 1
 2
 3

julia> allunique([a, a])
false
```
"""
function allunique(C)
    seen = Set{eltype(C)}()
    for x in C
        if in(x, seen)
            return false
        else
            push!(seen, x)
        end
    end
    true
end

allunique(::Set) = true

allunique(r::AbstractRange{T}) where {T} = (step(r) != zero(T)) || (length(r) <= 1)

function filter(f, s::Set)
    u = similar(s)
    for x in s
        if f(x)
            push!(u, x)
        end
    end
    return u
end
function filter!(f, s::Set)
    for x in s
        if !f(x)
            delete!(s, x)
        end
    end
    return s
end

const hashs_seed = UInt === UInt64 ? 0x852ada37cfe8e0ce : 0xcfe8e0ce
function hash(s::Set, h::UInt)
    h = hash(hashs_seed, h)
    for x in s
        h ⊻= hash(x)
    end
    return h
end

convert(::Type{Set{T}}, s::Set{T}) where {T} = s
convert(::Type{Set{T}}, x::Set) where {T} = Set{T}(x)

max_stars_repo_name: visr/julia

max_issues_repo_name: visr/julia

max_forks_repo_name: visr/julia


text:
# Much is taken from ExpressionUtils.jl but generalized to work with ExLike

using Core: MethodInstance
using Base: MethodList

const poppable_macro = (Symbol("@inline"), Symbol("@noinline"), Symbol("@propagate_inbounds"), Symbol("@eval"))

"""
    exf = funcdef_expr(ex)

Recurse, if necessary, into `ex` until the first function definition expression is found.

# Example

```jldoctest; setup=(using Revise), filter=r"#=.*=#"
julia> Revise.funcdef_expr(quote
       \"\"\"
       A docstring
       \"\"\"
       @inline foo(x) = 5
       end)
:(foo(x) = begin
          #= REPL[31]:5 =#
          5
      end)
```
"""
function funcdef_expr(ex)
    if ex.head == :macrocall
        if ex.args[1] isa GlobalRef && ex.args[1].name == Symbol("@doc")
            return funcdef_expr(ex.args[end])
        elseif ex.args[1] ∈ poppable_macro
            return funcdef_expr(ex.args[end])
        else
            io = IOBuffer()
            dump(io, ex)
            throw(ArgumentError(string("unrecognized macro expression:\n", String(take!(io)))))
        end
    end
    if is_trivial_block_wrapper(ex)
        return funcdef_expr(first(LineSkippingIterator(ex.args)))
    end
    if ex.head == :function || ex.head == :(=)
        return ex
    end
    dump(ex)
    throw(ArgumentError(string("expected function definition expression, got ", ex)))
end

function funcdef_body(ex)
    fex = funcdef_expr(ex)
    if fex.head == :function || fex.head == :(=)
        return fex.args[end]
    end
    throw(ArgumentError(string("expected function definition expression, got ", ex)))
end

"""
    sigex = get_signature(expr)

Extract the signature from an expression `expr` that defines a function.

If `expr` does not define a function, returns `nothing`.

# Examples

```jldoctest; setup = :(using Revise)
julia> Revise.get_signature(quote
       function count_different(x::AbstractVector{T}, y::AbstractVector{S}) where {S,T}
           sum(x .!= y)
       end
       end)
:(count_different(x::AbstractVector{T}, y::AbstractVector{S}) where {S, T})
```
"""
function get_signature(ex::E) where E <: ExLike
    if ex.head == :macrocall
        error("macros must be handled externally, got $ex")
    end
    if is_trivial_block_wrapper(ex)
        return get_signature(ex.args[end])
    end
    if ex.head == :function
        return ex.args[1]
    elseif ex.head == :(=) && isa(ex.args[1], E)
        ex = ex.args[1]::E
        if ex.head == :where || ex.head == :call
            return ex
        end
    end
    nothing
end

@noinline sigwarn(mod, sigtexs, def) = @warn "error processing module $mod signature expressions $sigtexs from $def"
@noinline sigwarn(mod, sigtexs, ::Nothing) = @warn "error processing module $mod signature expressions $sigtexs"

function sigex2sigts(mod::Module, sig::ExLike, def=nothing)
    # Generate the signature-types
    local sigtexs
    try
        sigtexs = sig_type_exprs(mod, sig)
    catch err
        sigwarn(mod, sig, def)
        rethrow(err)
    end
    local sigts
    try
        sigts = Any[Core.eval(mod, s) for s in sigtexs]
    catch err
        sigwarn(mod, sigtexs, def)
        rethrow(err)
    end
    return sigts
end

"""
    methsig = sigt2methsig(sig)

For a signature `sig`, try to return the signature `methsig` of a currently-defined method.
"""
function sigt2methsig(sig)
    ret = Base._methods_by_ftype(sig, -1, typemax(UInt))
    isempty(ret) && return sig
    methsig = ret[end][3].sig  # the last method returned is the least-specific that matches, and thus most likely to be type-equal
    @assert sig <: methsig && methsig <: sig
    return methsig
end

"""
    callex = get_callexpr(sigex::ExLike)

Return the "call" expression for a signature-expression `sigex`.
(This strips out `:where` statements.)

# Example

```jldoctest; setup=:(using Revise)
julia> Revise.get_callexpr(:(nested(x::A) where A<:AbstractVector{T} where T))
:(nested(x::A))
```
"""
function get_callexpr(sigex::ExLike)
    while sigex.head == :where
        sigex = sigex.args[1]
    end
    sigex.head == :call || throw(ArgumentError(string("expected call expression, got ", sigex)))
    return sigex
end

"""
    typexs = sig_type_exprs(mod::Module, sigex::Expr)

From a function signature-expression `sigex` (see [`get_signature`](@ref)), generate a list
`typexs` of concrete signature type expressions.
This list will have length 1 unless `sigex` has default arguments,
in which case it will produce one type signature per valid number of supplied arguments.

These type-expressions can be evaluated in the appropriate module to obtain a Tuple-type.
The `mod` argument must be supplied if `sigex` has macrocalls in the arguments.

# Examples

```jldoctest; setup=:(using Revise)
julia> Revise.sig_type_exprs(Main, :(foo(x::Int, y::String)))
1-element Array{Expr,1}:
 :(Tuple{(Core.Typeof)(foo), Int, String})

julia> Revise.sig_type_exprs(Main, :(foo(x::Int, y::String="hello")))
2-element Array{Expr,1}:
 :(Tuple{(Core.Typeof)(foo), Int})
 :(Tuple{(Core.Typeof)(foo), Int, String})

julia> Revise.sig_type_exprs(Main, :(foo(x::AbstractVector{T}, y) where T))
1-element Array{Expr,1}:
 :(Tuple{(Core.Typeof)(foo), AbstractVector{T}, Any} where T)
```
"""
function sig_type_exprs(mod::Module, sigex::Expr, @nospecialize(wheres))
    if sigex.head == :(::)
        # return type annotation
        sigex = sigex.args[1]
    end
    if sigex.head == :where
        return sig_type_exprs(mod, sigex.args[1], (sigex.args[2:end], wheres...))
    end
    typexs = Expr[_sig_type_exprs(mod, sigex, wheres)]
    # If the method has default arguments, generate one type signature
    # for each valid call. This replicates the syntactic sugar that defines
    # multiple methods from a single definition.
    while has_default_args(sigex)
        sigex = Expr(sigex.head, sigex.args[1:end-1]...)
        push!(typexs, _sig_type_exprs(mod, sigex, wheres))
    end
    return reverse!(typexs)  # method table is organized in increasing # of args
end
sig_type_exprs(mod::Module, sigex::Expr) = sig_type_exprs(mod, sigex, ())
sig_type_exprs(mod::Module, sigex::RelocatableExpr) = sig_type_exprs(mod, convert(Expr, sigex))

function _sig_type_exprs(mod::Module, ex, @nospecialize(wheres))
    Typeof = Core.Typeof
    fex = ex.args[1]
    if isa(fex, Expr) && fex.head == :(::)
        fexTex = fex.args[end]
    else
        fexTex = :($Typeof($fex))
    end
    sigex = Expr(:curly, :Tuple, fexTex, argtypeexpr(mod, ex.args[2:end]...)...)
    for w in wheres
        sigex = Expr(:where, sigex, w...)
    end
    sigex
end

"""
    typeex1, typeex2, ... = argtypeexpr(mod, ex...)

Return expressions that specify the types assigned to each argument in a method signature.
Returns `:Any` if no type is assigned to a specific argument. It also skips
keyword arguments.

`ex...` should be arguments `2:end` of a `:call` expression (i.e., skipping over the
function name).

# Examples

```jldoctest; setup=:(using Revise), filter=r"#=.*=#"
julia> sigex = :(varargs(x, rest::Int...))
:(varargs(x, rest::Int...))

julia> Revise.argtypeexpr(Revise.get_callexpr(sigex).args[2:end]...)
(:Any, :(Vararg{Int}))

julia> sigex = :(complexargs(w::Vector{T}, @nospecialize(x::Integer), y, z::String=""; kwarg::Bool=false) where T)
:(complexargs(w::Vector{T}, #= REPL[39]:1 =# @nospecialize(x::Integer), y, z::String=""; kwarg::Bool=false) where T)

julia> Revise.argtypeexpr(Revise.get_callexpr(sigex).args[2:end]...)
(:(Vector{T}), :Integer, :Any, :String)
```
"""
function argtypeexpr(mod::Module, ex::ExLike, @nospecialize(rest...))
    isempty(ex.args) && return argtypeexpr(mod, rest...)  # issue #171
    # Handle @nospecialize(x)
    a = ex.args[1]
    if ex.head == :macrocall && length(ex.args) >= 3 && (a == Symbol("@nospecialize") ||
            (a isa Expr && a.head == :. && length(a.args) > 1 && a.args[end] == QuoteNode(Symbol("@nospecialize"))))
        return (argtypeexpr(mod, ex.args[3])..., argtypeexpr(mod, rest...)...)
    end
    if ex.head == :meta && length(ex.args) >= 2 && ex.args[1] == :nospecialize
        return (argtypeexpr(mod, ex.args[2])..., argtypeexpr(mod, rest...)...)
    end
    # Handle other macros (e.g., #176)
    if ex.head == :macrocall
        exex = macroexpand(mod, ex)
        return argtypeexpr(mod, exex, rest...)
    end
    if ex.head == :...
        # Handle varargs (those expressed with dots rather than Vararg{T,N})
        @assert isempty(rest)
        @assert length(ex.args) == 1
        T = argtypeexpr(mod, ex.args[1])[1]
        return (:(Vararg{$T}),)
    end
    # Skip over keyword arguments
    ex.head == :parameters && return argtypeexpr(mod, rest...)
    # Handle default arguments
    (ex.head == :(=) || ex.head == :kw) && return (argtypeexpr(mod, ex.args[1])..., argtypeexpr(mod, rest...)...)
    # Handle a destructured argument like foo(x, (count, name))
    ex.head == :tuple && return (:Any, argtypeexpr(mod, rest...)...)
    # Should be a type specification, check and then return the type
    ex.head == :(::) || throw(ArgumentError("expected :(::) expression, got $ex"))
    1 <= length(ex.args) <= 2 || throw(ArgumentError("expected 1 or 2 args, got $(ex.args)"))
    return (ex.args[end], argtypeexpr(mod, rest...)...)
end
argtypeexpr(mod::Module, s::Symbol, rest...) = (:Any, argtypeexpr(mod, rest...)...)
argtypeexpr(mod::Module) = ()
argtypeexpr(ex::Union{Symbol,ExLike}, rest...) = argtypeexpr(Main, ex, rest...)

function has_default_args(sigex::Expr)
    a = sigex.args[end]
    return isa(a, Expr) && a.head == :kw
end

function is_trivial_block_wrapper(ex::ExLike)
    if ex.head == :block
        return length(ex.args) == 1 ||
            (length(ex.args) == 2 && (is_linenumber(ex.args[1]) || ex.args[1]===nothing))
    end
    false
end
is_trivial_block_wrapper(@nospecialize arg) = false

function is_linenumber(@nospecialize stmt)
    isa(stmt, LineNumberNode) || (isa(stmt, ExLike) && (stmt.head == :line))
end

function firstlineno(rex::ExLike)
    for a in rex.args
        if is_linenumber(a)
            isa(a, LineNumberNode) && return a.line
            return a.args[1]
        end
        if isa(a, ExLike)
            lineno = firstlineno(a)
            isa(lineno, Integer) && return lineno
        end
    end
    return nothing
end


text:
export env!, env

@doc """
    env!(S::GphysData[, chans=CC, v=V])
    env(S::GphysData)

In-place conversion of S.x[i] ==> Env(S.x[i]) (≡ |H(S.x[i])|, where H denotes
the Hilbert transform).

### Keywords
* chans=CC: only process channels in CC (with fs > 0.0).
* v=V: verbosity.
""" env!
function env!(S::GphysData;
  chans::ChanSpec=Int64[],
  v::Integer=KW.v)

  # preprocess data channels
  chans = mkchans(chans, S, keepirr=false)

  proc_str = string("env!(S, chans=", chans, ")")

  # Get groups
  GRPS = get_unique(S, ["eltype"], chans)

  # Initialize Y
  Nx = nx_max(S, chans)
  (Nx == 0) && return
  Y = Array{Float64,1}(undef, 2*max(Nx, 128))

  # Arrays to hold data and windows
  @inbounds for grp in GRPS
    c = grp[1]
    T = eltype(S.x[c])
    R = reinterpret(Complex{T}, Y)

    # Get views and window lengths of each segment
    c2 = T(2)
    (L, X) = get_views(S, grp)
    nL = length(L)
    nx = L[1]

    n2 = div(nx, 2)
    H = view(R, 1:nx)
    v1 = view(R, 1:n2 + 1)
    v2 = view(R, 2:n2 + (isodd(nx) ? 1 : 0))
    v3 = view(R, n2+1:nx)
    P = plan_rfft(X[1])
    nx_last = nx
    np = 0

    for i = 1:nL
      nx = L[i]
      x = X[i]

      # determine whether to update P, YS
      too_short = false
      # update P, nx, YS
      if nx < 128 || nx != nx_last
        # very short segments
        if nx < 128
          y = zeros(eltype(x), 128)
          copyto!(y, 1, X[i], 1, nx)
          X₀ = y
          nx₀ = nx
          x = view(y, :)
          too_short = true
          nx = length(x)
        end
        P = plan_rfft(x)
        n2 = div(nx, 2)
        H = view(R, 1:nx)
        v1 = view(R, 1:n2 + 1)
        v2 = view(R, 2:n2 + (isodd(nx) ? 1 : 0))
        v3 = view(R, n2+1:nx)
      end

      #  Compute envelope → adapted from DSP.hilbert for recycled H
      fill!(v3, zero(Complex{T}))
      mul!(v1, P, x)
      broadcast!(*, v2, v2, c2)
      ifft!(H)

      # overwrite x
      if too_short
        copyto!(X[i], 1, x, 1, nx₀)
      else
        broadcast!(abs, x, H)
      end

      # update nx_last
      nx_last = nx

    end
    # log processing to :notes
    proc_note!(S, grp, proc_str, "replaced :x with abs(H(:x))")
  end
  return nothing
end

function env!(C::GphysChannel;
  v::Integer=KW.v
  )

  C.fs == 0.0 && return

  # Arrays to hold data and windows
  if size(C.t, 1) == 2
    C.x .= abs.(DSP.hilbert(C.x))
  else
    # Initialize R
    T = eltype(C.x)
    Nx = nx_max(C)
    (Nx == 0) && return
    R = Array{Complex{T}, 1}(undef, max(Nx, 128))

    # Get views and window lengths of each segment
    c2 = T(2)
    (L, X) = get_views(C)
    nL = length(L)
    nx = L[1]
    n2 = div(nx, 2)
    H = view(R, 1:nx)
    v1 = view(R, 1:n2 + 1)
    v2 = view(R, 2:n2 + (isodd(nx) ? 1 : 0))
    v3 = view(R, n2+1:nx)
    P = plan_rfft(X[1])
    nx_last = nx
    np = 0

    for i = 1:nL
      nx = L[i]
      x = X[i]

      # determine whether to update P, YS
      too_short = false
      # update P, nx, YS
      if nx < 128 || nx != nx_last
        # very short segments
        if nx < 128
          y = copyto!(zeros(eltype(x), 128), x)
          X₀ = y
          nx₀ = nx
          nx = 128
          x = view(y, :)
          too_short = true
        end
        P = plan_rfft(x)
        n2 = div(nx, 2)
        H = view(R, 1:nx)
        v1 = view(R, 1:n2 + 1)
        v2 = view(R, 2:n2 + (isodd(nx) ? 1 : 0))
        v3 = view(R, n2+1:nx)
      end

      #  Compute envelope → adapted from DSP.hilbert for recycled H
      fill!(v3, zero(Complex{T}))
      mul!(v1, P, x)
      broadcast!(*, v2, v2, c2)
      ifft!(H)

      # overwrite x
      if too_short
        copyto!(X[i], 1, x, 1, nx₀)
      else
        broadcast!(abs, x, H)
      end

      # update nx_last
      nx_last = nx
    end
  end

  # log processing to :notes
  proc_note!(C, "env!(C)", "replaced :x with abs(H(:x))")
  return nothing
end

@doc (@doc env!)
function env(S::GphysData;
  chans::ChanSpec=Int64[],
  v::Integer=KW.v
  )

  U = deepcopy(S)
  env!(U, chans=chans, v=v)
  return U
end

function env(C::GphysChannel;
  v::Integer=KW.v
  )

  U = deepcopy(C)
  env!(U, v=v)
  return U
end


text:
a = 10;
b = 2;
c = a * b;
println(c);

max_stars_repo_name: fahimfarhan/legendary-coding-odyssey

max_issues_repo_name: fahimfarhan/legendary-coding-odyssey

max_forks_repo_name: fahimfarhan/legendary-coding-odyssey


text:
using StanSample, DiffEqBayesStan, MCMCChains
using OrdinaryDiffEq
using ParameterizedFunctions
using RecursiveArrayTools
using Distributions
using StatsPlots

ProjDir = @__DIR__

println("Four parameter case")

f1 = @ode_def begin
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a b c d
u0 = [1.0,1.0]
tspan = (0.0,10.0)
p = [1.5,1.0,3.0,1.0]
prob1 = ODEProblem(f1,u0,tspan,p)
sol = solve(prob1,Tsit5())
t = collect(range(1,stop=10,length=50))
randomized = VectorOfArray([(sol(t[i]) + .3randn(2)) for i in 1:length(t)])
data_array = convert(Array, randomized)

priors = [truncated(Normal(1.5,0.1),0.5,2),truncated(Normal(1.0,0.1),0.5,1.5),
          truncated(Normal(3.0,0.1),0.5,4),truncated(Normal(1.0,0.1),0.5,2)]

bayesian_result = stan_inference(prob1, t, data_array, priors;
  nchains=[4], num_samples=1000, num_warmup=1000,
  vars =(DiffEqBayesStan.StanODEData(), InverseGamma(4,1)),
  output_format=:dataframe)

if success(bayesian_result.return_code)
  df = bayesian_result.chains

  println("\ntheta means:\n",
    [mean(df.theta_1), mean(df.theta_2), mean(df.theta_3), mean(df.theta_4)])
  println("\ntheta vars:\n",
    [var(df.theta_1), var(df.theta_2), var(df.theta_3), var(df.theta_4)])
  println()

  dfs = read_samples(bayesian_result.model; output_format=:dataframes)
  for j in 1:4
    for i in 1:4
      if i == 1
        plot(dfs[i][:, "theta_$j"], lab="chain[$j]")
      else
        plot!(dfs[i][:, "theta_$j"], lab="chain[$j]")
      end
    end
    savefig(joinpath(ProjDir, "tmp", "theta_$j.png"))
  end

  sdf = read_summary(bayesian_result.model)
  sdf |> display
  println("\n")

end


text:
import Base.show
# Main GaussianProcess type

@doc """
# Description
Fits a Gaussian process to a set of training points. The Gaussian process is defined in terms of its mean and covaiance (kernel) functions, which are user defined. As a default it is assumed that the observations are noise free.

# Constructors:
    GP(X, y, m, k, logNoise)
    GP(; m=MeanZero(), k=SE(0.0, 0.0), logNoise=-1e8) # observation-free constructor

# Arguments:
* `X::Matrix{Float64}`: Input observations
* `y::Vector{Float64}`: Output observations
* `m::Mean`           : Mean function
* `k::kernel`         : Covariance function
* `logNoise::Float64` : Log of the standard deviation for the observation noise. The default is -1e8, which is equivalent to assuming no observation noise.

# Returns:
* `gp::GP`            : Gaussian process object, fitted to the training data if provided
""" ->
type GP
    m:: Mean                # Mean object
    k::Kernel               # Kernel object
    logNoise::Float64       # log standard deviation of observation noise

    # Observation data
    nobsv::Int              # Number of observations
    X::Matrix{Float64}      # Input observations
    y::Vector{Float64}      # Output observations
    data::KernelData        # Auxiliary observation data (to speed up calculations)
    dim::Int                # Dimension of inputs

    # Auxiliary data
    cK::AbstractPDMat       # (k + exp(2*obsNoise))
    alpha::Vector{Float64}  # (k + exp(2*obsNoise))⁻¹y
    mLL::Float64            # Marginal log-likelihood
    dmLL::Vector{Float64}   # Gradient marginal log-likelihood

    function GP(X::Matrix{Float64}, y::Vector{Float64}, m::Mean, k::Kernel, logNoise::Float64=-1e8)
        dim, nobsv = size(X)
        length(y) == nobsv || throw(ArgumentError("Input and output observations must have consistent dimensions."))
        gp = new(m, k, logNoise, nobsv, X, y, KernelData(k, X), dim)
        update_mll!(gp)
        return gp
    end

    GP(; m=MeanZero(), k=SE(0.0, 0.0), logNoise=-1e8) =  new(m, k, logNoise, 0)

end

# Creates GP object for 1D case
GP(x::Vector{Float64}, y::Vector{Float64}, meanf::Mean, kernel::Kernel, logNoise::Float64=-1e8) = GP(x', y, meanf, kernel, logNoise)

@doc """
# Description
Fits an existing Gaussian process to a set of training points.

# Arguments:
* `gp::GP`: Exiting Gaussian process object
* `X::Matrix{Float64}`: Input observations
* `y::Vector{Float64}`: Output observations

# Returns:
* `gp::GP`            : A Gaussian process fitted to the training data
""" ->
function fit!(gp::GP, X::Matrix{Float64}, y::Vector{Float64})
    length(y) == size(X,2) || throw(ArgumentError("Input and output observations must have consistent dimensions."))
    gp.X = X
    gp.y = y
    gp.data = KernelData(gp.k, X)
    gp.dim, gp.nobsv = size(X)
    update_mll!(gp)
    return gp
end

fit!(gp::GP, x::Vector{Float64}, y::Vector{Float64}) = fit!(gp, x', y)

# Update auxiliarly data in GP object after changes have been made
function update_mll!(gp::GP)
    μ = mean(gp.m,gp.X)
    Σ = cov(gp.k, gp.X, gp.data)
    gp.cK = PDMat(Σ + exp(2*gp.logNoise)*eye(gp.nobsv) + 1e-8*eye(gp.nobsv))
    gp.alpha = gp.cK \ (gp.y - μ)
    gp.mLL = -dot((gp.y - μ),gp.alpha)/2.0 - logdet(gp.cK)/2.0 - gp.nobsv*log(2π)/2.0 # Marginal log-likelihood
end

# Update gradient of marginal log likelihood

function update_mll_and_dmll!(gp::GP; noise::Bool=true, mean::Bool=true, kern::Bool=true)
    update_mll!(gp::GP)
    gp.dmLL = Array(Float64, noise + mean*num_params(gp.m) + kern*num_params(gp.k))

    # Calculate Gradient with respect to hyperparameters

    #Derivative wrt the observation noise
    if noise
        gp.dmLL[1] = exp(2*gp.logNoise)*trace((gp.alpha*gp.alpha' - gp.cK \ eye(gp.nobsv)))
    end

    #Derivative wrt to mean hyperparameters, need to loop over as same with kernel hyperparameters
    if mean
        Mgrads = grad_stack(gp.m, gp.X)
        for i in 1:num_params(gp.m)
            gp.dmLL[i+noise] = dot(Mgrads[:,i],gp.alpha)
        end
    end

    # Derivative of marginal log-likelihood with respect to kernel hyperparameters
    if kern
        Kgrads = grad_stack(gp.k, gp.X, gp.data)   # [dK/dθᵢ]
        for i in 1:num_params(gp.k)
            gp.dmLL[i+mean*num_params(gp.m)+noise] = trace((gp.alpha*gp.alpha' - gp.cK \ eye(gp.nobsv))*Kgrads[:,:,i])/2
        end
    end
end

@doc """
# Description
Calculates the posterior mean and variance of Gaussian Process at specified points

# Arguments:
* `gp::GP`: Gaussian Process object
* `X::Matrix{Float64}`:  matrix of points for which one would would like to predict the value of the process.
                       (each column of the matrix is a point)

# Keyword Arguments
* `full_cov::Bool`: indicates whether full covariance matrix should be returned instead of only variances (default is false)

# Returns:
* `(mu, Sigma)::(Vector{Float64}, Vector{Float64})`: respectively the posterior mean  and variances of the posterior
                                                    process at the specified points
""" ->
function predict(gp::GP, x::Matrix{Float64}; full_cov::Bool=false)
    size(x,1) == gp.dim || throw(ArgumentError("Gaussian Process object and input observations do not have consistent dimensions"))
    if full_cov
        return _predict(gp, x)
    else
        ## Calculate prediction for each point independently
            μ = Array(Float64, size(x,2))
            σ2 = similar(μ)
        for k in 1:size(x,2)
            m, sig = _predict(gp, x[:,k:k])
            μ[k] = m[1]
            σ2[k] = max(full(sig)[1,1], 0.0)
        end
        return μ, σ2
    end
end

# 1D Case for prediction
predict(gp::GP, x::Vector{Float64}; full_cov::Bool=false) = predict(gp, x'; full_cov=full_cov)

## compute predictions
function _predict(gp::GP, X::Array{Float64})
    n = size(X, 2)
    cK = cov(gp.k, X, gp.X)
    Lck = whiten(gp.cK, cK')
    mu = mean(gp.m,X) + cK*gp.alpha        # Predictive mean
    Sigma_raw = cov(gp.k, X) - Lck'Lck # Predictive covariance
    # Hack to get stable covariance
    Sigma = try PDMat(Sigma_raw) catch; PDMat(Sigma_raw+1e-8*sum(diag(Sigma_raw))/n*eye(n)) end
    return (mu, Sigma)
end

# Sample from the GP
function rand!(gp::GP, X::Matrix{Float64}, A::DenseMatrix)
    nobsv = size(X,2)
    n_sample = size(A,2)

    if gp.nobsv == 0
        # Prior mean and covariance
        μ = mean(gp.m, X);
        Σraw = cov(gp.k, X);
        Σ = try PDMat(Σraw) catch; PDMat(Σraw+1e-8*sum(diag(Σraw))/nobsv*eye(nobsv)) end
    else
        # Posterior mean and covariance
        μ, Σ = predict(gp, X; full_cov=true)
    end

    return broadcast!(+, A, μ, unwhiten!(Σ,randn(nobsv, n_sample)))
end

function rand(gp::GP, X::Matrix{Float64}, n::Int)
    nobsv=size(X,2)
    A = Array(Float64, nobsv, n)
    return rand!(gp, X, A)
end

# Sample from 1D GP
rand(gp::GP, x::Vector{Float64}, n::Int) = rand(gp, x', n)

# Generate only one sample from the GP and returns a vector
rand(gp::GP,X::Matrix{Float64}) = vec(rand(gp,X,1))
rand(gp::GP,x::Vector{Float64}) = vec(rand(gp,x',1))

function get_params(gp::GP; noise::Bool=true, mean::Bool=true, kern::Bool=true)
    params = Float64[]
    if noise; push!(params, gp.logNoise); end
    if mean;  append!(params, get_params(gp.m)); end
    if kern; append!(params, get_params(gp.k)); end
    return params
end

function set_params!(gp::GP, hyp::Vector{Float64}; noise::Bool=true, mean::Bool=true, kern::Bool=true)
    # println("mean=$(mean)")
    if noise; gp.logNoise = hyp[1]; end
    if mean; set_params!(gp.m, hyp[1+noise:noise+num_params(gp.m)]); end
    if kern; set_params!(gp.k, hyp[end-num_params(gp.k)+1:end]); end
end

function push!(gp::GP, X::Matrix{Float64}, y::Vector{Float64})
    warn("push! method is currently inefficient as it refits all observations")
    if gp.nobsv == 0
        GaussianProcesses.fit!(gp, X, y)
    elseif size(X,1) != size(gp.X,1)
        error("New input observations must have dimensions consistent with existing observations")
    else
        GaussianProcesses.fit!(gp, cat(2, gp.X, X), cat(1, gp.y, y))
    end
end

push!(gp::GP, x::Vector{Float64}, y::Vector{Float64}) = push!(gp, x', y)
push!(gp::GP, x::Float64, y::Float64) = push!(gp, [x], [y])
push!(gp::GP, x::Vector{Float64}, y::Float64) = push!(gp, reshape(x, length(x), 1), [y])

function show(io::IO, gp::GP)
    println(io, "GP object:")
    println(io, "  Dim = $(gp.dim)")
    println(io, "  Number of observations = $(gp.nobsv)")
    println(io, "  Mean function:")
    show(io, gp.m, 2)
    println(io, "  Kernel:")
    show(io, gp.k, 2)
    if (gp.nobsv == 0)
        println("  No observation data")
    else
        println(io, "  Input observations = ")
        show(io, gp.X)
        print(io,"\n  Output observations = ")
        show(io, gp.y)
        print(io,"\n  Variance of observation noise = $(exp(2*gp.logNoise))")
        print(io,"\n  Marginal Log-Likelihood = ")
        show(io, round(gp.mLL,3))
    end
end


text:
using FLOWMath
using Test

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


text:
f = open("jawiki-country.txt", "r")

for line in readlines(f)
    m = match(r"(^=+)", line)
    if !is(m, nothing)
     level = length(m.captures[1]) - 1
     println(replace(line, r"=|\s", "") * " $level")
   end
end


text:
#!/usr/bin/env julia

function write_toml(f::Function, names::String...)
    path = joinpath(names...) * ".toml"
    mkpath(dirname(path))
    open(path, "w") do io
        f(io)
    end
end

toml_key(str::String) = ismatch(r"[^\w-]", str) ? repr(str) : str
toml_key(strs::String...) = join(map(toml_key, [strs...]), '.')

prefix = joinpath(homedir(), ".julia", "registries", "Uncurated")

write_toml(prefix, "registry") do io
    repo = "https://github.com/JuliaRegistries/Uncurated.git"
    uuid = string(uuid5(uuid_registry, repo))
    println(io, "name = ", repr("Uncurated"))
    println(io, "uuid = ", repr(uuid))
    println(io, "repo = ", repr(repo))
    println(io, "\ndescription = \"\"\"")
    print(io, """
        Official uncurated Julia package registry where people can
        register any package they want without too much debate about
        naming and without enforced standards on documentation or
        testing. We nevertheless encourage documentation, testing and
        some amount of consideration when choosing package names.
        """)
    println(io, "\"\"\"")
    println(io, "\n[packages]")
    for (pkg, p) in sort!(collect(pkgs), by=(p->p.uuid.value)∘last)
        bucket = string(uppercase(first(pkg)))
        path = joinpath(bucket, pkg)
        println(io, p.uuid, " = { name = ", repr(pkg), ", path = ", repr(path), " }")
    end
end

buckets = Dict()
for (pkg, p) in pkgs
    bucket = string(uppercase(first(pkg)))
    push!(get!(buckets, bucket, []), (pkg, p))
end

include("utils.jl")
include("sha1map.jl")
const trees = sha1map(pkgs)

for (bucket, b_pkgs) in buckets, (pkg, p) in b_pkgs
    url = p.url
    uuid = string(p.uuid)
    startswith(url, "git://github.com") && (url = "https"*url[4:end])

    # package.toml
    write_toml(prefix, bucket, pkg, "package") do io
        println(io, "name = ", repr(pkg))
        println(io, "uuid = ", repr(uuid))
        println(io, "repo = ", repr(url))
    end

    # versions.toml
    write_toml(prefix, bucket, pkg, "versions") do io
        for (i, (ver, v)) in enumerate(sort!(collect(p.versions), by=first))
            i > 1 && println(io)
            println(io, "[", toml_key(string(ver)), "]")
            println(io, "hash-sha1 = ", repr(trees[uuid][v.sha1]))
        end
    end
    versions = sort!(collect(keys(p.versions)))

    function write_versions_data(f::Function, name::String; lt::Function=isless)
        data = Dict{VersionNumber,Dict{String,String}}()
        for (ver, v) in p.versions, (dep, d) in v.requires
            val = f(dep, d)
            val == nothing && continue
            haskey(data, ver) || (data[ver] = Dict{String,String}())
            data[ver][dep] = val
        end
        compressed = compress_versions_data(data, versions)
        !isempty(compressed) && write_toml(prefix, bucket, pkg, name) do io
            vers = unique(getindex.(compressed, 1))
            keys = sort!(unique(getindex.(compressed, 2)), lt=lt)
            what = (vers, keys)
            ord = (1, 2)
            for (i, x) in enumerate(what[ord[1]])
                i > 1 && println(io)
                println(io, "[", toml_key(x), "]")
                for y in what[ord[2]]
                    for t in compressed
                        t[ord[1]] == x && t[ord[2]] == y || continue
                        println(io, toml_key(y), " = ", t[3])
                    end
                end
            end
        end
    end

    # dependencies.toml
    write_versions_data("dependencies") do dep, d
        dep != "julia" ? repr(string(pkgs[dep].uuid)) : nothing
    end

    # compatibility.toml
    write_versions_data("compatibility", lt=packagelt) do dep, d
        versions_repr(compress_versions(
            d.versions, collect(keys(pkgs[dep].versions))
        ))
    end
end