text:
"""
    translation(C::SphericalHarmonicCoefficients,v::Array{Float64,1},x::Variable, y::Variable, z::Variable)
*Description:* Translation of the coefficients: Shifting the expansion point by v\\
\\
*Input:*
- `C`  - Coefficients
- `v`  - shift vector (length(v) = 3)
- `x, y, z` - Cartesian coordinates
*Output:*
- `cShifted` - Shifted coefficients, type: SphericalHarmonicCoefficients (cShifted.R = C.R, cShifted.solid = C.solid)

"""
function translation(C::SphericalHarmonicCoefficients,v::Array{Float64,1},x::Variable, y::Variable, z::Variable)

    cShifted = deepcopy(C)

    vx = v[1]
    vy = v[2]
    vz = v[3]

    # turn solid coefficients into spherical coefficients
    if C.solid
        spherical!(C)
    end

    # c[0,0]:
    sum = 0
    # Summand aus (I) - m > 0
    for λ=1:C.L
        for μ = 1:λ
            p = rlm(λ,μ,x,y,z)
            sum += ω₊(C[λ,μ],λ,μ,0,0,1)*p((x,y,z)=>(vx,vy,vz))
        end
    end
    # Summand aus (II) - m = 0
    for λ=0:C.L
        p = rlm(λ,0,x,y,z)
        sum += ω₀(C[λ,0],λ,0,0,1)*p((x,y,z)=>(vx,vy,vz))
    end
    # Summand aus (III) - m < 0
    for λ=1:C.L
        for μ = -λ:-1
            p = rlm(λ,μ,x,y,z)
            sum += ω₋(C[λ,μ],λ,μ,0,0,1)*p((x,y,z)=>(vx,vy,vz))
        end
    end
    cShifted[0,0] = sqrt(4*pi)*sum
    sum = 0

    # cShifted[l,0] (l > 0):
    for l=1:C.L
        # Summand aus (I) - m > 0
        for λ=l:C.L
            for μ = 1:λ-l
                p = rlm(λ-l,μ,x,y,z)
                sum += ω₊(C[λ,μ],λ,μ,l,0,1)*p((x,y,z)=>(vx,vy,vz))
            end
        end
        # Summand aus (II) - m = 0
        for λ=l:C.L
            p = rlm(λ-l,0,x,y,z)
            sum += ω₀(C[λ,0],λ,l,0,1)*p((x,y,z)=>(vx,vy,vz))
        end
        # Summand aus (III) - m < 0
        for λ=l:C.L
            for μ = -(λ-l):-1
                p = rlm(λ-l,μ,x,y,z)
                sum += ω₋(C[λ,μ],λ,μ,l,0,1)*p((x,y,z)=>(vx,vy,vz))
            end
        end
        cShifted[l,0] = sqrt(4*pi/(2*l+1))*sum
        sum = 0
    end

    # cShifted[l,m] for l ≠ 0, m > 0
    for l=1:C.L
        for m=1:l
            # 1: Summand 1 aus (I) - m > 0
            for λ=l:C.L
                for μ = m:m-(l-λ)
                    p = rlm(λ-l,μ-m,x,y,z)
                    sum += ω₊(C[λ,μ],λ,μ,l,m,1)*p((x,y,z)=>(vx,vy,vz))
                end
            end
            # 2: Summand 2 aus (I) - m > 0
            for λ=l+1:C.L
                for μ = max(1,m-(λ-l)):m-1
                    p = rlm(λ-l,abs(μ-m),x,y,z)
                    sum += ω₊(C[λ,μ],λ,μ,l,m,2)*p((x,y,z)=>(vx,vy,vz))
                end
            end
            # 3: Summand 3 aus (I) - m > 0
            for λ=l+1:C.L
                for μ = 1:-m-(l-λ)
                    p = rlm(λ-l,μ+m,x,y,z)
                    sum += ω₊(C[λ,μ],λ,μ,l,-m,3)*p((x,y,z)=>(vx,vy,vz))
                end
            end
            # 4: Summand aus (II) - m = 0
            for λ=m+l:C.L
                p = rlm(λ-l,m,x,y,z)
                sum += (ω₀(C[λ,0],λ,l,m,2)+ω₀(C[λ,0],λ,l,-m,3))*p((x,y,z)=>(vx,vy,vz))
            end
            # 5: Summand 1 aus (III) - m < 0
            for λ=l:C.L
                for μ = -m-(λ-l):-m-1 #min(-1,-m)
                    p = rlm(λ-l,μ+m,x,y,z)
                    sum += ω₋(C[λ,μ],λ,μ,l,-m,1)*p((x,y,z)=>(vx,vy,vz))
                end
            end
            # 6: Summand 2 aus (III) - m < 0
            for λ=l+1:C.L
                for μ = -m+1:min(-1,-m-(l-λ))
                    p = rlm(λ-l,-(μ+m),x,y,z)
                    sum -= ω₋(C[λ,μ],λ,μ,l,-m,2)*p((x,y,z)=>(vx,vy,vz))
                end
            end
            # 7: Summand 3 aus (III) - m < 0
            for λ=l+1:C.L
                for μ = m-(λ-l):-1
                    p = rlm(λ-l,μ-m,x,y,z)
                    sum += ω₋(C[λ,μ],λ,μ,l,m,3)*p((x,y,z)=>(vx,vy,vz))
                end
            end
            cShifted[l,m] = sqrt(4*pi/(2*l+1))*sum
            sum = 0
        end
    end

    # cShifted[l,m] for l ≠ 0, m < 0
    for l=1:C.L
        for m=-l:-1
            # 1: Summand 1 aus (I) - m > 0
            for λ=l:C.L
                for μ = -m+1:-m-(l-λ)
                    p = rlm(λ-l,-(μ+m),x,y,z)
                    sum -= ω₊(C[λ,μ],λ,μ,l,-m,1)*p((x,y,z)=>(vx,vy,vz))
                end
            end
            # 2: Summand 2 aus (I) - m > 0
            for λ=l+1:C.L
                for μ = max(1,-m-(λ-l)):-m-1
                    p = rlm(λ-l,-abs(μ+m),x,y,z)
                    sum += ω₊(C[λ,μ],λ,μ,l,-m,2)*p((x,y,z)=>(vx,vy,vz))
                end
            end
            # 3: Summand 3 aus (I) - m > 0
            for λ=l+1:C.L
                for μ = 1:m-(l-λ)
                    p = rlm(λ-l,-(μ-m),x,y,z)
                    sum += ω₊(C[λ,μ],λ,μ,l,m,3)*p((x,y,z)=>(vx,vy,vz))
                end
            end
            # 4: Summand aus (II) - m = 0
            for λ=l-m:C.L
                p = rlm(λ-l,m,x,y,z)
                sum += (ω₀(C[λ,0],λ,l,-m,2)+ω₀(C[λ,0],λ,l,m,3))*p((x,y,z)=>(vx,vy,vz))
            end
            # 5: Summand 1 aus (III) - m < 0
            for λ=l:C.L
                for μ = m-(λ-l):m#-1 #min(-1,m)
                    p = rlm(λ-l,abs(μ-m),x,y,z)
                    sum += ω₋(C[λ,μ],λ,μ,l,m,1)*p((x,y,z)=>(vx,vy,vz))
                end
            end
            # 6: Summand 2 aus (III) - m < 0
            for λ=l+1:C.L
                for μ = m+1:min(-1,m-(l-λ))
                    p = rlm(λ-l,μ-m,x,y,z)
                    sum += ω₋(C[λ,μ],λ,μ,l,m,2)*p((x,y,z)=>(vx,vy,vz))
                end
            end
            # 7: Summand 3 aus (III) - m < 0
            for λ=l+1:C.L
                for μ = -m-(λ-l):-1
                    p = rlm(λ-l,-(μ+m),x,y,z)
                    sum -= ω₋(C[λ,μ],λ,μ,l,-m,3)*p((x,y,z)=>(vx,vy,vz))
                end
            end
            cShifted[l,m] = sqrt(4*pi/(2*l+1))*sum
            sum = 0
        end
    end

    # turn spherical coefficients into solid coefficients
    if cShifted.solid
        cShifted.solid = false
        solid!(cShifted)
        solid!(C)
    end

    return cShifted
end

function ω₀(c,λ,l,m,num)
    ω = c*factorial(λ)*sqrt((2*λ+1)/(4*pi))/sqrt(4*factorial(l+m)*factorial(l-m)*factorial(λ-l+m)*factorial(λ-l-m))
    if num == 1 # m == 0
        ω *= 2
    else # num == 2 || num == 3
        ω *= (-1)^m
    end
    return ω
end

function ω₊(c,λ,μ,l,m,num)
    ω = c*sqrt(2*factorial(λ+μ)*factorial(λ-μ))*sqrt((2*λ+1)/(4*pi))/sqrt(4*factorial(l+m)*factorial(l-m)*factorial(λ-l+μ-m)*factorial(λ-l-(μ-m)))
    if num == 1
        if m == μ || m == 0
            ω *= sqrt(2)
        else
            ω *= 1
        end
    elseif num == 2
        ω *= (-1)^(μ-m)
    else # num == 3
        ω *= (-1)^m
    end
    return ω
end

function ω₋(c,λ,μ,l,m,num)
    ω = c*(-1)^μ*sqrt(2*factorial(λ+μ)*factorial(λ-μ))*sqrt((2*λ+1)/(4*pi))/sqrt(4*factorial(l+m)*factorial(l-m)*factorial(λ-l+μ-m)*factorial(λ-l-(μ-m)))
    if num == 1
        if m == μ || m == 0
            ω *= (-1)^μ*sqrt(2)
        else
            ω *= (-1)^μ
        end
    elseif num == 2
        ω *= (-1)^m
    else # num == 3
        ω *= (-1)^(μ-m)
    end
    return ω
end

function translateRlm(l::Int64, m::Int64,vx,vy,vz, x::Variable, y::Variable, z::Variable)

    Rlmt = 0;

    # two cases m≥0 and m<0, each with three summands:
    if m >= 0
    for λ=0:l
        for μ=max(0,λ-l+m):min(λ,m)
            if m == 0
                p = rlm(l-λ,0,x,y,z);
                mult = 1/(factorial(λ)*factorial(l-λ));
                sum1 = rlm(λ,0,x,y,z)*p((x,y,z)=>(vx,vy,vz));
                Rlmt += mult*sum1;
            elseif μ == 0 # => m-μ = m > 0
                p = rlm(l-λ,m,x,y,z);
                mult = 1/(factorial(λ)*sqrt(2*factorial(l-λ+m)*factorial(l-λ-m)));
                sum1 = rlm(λ,0,x,y,z)*p((x,y,z)=>(vx,vy,vz));
                Rlmt += mult*sum1;
            elseif μ == m # => m-μ = 0
                p = rlm(l-λ,0,x,y,z);
                mult = 1/(factorial(l-λ)*sqrt(2*factorial(λ+m)factorial(λ-m)));
                sum1 = rlm(λ,m,x,y,z)*p((x,y,z)=>(vx,vy,vz));
                Rlmt += mult*sum1;
            else # m-μ > 0, μ > 0
                p = rlm(l-λ,m-μ,x,y,z);
                q = rlm(l-λ,-(m-μ),x,y,z);
                mult = 1/sqrt(4*factorial(λ+μ)*factorial(λ-μ)*factorial(l-λ+m-μ)*factorial(l-λ-(m-μ)));
                sum1 = rlm(λ,μ,x,y,z)*p((x,y,z)=>(vx,vy,vz));
                sum2 = rlm(λ,-μ,x,y,z)*q((x,y,z)=>(vx,vy,vz));
                Rlmt += mult*(sum1 - sum2);
            end
        end
    end

    for λ=m+1:l-1
        for μ=m+1:min(λ,-λ+l+m) # μ > 0, m-μ < 0
            p = rlm(l-λ,abs(m-μ),x,y,z);
            q = rlm(l-λ,-abs(m-μ),x,y,z);
            mult = (-1)^(m-μ)/sqrt(4*factorial(λ+μ)*factorial(λ-μ)*factorial(l-λ+(m-μ))*factorial(l-λ-(m-μ)));
            sum1 = rlm(λ,μ,x,y,z)*p((x,y,z)=>(vx,vy,vz));
            sum2 = rlm(λ,-μ,x,y,z)*q((x,y,z)=>(vx,vy,vz));
            Rlmt += mult*(sum1 + sum2);
        end
    end

    for λ=1:l-m-1
        for μ=max(-λ,λ-l+m):-1 # μ < 0, m-μ > 0
            p = rlm(l-λ,m-μ,x,y,z);
            q = rlm(l-λ,-(m-μ),x,y,z);
            mult = (-1)^(μ)/sqrt(4*factorial(λ+μ)*factorial(λ-μ)*factorial(l-λ+m-μ)*factorial(l-λ-(m-μ)));
            sum1 = rlm(λ,abs(μ),x,y,z)*p((x,y,z)=>(vx,vy,vz));
            sum2 = rlm(λ,-abs(μ),x,y,z)*q((x,y,z)=>(vx,vy,vz));
            Rlmt += mult*(sum1 + sum2);
        end
    end

    Rlmt *= sqrt(factorial(l-m)*factorial(l+m));
    Rlmt *= (m == 0) ? 1 : sqrt(2);

    else # m < 0
        for λ=0:l
            for μ=max(m,-λ):min(0,-λ+l+m)
                if μ == 0 # => m-μ = m < 0
                    p = rlm(l-λ,-abs(m),x,y,z);
                    mult = (-1)^m/(factorial(λ)*sqrt(2*factorial(l-λ+m)*factorial(l-λ-m)));
                    sum1 = rlm(λ,0,x,y,z)*p((x,y,z)=>(vx,vy,vz));
                    Rlmt += (-1)*mult*sum1;
                elseif μ == m # => m-μ = 0
                    p = rlm(l-λ,0,x,y,z);
                    mult = (-1)^m/(factorial(l-λ)*sqrt(2*factorial(λ+m)factorial(λ-m)));
                    sum1 = rlm(λ,-abs(m),x,y,z)*p((x,y,z)=>(vx,vy,vz));
                    Rlmt += (-1)*mult*sum1;
                else # μ < 0, m-μ < 0
                    p = rlm(l-λ,abs(m-μ),x,y,z);
                    q = rlm(l-λ,-abs(m-μ),x,y,z);
                    mult = (-1)^m/sqrt(4*factorial(λ+μ)*factorial(λ-μ)*factorial(l-λ+(m-μ))*factorial(l-λ-(m-μ)));
                    sum1 = rlm(λ,-abs(μ),x,y,z)*p((x,y,z)=>(vx,vy,vz));
                    sum2 = rlm(λ,abs(μ),x,y,z)*q((x,y,z)=>(vx,vy,vz));
                    Rlmt += (-1)*mult*(sum1 + sum2);
                end
            end
        end

        for λ=-m+1:l-1
            for μ=max(-λ,λ-l+m):m-1 # μ < 0, m-μ > 0
                p = rlm(l-λ,-(m-μ),x,y,z);
                q = rlm(l-λ,(m-μ),x,y,z);
                mult = (-1)^(μ)/sqrt(4*factorial(λ+μ)*factorial(λ-μ)*factorial(l-λ+m-μ)*factorial(l-λ-(m-μ)));
                sum1 = rlm(λ,abs(μ),x,y,z)*p((x,y,z)=>(vx,vy,vz));
                sum2 = rlm(λ,-abs(μ),x,y,z)*q((x,y,z)=>(vx,vy,vz));
                Rlmt += mult*(sum1 - sum2);
            end
        end

        for λ=1:l+m-1
            for μ=1:min(λ,-λ+l+m) # μ > 0, m-μ < 0
                p = rlm(l-λ,abs(m-μ),x,y,z);
                q = rlm(l-λ,-abs(m-μ),x,y,z);
                mult = (-1)^(m-μ)/sqrt(4*factorial(λ+μ)*factorial(λ-μ)*factorial(l-λ+(m-μ))*factorial(l-λ-(m-μ)));
                sum1 = rlm(λ,-μ,x,y,z)*p((x,y,z)=>(vx,vy,vz));
                sum2 = rlm(λ,μ,x,y,z)*q((x,y,z)=>(vx,vy,vz));
                Rlmt += mult*(sum1 - sum2);
            end
        end

        Rlmt *= (-1)*(-1)^m*sqrt(2*factorial(l-abs(m))*factorial(l+abs(m)));
    end
    return Rlmt
end

# Error propagation
function errorTranslation(C::SphericalHarmonicCoefficients,v,x::Variable, y::Variable, z::Variable)

    cShifted = deepcopy(C)

    vx = v[1]
    vy = v[2]
    vz = v[3]

    # turn solid coefficients into spherical coefficients
    if C.solid
        spherical!(C)
    end

    # cShifted[0,0]:
    sum = 0
    # Summand aus (I) - m > 0
    for λ=1:C.L
        for μ = 1:λ
            p = rlm(λ,μ,x,y,z)
            sum += abs(ω₊(C[λ,μ],λ,μ,0,0,1)*p((x,y,z)=>(vx,vy,vz)))
        end
    end
    # Summand aus (II) - m = 0
    for λ=0:C.L
        p = rlm(λ,0,x,y,z)
        sum += abs(ω₀(C[λ,0],λ,0,0,1)*p((x,y,z)=>(vx,vy,vz)))
    end
    # Summand aus (III) - m < 0
    for λ=1:C.L
        for μ = -λ:-1
            p = rlm(λ,μ,x,y,z)
            sum += abs(ω₋(C[λ,μ],λ,μ,0,0,1)*p((x,y,z)=>(vx,vy,vz)))
        end
    end
    cShifted[0,0] = sqrt(4*pi)*sum
    sum = 0

    # cShifted[l,0] (l > 0):
    for l=1:C.L
        # Summand aus (I) - m > 0
        for λ=l:C.L
            for μ = 1:λ-l
                p = rlm(λ-l,μ,x,y,z)
                sum += abs(ω₊(C[λ,μ],λ,μ,l,0,1)*p((x,y,z)=>(vx,vy,vz)))
            end
        end
        # Summand aus (II) - m = 0
        for λ=l:C.L
            p = rlm(λ-l,0,x,y,z)
            sum += abs(ω₀(C[λ,0],λ,l,0,1)*p((x,y,z)=>(vx,vy,vz)))
        end
        # Summand aus (III) - m < 0
        for λ=l:C.L
            for μ = -(λ-l):-1
                p = rlm(λ-l,μ,x,y,z)
                sum += abs(ω₋(C[λ,μ],λ,μ,l,0,1)*p((x,y,z)=>(vx,vy,vz)))
            end
        end
        cShifted[l,0] = sqrt(4*pi/(2*l+1))*sum
        sum = 0
    end

    # cShifted[l,m] for l ≠ 0, m > 0
    for l=1:C.L
        for m=1:l
            # 1: Summand 1 aus (I) - m > 0
            for λ=l:C.L
                for μ = m:m-(l-λ)
                    p = rlm(λ-l,μ-m,x,y,z)
                    sum += abs(ω₊(C[λ,μ],λ,μ,l,m,1)*p((x,y,z)=>(vx,vy,vz)))
                end
            end
            # 2: Summand 2 aus (I) - m > 0
            for λ=l+1:C.L
                for μ = max(1,m-(λ-l)):m-1
                    p = rlm(λ-l,abs(μ-m),x,y,z)
                    sum += abs(ω₊(C[λ,μ],λ,μ,l,m,2)*p((x,y,z)=>(vx,vy,vz)))
                end
            end
            # 3: Summand 3 aus (I) - m > 0
            for λ=l+1:C.L
                for μ = 1:-m-(l-λ)
                    p = rlm(λ-l,μ+m,x,y,z)
                    sum += abs(ω₊(C[λ,μ],λ,μ,l,-m,3)*p((x,y,z)=>(vx,vy,vz)))
                end
            end
            # 4: Summand aus (II) - m = 0
            for λ=m+l:C.L
                p = rlm(λ-l,m,x,y,z)
                sum += abs((ω₀(C[λ,0],λ,l,m,2)+ω₀(C[λ,0],λ,l,-m,3))*p((x,y,z)=>(vx,vy,vz)))
            end
            # 5: Summand 1 aus (III) - m < 0
            for λ=l:C.L
                for μ = -m-(λ-l):-m-1 #min(-1,-m)
                    p = rlm(λ-l,μ+m,x,y,z)
                    sum += abs(ω₋(C[λ,μ],λ,μ,l,-m,1)*p((x,y,z)=>(vx,vy,vz)))
                end
            end
            # 6: Summand 2 aus (III) - m < 0
            for λ=l+1:C.L
                for μ = -m+1:min(-1,-m-(l-λ))
                    p = rlm(λ-l,-(μ+m),x,y,z)
                    sum += abs(ω₋(C[λ,μ],λ,μ,l,-m,2)*p((x,y,z)=>(vx,vy,vz)))
                end
            end
            # 7: Summand 3 aus (III) - m < 0
            for λ=l+1:C.L
                for μ = m-(λ-l):-1
                    p = rlm(λ-l,μ-m,x,y,z)
                    sum += abs(ω₋(C[λ,μ],λ,μ,l,m,3)*p((x,y,z)=>(vx,vy,vz)))
                end
            end
            cShifted[l,m] = sqrt(4*pi/(2*l+1))*sum
            sum = 0
        end
    end

    # cShifted[l,m] for l ≠ 0, m < 0
    for l=1:C.L
        for m=-l:-1
            # 1: Summand 1 aus (I) - m > 0
            for λ=l:C.L
                for μ = -m+1:-m-(l-λ)
                    p = rlm(λ-l,-(μ+m),x,y,z)
                    sum += abs(ω₊(C[λ,μ],λ,μ,l,-m,1)*p((x,y,z)=>(vx,vy,vz)))
                end
            end
            # 2: Summand 2 aus (I) - m > 0
            for λ=l+1:C.L
                for μ = max(1,-m-(λ-l)):-m-1
                    p = rlm(λ-l,-abs(μ+m),x,y,z)
                    sum += abs(ω₊(C[λ,μ],λ,μ,l,-m,2)*p((x,y,z)=>(vx,vy,vz)))
                end
            end
            # 3: Summand 3 aus (I) - m > 0
            for λ=l+1:C.L
                for μ = 1:m-(l-λ)
                    p = rlm(λ-l,-(μ-m),x,y,z)
                    sum += abs(ω₊(C[λ,μ],λ,μ,l,m,3)*p((x,y,z)=>(vx,vy,vz)))
                end
            end
            # 4: Summand aus (II) - m = 0
            for λ=l-m:C.L
                p = rlm(λ-l,m,x,y,z)
                sum += abs((ω₀(C[λ,0],λ,l,-m,2)+ω₀(C[λ,0],λ,l,m,3))*p((x,y,z)=>(vx,vy,vz)))
            end
            # 5: Summand 1 aus (III) - m < 0
            for λ=l:C.L
                for μ = m-(λ-l):m#-1 #min(-1,m)
                    p = rlm(λ-l,abs(μ-m),x,y,z)
                    sum += abs(ω₋(C[λ,μ],λ,μ,l,m,1)*p((x,y,z)=>(vx,vy,vz)))
                end
            end
            # 6: Summand 2 aus (III) - m < 0
            for λ=l+1:C.L
                for μ = m+1:min(-1,m-(l-λ))
                    p = rlm(λ-l,μ-m,x,y,z)
                    sum += abs(ω₋(C[λ,μ],λ,μ,l,m,2)*p((x,y,z)=>(vx,vy,vz)))
                end
            end
            # 7: Summand 3 aus (III) - m < 0
            for λ=l+1:C.L
                for μ = -m-(λ-l):-1
                    p = rlm(λ-l,-(μ+m),x,y,z)
                    sum += abs(ω₋(C[λ,μ],λ,μ,l,-m,3)*p((x,y,z)=>(vx,vy,vz)))
                end
            end
            cShifted[l,m] = sqrt(4*pi/(2*l+1))*sum
            sum = 0
        end
    end

    # turn spherical coefficients into solid coefficients
    if cShifted.solid
        cShifted.solid = false
        solid!(cShifted)
        solid!(C)
    end

    return cShifted
end


text:
function DragForce(arg0::Atmosphere, arg1::DragSensitive)
    return DragForce((Atmosphere, DragSensitive), arg0, arg1)
end

function acceleration(obj::DragForce, arg0::FieldSpacecraftState, arg1::Vector{RealFieldElement})
    return jcall(obj, "acceleration", FieldVector3D, (FieldSpacecraftState, Vector{RealFieldElement}), arg0, arg1)
end

function acceleration(obj::DragForce, arg0::SpacecraftState, arg1::Vector{jdouble})
    return jcall(obj, "acceleration", Vector3D, (SpacecraftState, Vector{jdouble}), arg0, arg1)
end

function get_atmosphere(obj::DragForce)
    return jcall(obj, "getAtmosphere", Atmosphere, ())
end

function get_events_detectors(obj::DragForce)
    return jcall(obj, "getEventsDetectors", Stream, ())
end

function get_field_events_detectors(obj::DragForce, arg0::Field)
    return jcall(obj, "getFieldEventsDetectors", Stream, (Field,), arg0)
end

function get_parameters_drivers(obj::DragForce)
    return jcall(obj, "getParametersDrivers", Vector{ParameterDriver}, ())
end

function get_spacecraft(obj::DragForce)
    return jcall(obj, "getSpacecraft", DragSensitive, ())
end

ext: jl


text:
function getfreenodes(n, dirichletnodes)
	isfreenode = fill(true, n)
	isfreenode[dirichletnodes] .= false
	nodei2freenodei = fill(-1, length(isfreenode))
	freenodei2nodei = Array{Int}(undef, sum(isfreenode))
	j = 1
	for i = 1:length(isfreenode)
		if isfreenode[i]
			nodei2freenodei[i] = j
			freenodei2nodei[j] = i
			j += 1
		end
	end
	return isfreenode, nodei2freenodei, freenodei2nodei
end

function addboundaryconditions(hfree, dirichletnodes, dirichleths, isfreenode, nodei2freenodei)
	return map(i->isfreenode[i] ? hfree[nodei2freenodei[i]] : dirichleths[i], 1:length(dirichleths))
end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id$
#*--------------------------------------------------------------------
type Shell_Side_Main
	Shell_Side_Main()=begin
		new(
			length(Dict{Symbol,Any}(
				:Brief=>"Inside Shell Diameter",
				:Lower=>1E-6
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Shellside Fouling Resistance",
				:Unit=>"m^2*K/kW",
				:Symbol=>"Rf_{shell}",
				:Default=>1E-6 ,
				:Lower=>0
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Inlet Nozzle Inside Diameter",
				:Lower=>1E-6
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Outlet Nozzle Inside Diameter",
				:Lower=>1E-6
			)),
			Shell_Pdrop(Dict{Symbol,Any}(
				:Brief=>"Shell Side Pressure Drop",
				:Symbol=>" "
			)),
			Shell_Heat_Transfer(Dict{Symbol,Any}(
				:Brief=>"Shell Side Heat Transfer",
				:Symbol=>" "
			)),
			Physical_Properties(Dict{Symbol,Any}(
				:Brief=>"ShellSide Properties",
				:Symbol=>" "
			)),
			[:ShellID,:Fouling,:InletNozzleID,:OutletNozzleID,],
			[:PressureDrop,:HeatTransfer,:Properties,]
		)
	end
	ShellID::length
	Fouling::positive
	InletNozzleID::length
	OutletNozzleID::length
	PressureDrop::Shell_Pdrop
	HeatTransfer::Shell_Heat_Transfer
	Properties::Physical_Properties
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Shell_Side_Main
function atributes(in::Shell_Side_Main,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="Main variables in the Shell Side section of a shell and tube heat exchanger."
	drive!(fields,_)
	return fields
end
Shell_Side_Main(_::Dict{Symbol,Any})=begin
	newModel=Shell_Side_Main()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
@setup function()
    local len = 50
    tentacle = Segment(width/2, height, len, 0, 0)
    local current = tentacle
    local t = 0.0
    for i in 1:20
        t += 0.1
        len *= 0.75
        next = Segment(current, 20, 0, t)
        current.child = next
        current = next
    end

    xoff = 0
end

function update(dt)
    background("#111")

    next = tentacle
    while !isnothing(next)
        wiggle!(next)
        update!(next)
        draw(next)
        next = next.child
    end
end

max_stars_repo_name: jorge-brito/LynxGallery

max_issues_repo_name: jorge-brito/LynxGallery

max_forks_repo_name: jorge-brito/LynxGallery


text:
using CSV, Primes, DataFrames
a=[1 3 7 9];

# 6
P=Int64.(zeros(1,6));
X=[1 2 3 4 5 6];

for i=1:4
    for k=0:100
        for j=1:100
            for l=1:6
                P[l]=a[i]+10*k+10*j*(l-1);
            end
            if [isprime(P[1]) isprime(P[2]) isprime(P[3]) isprime(P[4]) isprime(P[5]) isprime(P[6]) ]==[Bool(1) Bool(1)  Bool(1) Bool(1) Bool(1) Bool(1)]
               X=[X;P];
            end
        end
    end
end

Y=DataFrame(X,:auto);
CSV.write("gt06.csv", Y);

# 7
P=Int64.(zeros(1,7));
X=[1 2 3 4 5 6 7];

for i=1:4
    for k=0:200
        for j=1:200
            for l=1:7
                P[l]=a[i]+10*k+10*j*(l-1);
            end
            if [isprime(P[1]) isprime(P[2]) isprime(P[3]) isprime(P[4]) isprime(P[5]) isprime(P[6]) isprime(P[7])]==[Bool(1) Bool(1)  Bool(1) Bool(1) Bool(1) Bool(1) Bool(1)]
               X=[X;P];
            end
        end
    end
end

Y=DataFrame(X,:auto);
CSV.write("gt07.csv", Y);

# 8
P=Int64.(zeros(1,8));
X=[1 2 3 4 5 6 7 8];

for i=1:4
    for k=0:500
        for j=1:500
            for l=1:8
                P[l]=a[i]+10*k+10*j*(l-1);
            end
            if [isprime(P[1]) isprime(P[2]) isprime(P[3]) isprime(P[4]) isprime(P[5]) isprime(P[6]) isprime(P[7]) isprime(P[8])]==[Bool(1) Bool(1)  Bool(1) Bool(1) Bool(1) Bool(1) Bool(1) Bool(1)]
               X=[X;P];
            end
        end
    end
end

Y=DataFrame(X,:auto);
CSV.write("gt08.csv", Y);

# 9
P=Int64.(zeros(1,9));
X=[1 2 3 4 5 6 7 8 9];

for i=1:4
    for k=0:1500
        for j=1:1500
            for l=1:9
                P[l]=a[i]+10*k+10*j*(l-1);
            end
            if [isprime(P[1]) isprime(P[2]) isprime(P[3]) isprime(P[4]) isprime(P[5]) isprime(P[6]) isprime(P[7]) isprime(P[8]) isprime(P[9])]==[Bool(1) Bool(1)  Bool(1) Bool(1) Bool(1) Bool(1) Bool(1) Bool(1) Bool(1)]
               X=[X;P];
            end
        end
    end
end

Y=DataFrame(X,:auto);
CSV.write("gt09.csv", Y);

# 10
P=Int64.(zeros(1,10));
X=[1 2 3 4 5 6 7 8 9 10];

for i=1:4
    for k=0:6000
        for j=1:6000
            for l=1:10
                P[l]=a[i]+10*k+10*j*(l-1);
            end
            if [isprime(P[1]) isprime(P[2]) isprime(P[3]) isprime(P[4]) isprime(P[5]) isprime(P[6]) isprime(P[7]) isprime(P[8]) isprime(P[9]) isprime(P[10])]==[Bool(1) Bool(1)  Bool(1) Bool(1) Bool(1) Bool(1) Bool(1) Bool(1) Bool(1) Bool(1)]
               X=[X;P];
            end
        end
    end
end

Y=DataFrame(X,:auto);
CSV.write("gt10.csv", Y);


text:
using  Plots
gr()
function euler(x₀, δ_t,nrsteps)
    x = x₀
    for i in 1:nrsteps
        dxdt = (x-3)*(x+6)*x
        x = x+ δ_t * dxdt
    end
    return x
end

# Each function should only have one task to do
function eulerarray(x0, deltat,nrsteps)
    x = Float64[] # Initialize an empty array
    append!(x, x0)
    for i in 1:nrsteps
        dxdt = (x[i] - 3) * (x[i] + 6) * x[i]
        temp = x[i] + deltat * dxdt
        append!(x,temp)
    end
    return x
end
"""
Demo function for Euler method
    - returns a plot
"""
function eulerarraydemo(x0,deltat,nrsteps)
    x = eulerarray(x0,deltat,nrsteps)
    plot(x, xlabel = "Time t", ylabel = "x",
            title = "Euler Example",
            legend=false)
end

function demo()
    println("x0? (try 2.5) :")
    x0 = chomp(readline())
    x0 = parse(Float64,x0)
    println("delta? (try 0.01) :")
    delta = chomp(readline())
    delta = parse(Float64,delta)
    println("n: (try 100)?")
    n = chomp(readline())
    n = parse(Int,n)
    eulerarraydemo(x0,delta,n)
end

println("Call demo() for an interactive session")
#demo()
#savefig(plot,"eulertest.png")

max_stars_repo_name: coinslab/ComputationalCognitiveModeling

max_issues_repo_name: coinslab/ComputationalCognitiveModeling

max_forks_repo_name: coinslab/ComputationalCognitiveModeling


text:
using AbstractPlotting.PlotUtils, AbstractPlotting.Colors

################################################################################
#                              Colormap reference                              #
################################################################################

function colors_svg(cs, w, h)
    n = length(cs)
    ws = min(w / n, h)
    html = """
    <?xml version="1.0" encoding="UTF-8"?>
    <!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN"
     "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
     <svg xmlns="http://www.w3.org/2000/svg" version="1.1"
          width="$(n * ws)mm" height="$(h)mm"
          viewBox="0 0 $n 1" preserveAspectRatio="none"
          shape-rendering="crispEdges" stroke="none">
    """
    for (i, c) in enumerate(cs)
        html *= """
        <rect width="$(ws)mm" height="$(h)mm" x="$(i-1)" y="0" fill="#$(hex(convert(RGB, c)))" />
        """
    end
    html *= "</svg>"
    return html
end

function generate_colorschemes_table(ks)
    extra_dir = get(ENV, "CI", "false") == "true" ? "../" : ""
    html = "<link type=\"text/css\" rel=\"stylesheet\" href=\"$(extra_dir)../assets/tables.css\" />

"
    for header in ["NAME", "Categorical variant", "Continuous variant"]
        html *= "$header	"
    end
    html *= "
"
    w, h = 70, 5
    for k in ks
        grad = cgrad(k)
        p = color_list(grad)
        cg = grad[range(0, 1, length = 100)]
        cp = length(p) <= 100 ? p : cg
        # cp7 = color_list(palette(k, 7))

        html *= "
:$k	"
        html *= colors_svg(cp, w, h)
        html *= "	"
        html *= colors_svg(cg, w, h)
        # html *= "	"
        # html *= colors_svg(cp7, 35, h)
        html *= "
"
    end
    html *= "
"
    return html
end

function generate_colorschemes_markdown(; GENDIR = joinpath(dirname(@__DIR__), "docs", "src", "generated"))
    md = open(joinpath(GENDIR, "colors.md"), "w")

    for line in readlines(normpath(@__DIR__, "..", "docs", "src", "colors.md"))
        write(md, line)
        write(md, "\n")
    end

    write(md, """
    ## misc
    These colorschemes are not defined or provide different colors in ColorSchemes.jl
    They are kept for compatibility with the old behaviour of Makie, before v0.10.
    """)
    write(md, "```@raw html\n")
    write(
        md,
        generate_colorschemes_table(
            [:default; sort(collect(keys(PlotUtils.MISC_COLORSCHEMES)))]
        )
    )
    write(md, "\n```\n\nThe following colorschemes are defined by ColorSchemes.jl.\n\n")
    for cs in ["cmocean", "scientific", "matplotlib", "colorbrewer", "gnuplot", "colorcet", "seaborn", "general"]
        ks = sort([k for (k, v) in PlotUtils.ColorSchemes.colorschemes if occursin(cs, v.category)])
        write(md, "\n\n## $cs\n\n```@raw html\n")
        write(md, generate_colorschemes_table(ks))
        write(md, "\n```\n\n")
    end

    close(md)
end


text:
using GeometryTypes: Point, Vec

import ConvexBodyProximityQueries.support
using GeometryTypes: HyperSphere, HyperRectangle, HyperCube

function ConvexBodyProximityQueries.support(sphere::HyperSphere{N, T}, dir::AbstractVector) where {N, T}
    SVector{N}(sphere.center + sphere.r*normalize(dir, 2))
end

@generated function ConvexBodyProximityQueries.support(rect::HyperRectangle{N, T}, dir::AbstractVector) where {N, T}
    exprs = Array{Expr}(undef, (N,))
    for i = 1:N
        exprs[i] = :(rect.widths[$i]*(dir[$i] ≥ 0.0 ? 1.0 : -1.0)/2.0 + rect.origin[$i])
    end

    return quote
        Base.@_inline_meta
        @inbounds elements = tuple($(exprs...))
        @inbounds return SVector{N, T}(elements)
    end
end

@generated function ConvexBodyProximityQueries.support(cube::HyperCube{N, T}, dir::AbstractVector) where {N, T}
    exprs = Array{Expr}(undef, (N,))
    for i = 1:N
        exprs[i] = :(cube.width*(dir[$i] ≥ 0.0 ? 1.0 : -1.0)/2.0 + cube.origin[$i])
    end

    return quote
        Base.@_inline_meta
        @inbounds elements = tuple($(exprs...))
        @inbounds return SVector{N, T}(elements)
    end
end

@testset "gjk (2D)" begin
    vec2 = SVector{2}([1.0, 0.0])
    @testset "array of vertices" begin
        rect2 = SMatrix{3,4}([0.0 1.0 1.0 0.0; 0.0 0.0 2.0 2.0; 1.0 1.0 1.0 1.0])
        transform2(θ,x,y) = SMatrix{2,3}([cos(θ) -sin(θ) x; sin(θ) cos(θ) y])

        polyA = transform2(0,0,0) * rect2
        polyB = transform2(0,3,0) * rect2
        ret = closest_points(polyA, polyB, vec2)
        @test collision_detection(polyA, polyB, vec2) == false
        @test all(ret[1] .≈ [1.0, 0.0])
        @test all(ret[2] .≈ [3.0, 0.0])

        polyA = transform2(0,0,0) * rect2
        polyB = transform2(0,1,0) * rect2
        ret = closest_points(polyA, polyB, vec2)
        @test collision_detection(polyA, polyB, vec2) == true

        polyA = transform2(0,0,0) * rect2
        polyB = transform2(0,0.5,0) * rect2
        ret = closest_points(polyA, polyB, vec2)
        @test collision_detection(polyA, polyB, vec2) == true

        polyA = transform2(0,0,0) * rect2
        polyB = transform2(π/4,3,0) * rect2
        ret = closest_points(polyA, polyB, vec2)
        @test collision_detection(polyA, polyB, vec2) == false
        @test all(ret[1] .≈ [1.0, 1.4142135623730951])
        @test all(ret[2] .≈ [1.585786437626905, 1.4142135623730954])

        polyA = transform2(0,0,0) * rect2
        polyB = transform2(π/4,1,0) * rect2
        ret = closest_points(polyA, polyB, vec2)
        @test collision_detection(polyA, polyB, vec2) == true

        polyA = transform2(0,0,0) * rect2
        polyB = transform2(π/4,0.5,0) * rect2
        ret = closest_points(polyA, polyB, vec2)
        @test collision_detection(polyA, polyB, vec2) == true

        polyA = transform2(-π/4,0,0) * rect2
        polyB = transform2(π/4,3,0) * rect2
        ret = closest_points(polyA, polyB, vec2)
        @test collision_detection(polyA, polyB, vec2) == false
        @test all(ret[1] .≈ [1.914213562373095, 0.9142135623730953])
        @test all(ret[2] .≈ [2.0, 1.0])

        polyA = transform2(-π/4,0,0) * rect2
        polyB = transform2(π/4,1,0) * rect2
        ret = closest_points(polyA, polyB, vec2)
        @test collision_detection(polyA, polyB, vec2) == true

        polyA = transform2(-π/4,0,0) * rect2
        polyB = transform2(π/4,0.5,0) * rect2
        ret = closest_points(polyA, polyB, vec2)
        @test collision_detection(polyA, polyB, vec2) == true
    end # array of vertices

    @testset "sphere" begin
        sphereA = HyperSphere(Point(0.0, 0.0), 1.0)
        sphereB = HyperSphere(Point(3.0, 0.0), 1.0)
        ret = closest_points(sphereA, sphereB, vec2)
        @test collision_detection(sphereA, sphereB, vec2) == false
        @test all(ret[1] .≈ [1.0, 0.0])
        @test all(ret[2] .≈ [2.0, 0.0])

        sphereA = HyperSphere(Point(0.0, 0.0), 1.0)
        sphereB = HyperSphere(Point(2.0, 0.0), 1.0)
        ret = closest_points(sphereA, sphereB, vec2)
        @test collision_detection(sphereA, sphereB, vec2) == true

        sphereA = HyperSphere(Point(0.0, 0.0), 1.0)
        sphereB = HyperSphere(Point(1.0, 0.0), 1.0)
        ret = closest_points(sphereA, sphereB, vec2)
        @test collision_detection(sphereA, sphereB, vec2) == true
    end # sphere

    @testset "rectangle" begin
        rectangleA = HyperRectangle(Vec(0.0, 0.0), Vec(1.0, 2.0))
        rectangleB = HyperRectangle(Vec(3.0, 0.0), Vec(1.0, 2.0))
        ret = closest_points(rectangleA, rectangleB, vec2)
        @test collision_detection(rectangleA, rectangleB, vec2) == false
        @test all(ret[1] .≈ [0.5, 1.0])
        @test all(ret[2] .≈ [2.5, 1.0])

        rectangleA = HyperRectangle(Vec(0.0, 0.0), Vec(1.0, 2.0))
        rectangleB = HyperRectangle(Vec(1.0, 0.0), Vec(1.0, 2.0))
        ret = closest_points(rectangleA, rectangleB, vec2)
        @test collision_detection(rectangleA, rectangleB, vec2) == true

        rectangleA = HyperRectangle(Vec(0.0, 0.0), Vec(1.0, 2.0))
        rectangleB = HyperRectangle(Vec(0.5, 0.0), Vec(1.0, 2.0))
        ret = closest_points(rectangleA, rectangleB, vec2)
        @test collision_detection(rectangleA, rectangleB, vec2) == true
    end # rectangle

    @testset "cube" begin
        cubeA = HyperCube(Vec(0.0, 0.0), 1.0)
        cubeB = HyperCube(Vec(3.0, 0.0), 1.0)
        ret = closest_points(cubeA, cubeB, vec2)
        @test collision_detection(cubeA, cubeB, vec2) == false
        @test all(ret[1] .≈ [0.5, 0.5])
        @test all(ret[2] .≈ [2.5, 0.5])

        cubeA = HyperCube(Vec(0.0, 0.0), 1.0)
        cubeB = HyperCube(Vec(1.0, 0.0), 1.0)
        ret = closest_points(cubeA, cubeB, vec2)
        @test collision_detection(cubeA, cubeB, vec2) == true

        cubeA = HyperCube(Vec(0.0, 0.0), 1.0)
        cubeB = HyperCube(Vec(0.5, 0.0), 1.0)
        ret = closest_points(cubeA, cubeB, vec2)
        @test collision_detection(cubeA, cubeB, vec2) == true
    end # cube
end # gjk

@testset "gjk (3D)" begin
    vec3 = SVector{3}([1.0, 0.0, 0.0])
    @testset "array of vertices" begin
        # TODO (add several tetrahedron configurations)
    end # array of vertices

    @testset "sphere" begin
        sphereA = HyperSphere(Point(0.0, 0.0, 0.0), 1.0)
        sphereB = HyperSphere(Point(3.0, 3.0, 0.0), 1.0)
        ret = closest_points(sphereA, sphereB, vec3)
        @test collision_detection(sphereA, sphereB, vec3) == false
        @test all(ret[1] .≈ [0.7071067841796789, 0.7071067781934159, 0.0])
        @test all(ret[2] .≈ [2.2928932158203210, 2.2928932218065840, 0.0])

        sphereA = HyperSphere(Point(0.0, 0.0, 0.0), 1.0)
        sphereB = HyperSphere(Point(2.0, 0.0, 0.0), 1.0)
        ret = closest_points(sphereA, sphereB, vec3)
        @test collision_detection(sphereA, sphereB, vec3) == true

        sphereA = HyperSphere(Point(0.0, 0.0, 0.0), 1.0)
        sphereB = HyperSphere(Point(1.0, 0.0, 0.0), 1.0)
        ret = closest_points(sphereA, sphereB, vec3)
        @test collision_detection(sphereA, sphereB, vec3) == true
    end # sphere

    @testset "rectangle" begin
        rectangleA = HyperRectangle(Vec(0.0, 0.0, 0.0), Vec(1.0, 2.0, 3.0))
        rectangleB = HyperRectangle(Vec(3.0, 3.0, 0.0), Vec(1.0, 2.0, 3.0))
        ret = closest_points(rectangleA, rectangleB, vec3)
        @test collision_detection(rectangleA, rectangleB, vec3) == false
        @test all(ret[1] .≈ [0.5, 1.0, 1.5])
        @test all(ret[2] .≈ [2.5, 2.0, 1.5])

        rectangleA = HyperRectangle(Vec(0.0, 0.0, 0.0), Vec(1.0, 2.0, 3.0))
        rectangleB = HyperRectangle(Vec(1.0, 2.0, 0.0), Vec(1.0, 2.0, 3.0))
        ret = closest_points(rectangleA, rectangleB, vec3)
        @test collision_detection(rectangleA, rectangleB, vec3) == true

        rectangleA = HyperRectangle(Vec(0.0, 0.0, 0.0), Vec(1.0, 2.0, 3.0))
        rectangleB = HyperRectangle(Vec(0.5, 0.0, 0.0), Vec(1.0, 2.0, 3.0))
        ret = closest_points(rectangleA, rectangleB, vec3)
        @test collision_detection(rectangleA, rectangleB, vec3) == true
    end # rectangle

    @testset "cube" begin
        cubeA = HyperCube(Vec(0.0, 0.0, 0.0), 1.0)
        cubeB = HyperCube(Vec(3.0, 3.0, 0.0), 1.0)
        ret = closest_points(cubeA, cubeB, vec3)
        @test collision_detection(cubeA, cubeB, vec3) == false
        @test all(ret[1] .≈ [0.5, 0.5, 0.5])
        @test all(ret[2] .≈ [2.5, 2.5, 0.5])

        cubeA = HyperCube(Vec(0.0, 0.0, 0.0), 1.0)
        cubeB = HyperCube(Vec(1.0, 1.0, 0.0), 1.0)
        ret = closest_points(cubeA, cubeB, vec3)
        @test collision_detection(cubeA, cubeB, vec3) == true

        cubeA = HyperCube(Vec(0.0, 0.0, 0.0), 1.0)
        cubeB = HyperCube(Vec(0.5, 0.0, 0.0), 1.0)
        ret = closest_points(cubeA, cubeB, vec3)
        @test collision_detection(cubeA, cubeB, vec3) == true
    end # cube
end # gjk


text:
## TEAVAR TE
function TEAVAR(GRB_ENV, edges, capacity, flows, demand, beta, T, Tf, scenarios, scenario_probs; average=false)
    printstyled("\n** solving TEAVAR LP at beta $(beta)..\n", color=:yellow)
    nedges = length(edges)
    nflows = length(flows)
    ntunnels = length(T)
    nscenarios = length(scenarios)
    p = scenario_probs

    tunnel_num = 0
    for x in 1:length(Tf)
        if size(Tf[x],1) > tunnel_num
            tunnel_num = size(Tf[x],1)
        end
    end

    #CREATE TUNNEL SCENARIO MATRIX
    X  = ones(nscenarios,ntunnels)
    for s in 1:nscenarios
        for t in 1:ntunnels
            if size(T[t],1) == 0
                X[s,t] = 0
            else
                for e in 1:nedges
                    if scenarios[s][e] == 0
                        back_edge = findfirst(x -> x == (edges[e][2],edges[e][1], edges[e][3]), edges)
                        if in(e, T[t]) || in(back_edge, T[t])
                        # if in(e, T[t])
                            X[s,t] = 0
                        end
                    end
                end
            end
        end
    end

    #CREATE TUNNEL EDGE MATRIX
    L = zeros(ntunnels, nedges)
    for t in 1:ntunnels
        for e in 1:nedges
            if in(e, T[t])
                L[t,e] = 1
            end
        end
    end

    model = Model(() -> Gurobi.Optimizer(GRB_ENV))
    set_optimizer_attribute(model, "OutputFlag", 0)
    set_optimizer_attribute(model, "Threads", 32)

    model_def_time = @elapsed begin

    @variable(model, a[1:nflows, 1:tunnel_num] >= 0)
    @variable(model, alpha >= 0)
    @variable(model, umax[1:nscenarios] >= 0)
    @variable(model, u[1:nscenarios, 1:nflows] >= 0)

    for e in 1:nedges
        @constraint(model, sum(a[f,t] * L[Tf[f][t],e] for f in 1:nflows for t in 1:size(Tf[f],1)) <= capacity[e])   #overlapping flows cannot add up to the capacity of that link
    end

    # FLOW LEVEL LOSS
    @expression(model, satisfied[s=1:nscenarios, f=1:nflows], sum(a[f,t] * X[s,Tf[f][t]] for t in 1:size(Tf[f],1)) / demand[f])

    for s in 1:nscenarios
        for f in 1:nflows
            @constraint(model, u[s,f] == 1 - satisfied[s,f])
        end
    end

    for s in 1:nscenarios
        if average
            @constraint(model, umax[s] + alpha >= (sum(u[s,f] for f in 1:nflows)) / nflows)
        else
            for f in 1:nflows
                @constraint(model, umax[s] >= u[s,f] - alpha)
            end
        end
    end

    @objective(model, Min, alpha + (1 / (1 - beta)) * sum((p[s] * umax[s] for s in 1:nscenarios)))
    end  # model_def_time end

    optimize!(model)
    solve_runtime = solve_time(model)
    opt_runtime = solve_runtime + model_def_time

    return value.(alpha), objective_value(model), value.(a), value.(umax), value.(u), solve_runtime, opt_runtime
end

max_stars_repo_name: hipersys-team/arrow

max_issues_repo_name: arrow-bow/arrow

max_forks_repo_name: arrow-bow/arrow


text:
# dual only versus dual only no scaling

using JLD

include("perfplots.jl")

timeshift = 10.0

# produced by process_csv.jl in perfprofile mode
jldsource = joinpath(dirname(@__FILE__),"subpscaling_perf.jld")

d = load(jldsource)
time_table = d["time_table"] + timeshift
itercount_table = d["itercount_table"]
solvers = d["solvers"]
@assert solvers == ["PAJ_CPLEX_MOSEK_dualonly","PAJ_CPLEX_MOSEK_dualonlynoscale"]
instance_names = d["instance_names"]

Plots.scalefontsizes(2.3)

performance_profile(time_table, ["scaling","no scaling"],logscale=false, ymax=1.0, xmax=4, linewidth=3, xticks=[1,2,3,4])
Plots.savefig(joinpath(dirname(@__FILE__),"scaletime.eps"))

performance_profile(itercount_table, ["scaling","no scaling"],logscale=false, ymax=1.0, xmax=4, linewidth=3, xticks=[1,2,3,4])
Plots.savefig(joinpath(dirname(@__FILE__),"scaleiter.eps"))

max_stars_repo_name: mlubin/PajaritoSupplement

max_issues_repo_name: mlubin/PajaritoSupplement

max_forks_repo_name: mlubin/PajaritoSupplement


text:
using Vulkan
using Base.Test

const api = vk.api

err = api.VkResult(0)

toversion(version::Cuint) = VersionNumber(version >> 22,  (version >> 12) & 0x3ff, version & 0xfff)

count = Ref{Cuint}(0)
# Scan layers
err = api.vkEnumerateInstanceLayerProperties(count, C_NULL)
assert(err == api.VK_SUCCESS)
global_layer_properties = Array(api.VkLayerProperties, count[])
err = api.vkEnumerateInstanceLayerProperties(count, global_layer_properties)
assert(err == api.VK_SUCCESS)

function Base.show(io::IO, lp::api.VkLayerProperties)
	println(io, "Layer Properties: ")
	println(io, "    Layer Name: ", bytestring(filter(x->x!=0, UInt8[lp.layerName...])))
	println(io, "    Spec Version: ", toversion(lp.specVersion))
	println(io, "    Implementation Version: ", toversion(lp.implementationVersion))
	println(io, "    description: ", bytestring(filter(x->x!=0, UInt8[lp.description...])))
end
for elem in global_layer_properties
	println(elem)
end

appname = b"vulkaninfo"

app_info = Ref(api.VkApplicationInfo(
    api.VK_STRUCTURE_TYPE_APPLICATION_INFO,
    C_NULL,
    pointer(appname),
    1,
    pointer(appname),
    1,
    api.VK_VERSION,
))

inst_info = Ref(api.VkInstanceCreateInfo(
        api.VK_STRUCTURE_TYPE_INSTANCE_CREATE_INFO,
        C_NULL,
        UInt32(0),
        Base.unsafe_convert(Ptr{api.VkApplicationInfo}, app_info),
        0,
        C_NULL,
        0,
        C_NULL,
))

instance = Ref{api.VkInstance}(C_NULL)

err = api.vkCreateInstance(inst_info, C_NULL, instance)
println(err)
println(instance)

gpu_count = Ref{Cuint}(0)
err = api.vkEnumeratePhysicalDevices(instance[], gpu_count, C_NULL)
println(err)
devices = Array(api.VkPhysicalDevice, gpu_count[])

err = api.vkEnumeratePhysicalDevices(instance[], gpu_count, devices)

println(err)

deviceprops = Ref{api.VkPhysicalDeviceProperties}()
err = api.vkGetPhysicalDeviceProperties(devices[], deviceprops)

println(err)

function Base.show(io::IO, pdsp::api.VkPhysicalDeviceSparseProperties)
	println(io, "Physical device Sparse Properties: ")
	for name in fieldnames(pdsp)
		println(io, "    ", name, " ", pdsp.(name) == 1)
	end
end
function Base.show(io::IO, pdp::api.VkPhysicalDeviceProperties)
	println(io, "Physical Device Properties: ")
	println(io, "    API Version: ", toversion(pdp.apiVersion))
	println(io, "    Driver Version: ", toversion(pdp.driverVersion))

	println(io, "    Vendor ID ", pdp.vendorID)
	println(io, "    Device ID: ", pdp.deviceID)
	println(io, "    Device Type: ", pdp.deviceType)
	println(io, "    Device Name: ", bytestring(filter(x->x!=0, UInt8[pdp.deviceName...])))
	println(io, "    Pipeline Cache UUID: ", pdp.pipelineCacheUUID)
	println(io, "    Limits: ", pdp.limits)
	println(io, "    Sparse Properties: \n    ", pdp.sparseProperties)
end

limitshow(x::Cuint) = Int(x)
limitshow(x::NTuple) = "<"*join(map(limitshow, x), " ")*">"
limitshow(x) = x

function Base.show(io::IO, pdl::api.VkPhysicalDeviceLimits)
	println(io, "Physical Device Limits: ")
	for name in fieldnames(pdl)
		println(io, "    ", name, " ", limitshow(pdl.(name)))
	end
end

println(deviceprops[])
queue_count = Ref{Cuint}(0)

api.vkGetPhysicalDeviceQueueFamilyProperties(devices[], queue_count, C_NULL)
queueprops = Array(api.VkQueueFamilyProperties, queue_count[])
println(queue_count[])
api.vkGetPhysicalDeviceQueueFamilyProperties(devices[], queue_count, queueprops)
println(queueprops)

memprops = Ref{api.VkPhysicalDeviceMemoryProperties}()
api.vkGetPhysicalDeviceMemoryProperties(devices[], memprops)

println(memprops[])
devicefeatures = Ref{api.VkPhysicalDeviceFeatures}()
api.vkGetPhysicalDeviceFeatures(devices[], devicefeatures)

function Base.show(io::IO, df::api.VkPhysicalDeviceFeatures)
	println(io, "Physical device features: ")
	for name in fieldnames(df)
		println(io, "    ", name, ": ", df.(name) != 1 ? "un" : "", "supported")
	end
end
println(devicefeatures[])

# app_dev_init(&gpu->dev, gpu);
# app_dev_init_formats(&gpu->dev);

include("windowing.jl")

connection, scr = XCB.connect(C_NULL)
println(connection)
setup = XCB.get_setup(connection)
println(setup)
iter = Ref(XCB.setup_roots_iterator(setup))
println(iter)

for i=scr:-1:0
    XCB.screen_next(iter)
end
screen = unsafe_load(iter[].data, 1)

window = XCB.generate_id(connection)
value_mask = UInt32(XCB.CW_BACK_PIXEL) | UInt32(XCB.CW_EVENT_MASK)
value_list = zeros(UInt32, 32)
value_list[1] = screen.black_pixel
value_list[2] = (
	UInt32(XCB.EVENT_MASK_KEY_RELEASE) |
	UInt32(XCB.EVENT_MASK_EXPOSURE) |
    UInt32(XCB.EVENT_MASK_STRUCTURE_NOTIFY)
)
println(screen)

surface = Ref{api.VkSurfaceKHR}(C_NULL)
XCB.create_window(
	connection, XCB.COPY_FROM_PARENT, window,
	screen.root, 0, 0, 512, 512, 0,
	XCB.WINDOW_CLASS_INPUT_OUTPUT, screen.root_visual,
	value_mask, value_list
)

@windows_only begin
	createInfo = Ref(api.VkWin32SurfaceCreateInfoKHR(
		VK_STRUCTURE_TYPE_WIN32_SURFACE_CREATE_INFO_KHR,
    	C_NULL,
    	0,
    	connection,
    	window
	))

    err = vkCreateWin32SurfaceKHR(instance[], createInfo, C_NULL, surface);
end

@unix_only begin
	createInfo = Ref(api.VkXcbSurfaceCreateInfoKHR(
		api.VK_STRUCTURE_TYPE_XCB_SURFACE_CREATE_INFO_KHR,
		C_NULL,
		0,
		connection,
		window
	))

	err = api.vkCreateXcbSurfaceKHR(instance[], createInfo, C_NULL, surface)
	println(err)
end


text:
using Test
using Statistics: mean
using Transducers
using Transducers: finaltype, simple_transduce, Reduced

inc(x) = x + oneunit(x)

isfiniteiter(iter) =
    Base.IteratorSize(iter) isa Union{Base.HasLength, Base.HasShape}

function iterator_variants(xs)
    iters = Any[xs]
    if !(xs isa Array) && isfiniteiter(xs)
        push!(iters, collect(xs))
    end
    push!(iters, Base.Generator(identity, xs))
    return iters
end


text:
using RunTests
using Base.Test

@testmodule ErrorAnalysisTest begin
  using MLLibrary

  function test_confusionMatrix()
    testset = [1 0 1 1 1 1 1 0 0 0 0]
    predset = [1 1 0 1 1 1 0 1 0 0 0]

    @test MLLibrary.Utils.ErrorAnalysis.generateConfusionMatrix(testset,predset)==[0 0 1;
                    0 3 2;
                    1 2 4;]

    testset = [1 2 3 4 1 2 4 1 3 2 1 2]
    predset = [1 1 4 1 1 2 3 2 3 2 1 4]

    @test MLLibrary.Utils.ErrorAnalysis.generateConfusionMatrix(testset,predset)==[0 1 2 3 4;
                    1 3 1 0 0;
                    2 1 2 0 1;
                    3 0 0 1 1;
                    4 1 0 1 0;]
  end

  function test_generateWeightedErrorRate()
    testset = [1 0 1 1 1 1 1 0 0 0 0]
    predset = [1 1 0 1 1 1 0 1 0 0 0]

    confusionmatrix =  MLLibrary.Utils.ErrorAnalysis.generateConfusionMatrix(testset,predset)

    @test round(MLLibrary.Utils.ErrorAnalysis.calculateWeightedSuccessRate(confusionmatrix),2) == .63

    testset = [1 2 3 4 1 2 4 1 3 2 1 2]
    predset = [1 1 4 1 1 2 3 2 3 2 1 4]

    confusionmatrix= MLLibrary.Utils.ErrorAnalysis.generateConfusionMatrix(testset,predset)

    @test round(MLLibrary.Utils.ErrorAnalysis.calculateWeightedSuccessRate(confusionmatrix),2) == .44

  end

  function test_generateErrorRate()
    testset = [1 0 1 1 1 1 1 0 0 0 0]
    predset = [1 1 0 1 1 1 0 1 0 0 0]

    @test round(MLLibrary.Utils.ErrorAnalysis.calculateErrorRate(testset,predset),2)==.36
    #
    testset = [1 2 3 4 1 2 4 1 3 2 1 2]
    predset = [1 1 4 1 1 2 3 2 3 2 1 4]

    @test MLLibrary.Utils.ErrorAnalysis.calculateErrorRate(testset,predset)==.5

  end

  function test_calculateRootMeanStandardError()
    y=[1;2;3;4;5;6;1;2;3]
    observed = [0;4;2;1;5;2;1;6;1]

    @test round(MLLibrary.Utils.ErrorAnalysis.calculateRootMeanStandardError(observed,y),2)==2.38

  end

  function test_calculateWCSS()
    centroids=[1 2 3 4;5 6 7 8]
    data=[1 2 1 1;4 1 2 2; 5 2 1 4;3 3 2 1]
    y=[1;1;2;2]
    wcss = MLLibrary.Utils.ErrorAnalysis.calculateWCSS(centroids,data,y)
    println(wcss)

    @test wcss==183
  end

end

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
using ContactImplicitMPC

using Test
using ForwardDiff
using JLD2
using QDLDL
using Symbolics
using StaticArrays
using LinearAlgebra
using Random
using SparseArrays
using BenchmarkTools
using Quaternions

# Solver
include("solver/qr.jl")
include("solver/ldl.jl")
include("solver/lu.jl")
include("solver/schur.jl")
include("solver/ip_solver.jl")

# Dynamics
include("dynamics/lagrangian.jl") # need to fix d_fast
include("dynamics/model.jl")
include("dynamics/particle.jl")
include("dynamics/quadruped.jl")
include("dynamics/quaternion.jl")

# Simulator
include("simulator/rotations.jl")
include("simulator/environment.jl")
include("simulator/trajectory.jl")
include("simulator/simulator.jl")
include("simulator/particle.jl")
# include("simulator/hopper_2D.jl") #TODO: set tests to raibert model
include("simulator/hopper_3D.jl")
include("simulator/quadruped.jl")
include("simulator/open_loop.jl")
include("simulator/flamingo.jl")

# Controller
include("controller/objective.jl")
include("controller/linearized_step.jl")
include("controller/implicit_dynamics.jl")
include("controller/linearized_solver.jl")
include("controller/newton.jl")
include("controller/newton_structure_solver.jl") # fails on github actions

# MPC examples
include("controller/mpc_quadruped.jl")
include("controller/mpc_flamingo.jl")


text:
abstract type AbstractGraphFormat end

"""
    loadgraph(file, gname="graph", format=LGFormat())

Read a graph named `gname` from `file` in the format `format`.

### Implementation Notes
`gname` is graph-format dependent and is only used if the file contains
multiple graphs; if the file format does not support multiple graphs, this
value is ignored. The default value may change in the future.
"""
function loadgraph(fn::AbstractString, gname::AbstractString, format::AbstractGraphFormat)
    open(fn, "r") do io
        loadgraph(auto_decompress(io), gname, format)
    end
end
loadgraph(fn::AbstractString) = loadgraph(fn, "graph", LGFormat())
loadgraph(fn::AbstractString, gname::AbstractString) = loadgraph(fn, gname, LGFormat())
loadgraph(fn::AbstractString, format::AbstractGraphFormat) = loadgraph(fn, "graph", format)

"""
    loadgraphs(file, format=LGFormat())

Load multiple graphs from `file` in the format `format`.
Return a dictionary mapping graph name to graph.

### Implementation Notes
For unnamed graphs the default name \"graph\" will be used. This default
may change in the future.
"""
function loadgraphs(fn::AbstractString, format::AbstractGraphFormat)
    open(fn, "r") do io
        loadgraphs(auto_decompress(io), format)
    end
end

loadgraphs(fn::AbstractString) = loadgraphs(fn, LGFormat())

function auto_decompress(io::IO)
    format = :raw
    mark(io)
    if !eof(io)
        b1 = read(io, UInt8)
        if !eof(io)
            b2 = read(io, UInt8)
            if (b1, b2) == (0x1f, 0x8b)  # check magic bytes
                format = :gzip
            end
        end
    end
    reset(io)
    if format == :gzip
        io = InflateGzipStream(io)
    end
    return io
end

"""
    savegraph(file, g, gname="graph", format=LGFormat)

Saves a graph `g` with name `gname` to `file` in the format `format`.
Return the number of graphs written.

### Implementation Notes
The default graph name assigned to `gname` may change in the future.
"""
function savegraph(fn::AbstractString, g::AbstractGraph, gname::AbstractString,
        format::AbstractGraphFormat; compress=nothing
    )
    if compress !== nothing
        if !compress
            @info("Note: the `compress` keyword is no longer supported in LightGraphs. Saving uncompressed.")
        else
            Base.depwarn("Saving compressed graphs is no longer supported in LightGraphs. Use `LGCompressedFormat()` from the `GraphIO.jl` package instead. Saving uncompressed.", :savegraph)
        end
    end
    io = open(fn, "w")
    try
        return savegraph(io, g, gname, format)
    catch
        rethrow()
    finally
        close(io)
    end
end

# without graph name
savegraph(fn::AbstractString, g::AbstractGraph, format::AbstractGraphFormat; compress=nothing) =
    savegraph(fn, g, "graph", format, compress=compress)

# without format - default to LGFormat()
savegraph(fn::AbstractString, g::AbstractSimpleGraph; compress=nothing) = savegraph(fn, g, "graph", LGFormat(), compress=compress)

"""
    savegraph(file, g, d, format=LGFormat)

Save a dictionary of `graphname => graph` to `file` in the format `format`.
Return the number of graphs written.

### Implementation Notes
Will only work if the file format supports multiple graph types.
"""
function savegraph(fn::AbstractString, d::Dict{T,U},
    format::AbstractGraphFormat; compress=nothing) where T <: AbstractString where U <: AbstractGraph
    compress === nothing ||
    Base.depwarn("Saving compressed graphs is no longer supported in LightGraphs. Use `LGCompressedFormat()` from the `GraphIO.jl` package instead. Saving uncompressed.", :savegraph)
        io = open(fn, "w")

    try
        return savegraph(io, d, format)
    catch
        rethrow()
    finally
        close(io)
    end
end

savegraph(fn::AbstractString, d::Dict) = savegraph(fn, d, LGFormat())


text:
using FFTW,LinearAlgebra,Plots,SpecialFunctions;pic=()->(N=64;P=32N;dt=1/6N;
T=2^13;W=32π^2/3;w=W/P*N;ξ=im*zeros(N);x=(bitreverse.(0:P-1).+2.0^63)/2.0^64;
v=collect(1:P.>P/2).*2 .-1.;r=zeros(N);X=copy(x);V=copy(v);l=1e-14;
E=zeros(3,N);F=copy(E);ik=2π*im*vcat(1,1:N/2,-N/2+1:-1);D=zeros(T,4);
d(y)=(i=Int(mod1(ceil(y*N),N));o=ceil(y*N)-y*N;((i,1-o),(mod1(i-1,N),o)));
ρ(x,y=x)=(r.*=0;for l∈d.((x.+y)./2);for k∈l;r[k[1]]+=k[2]*w;end;end;r);
@time @gif for t ∈1:T;X.=x;V.=v;F.*=NaN;
 E[1,:].=real.(ifft((ξ.=fft(ρ(X,X))./ik;ξ[1]*=0;ξ)));
 for _∈0:9;≈(F,E,rtol=l,atol=0)&&break;F.=E;
 x.=X.+(v.+V)/2*dt;
 @inbounds for j∈1:P;v[j]=V[j]+sum(k->E[1,k[1]]*k[2],d((X[j]+X[j])/2))*dt/6;end;
 E[2,:].=real.(ifft((ξ.=fft(ρ(X,x))./ik;ξ[1]*=0;ξ)));
 @inbounds for j∈1:P;v[j]=v[j]+sum(k->E[2,k[1]]*k[2],d((X[j]+x[j])/2))*4dt/6;end;
 E[3,:].=real.(ifft((ξ.=fft(ρ(x,x))./ik;ξ[1]*=0;ξ)));
 @inbounds for j∈1:P;v[j]=v[j]+sum(k->E[3,k[1]]*k[2],d((x[j]+x[j])/2))*dt/6;end;
 end;x.=mod.(x,1);
 D[t,1:2].=(sum(E[end,:].^2)/N,sum(v.^2)*W/P)./2;D[t,3:4].=sum.((D[t,1:2],v/P));
 scatter(x[1:P÷2],v[1:P÷2],c=1);scatter!(x[P÷2+1:P],v[P÷2+1:P],ylims=(-3,3),c=2);
 D[t,1:3].*=2/W;plot!((1:t)./T,D[1:t,:],legend=0==1,xlims=(0,1)); @show t/T
 plot!((0.5:N)/N,ρ(x[1:P÷2])./W,c=1);plot!((0.5:N)/N,ρ(x[P÷2+1:P])./W,c=2)
end every 8;return (D,N,P,dt,T,W,w););(D,N,P,dt,T,W,w)=pic()
la(x)=log10(abs(x));t=(1:T).*dt;plot(t,la.(D[:,1]),label="electric field energy")
plot!(t,la.(1 .-D[:,3]),label="total energy change");yticks!(-35:1)
plot!(t,la.(D[:,4]),label="total momentum",legend=:bottomright)
γ(x)=imag(sqrt(Complex(x^2+1-sqrt(4*x^2+1))))*sqrt(W/2)/log(10);ylims!(-35,1)
plot!(t,2*γ(2π/sqrt(W/2)).*(t.-T÷8*dt).+la.(D[T÷8,1]),label="predicted")
xlabel!("Time");ylabel!("log 10");savefig("AreaFixedPointQuietSimpson13.pdf")


text:
abstract type NeutralElement end
abstract type HasNeutralElement <: NeutralElement end
struct HasZero <: HasNeutralElement end
struct HasOne <: HasNeutralElement end
struct HasZeroAndOne <: HasNeutralElement end
struct HasNoNeutralElement <: NeutralElement end

abstract type Idempotency end
struct IsIdempotent <: Idempotency end

NeutralElement(::Type{T}) where T<:Number = HasZeroAndOne()
NeutralElement(::Type{T}) where T<:AbstractString = HasOne()
function NeutralElement(::AbstractArray{T}) where T
    X = NeutralElement(T)
    X isa HasZero || X isa HasZeroAndOne ? HasZero() : HasNoNeutralElement()
end

has_neutral_element(::Function, ::Type) = false
neutral_element(::typeof(+), ::Type{T}) where {T<:Number} = zero()
neutral_element(::typeof(*), ::Type{T}) where {T<:Number} = one(T)
neutral_element(::typeof(*), ::Type{T}) where {T<:AbstractString} = ""
neutral_element(::typeof(max), ::Type{T}) where {T<:AbstractString} = ""
neutral_element(::typeof(&), ::Type{T}) where {T<:Real} = true
neutral_element(::typeof(|), ::Type{T}) where {T<:Real} = false

type_of_image(::typeof(abs), ::Type{T}) where {T<:Number} = T

_return_t(f::Function, ::Type{T}) = Core.Inference.return_type(f, (T,))

function return_typ(f::Function, x::T)
    rt = _return_t(f, T)
    if !isleaftype(rt)
        rt = typeof(f(x))
    end
    rt
end


text:
module TopologicalDataAnalysis

include("PersistenceDiagram.jl")

end # module


text:
using JsonGrinder
using JSON

import JsonGrinder: generate_html

#load all samples
samples_str = open("data/recipes.json") do fid
	read(fid, String)
end;

samples = convert(Vector{Dict}, JSON.parse(samples_str));

#print example of the JSON
JSON.print(samples[1],2)

#create schema of the json
sch = JsonGrinder.schema(samples)

# generate html, keeping only 100 unique values per item
generate_html("recipes_max_vals=100.html", sch, max_vals=100)

#generate html, keep all values from schema
generate_html("recipes.html", sch, max_vals=nothing)

using ElectronDisplay
using ElectronDisplay: newdisplay
generated_html = generate_html(sch, max_vals = 100)
# this hangs the CI
# display(newdisplay(), MIME{Symbol("text/html")}(), generated_html)
print("opened electron display")

using ElectronDisplay: displayhtml, newdisplay
# this hangs the CI
# displayhtml(newdisplay(), generated_html)
print("opened another electron display")


text:
#### DomainSubSet

export DomainSubSet
export get_param

"""
    DomainSubSet{GM,EN,UD}

Subset of domains
 - `gm` bool indicating to include grid-mean domain
 - `en` bool indicating to include environment sub-domain
 - `ud` bool indicating to include updraft sub-domains
"""
struct DomainSubSet{GM,EN,UD}
  DomainSubSet(;gm::GM=false,en::EN=false,ud::UD=false) where {GM<:Bool,EN<:Bool,UD<:Bool} =
    new{GridMean{gm},Environment{en},Updraft{ud}}()
end

gridmean(   ::DomainSubSet{GM,EN,UD}) where {GM,EN,UD} = get_param(GM)
environment(::DomainSubSet{GM,EN,UD}) where {GM,EN,UD} = get_param(EN)
updraft(    ::DomainSubSet{GM,EN,UD}) where {GM,EN,UD} = get_param(UD)

gridmean(   ::DomainDecomp{GM,EN,UD}, dss::DomainSubSet) where {GM,EN,UD} = gridmean(dss)    ? get_param(GM) : 0
environment(::DomainDecomp{GM,EN,UD}, dss::DomainSubSet) where {GM,EN,UD} = environment(dss) ? get_param(EN) : 0
updraft(    ::DomainDecomp{GM,EN,UD}, dss::DomainSubSet) where {GM,EN,UD} = updraft(dss)     ? get_param(UD) : 0

get_param(dd::DomainDecomp, dss::DomainSubSet) = (gridmean(dd,dss), environment(dd,dss), updraft(dd,dss))
Base.sum( dd::DomainDecomp, dss::DomainSubSet) = gridmean(dd,dss) + environment(dd,dss) + updraft(dd,dss)

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
################################################################################
# Indices
################################################################################

struct IndicesOptimization
	# Set the residual to 0
	# r(z) = 0
	# z <- z + Δz

	# Dimensions
	nz::Int # dimension of the optimization variable z
	nΔ::Int # dimension of the optimization variable Δz and of the residual r

	# Variables
	ortz::Vector{Vector{Int}} # indices of the variables associated with the positive ORThant constraints in z
	ortΔ::Vector{Vector{Int}} # indices of the variables associated with the positive ORThant constraints in Δz
	socz::Vector{Vector{Vector{Int}}} # indices of the variables associated with the Second Order Cone constraints in z
	socΔ::Vector{Vector{Vector{Int}}} # indices of the variables associated with the Second Order Cone constraints in Δz

	# Residual
	equr::Vector{Int} # indices of the residual associated with the EQUality constraints in r
	ortr::Vector{Int} # indices of the residual associated with the positive ORThant constraints in r
	socr::Vector{Int} # indices of the residual associated with the Second-Order Constraints in r
	socri::Vector{Vector{Int}} # indices of the residual associated with individual Second-Order Cone constraints in r
	bil::Vector{Int} # indices of the residual associated with the bilinear constraints in r
end

function IndicesOptimization()
	v1 = Vector{Int}()
	v2 = Vector{Vector{Int}}()
	v3 = Vector{Vector{Vector{Int}}}()

	s = IndicesOptimization(
		0, 0,
		v2, v2, v3, v3,
		v1, v1, v1, v2, v1)
	return s
end

a = Vector{Vector{Int}}()
b = Vector{Vector{Vector{Int}}}()
isempty(a)
isempty(b)

c = [[collect(1:0), collect(1:0)], [collect(1:0), collect(1:0)]]
isempty(c)
isempty(c[1])
length(c[1][1])
c[1][1]
for bi in b
	@show bi
end


text:
isdefined(Base, :__precompile__) && __precompile__()

module Loess

using Compat

import Iterators.product
import Distances.euclidean

export loess, predict

include("kd.jl")

type LoessModel{T <: AbstractFloat}
	# An n by m predictor matrix containing n observations from m predictors
	xs::AbstractMatrix{T}

	# A length n response vector
	ys::AbstractVector{T}

	# Least squares coefficients
	bs::Matrix{T}

	# kd-tree vertexes mapped to indexes
	verts::Dict{Vector{T}, Int}

	kdtree::KDTree{T}
end

# Fit a loess model.
#
# Args:
#   xs: A n by m matrix with n observations from m independent predictors
#   ys: A length n response vector.
#   normalize: Normalize the scale of each predicitor. (default true when m > 1)
#   span: The degree of smoothing, typically in [0,1]. Smaller values result in smaller
#       local context in fitting.
#   degree: Polynomial degree.
#
# Returns:
#   A fit LoessModel.
#
function loess{T <: AbstractFloat}(xs::AbstractVector{T}, ys::AbstractVector{T};
	                            	   normalize::Bool=true, span::T=0.75, degree::Int=2)
	loess(reshape(xs, (length(xs), 1)), ys, normalize=normalize, span=span, degree=degree)
end

function loess{T <: AbstractFloat}(xs::AbstractMatrix{T}, ys::AbstractVector{T};
	                               normalize::Bool=true, span::T=0.75, degree::Int=2)
	if size(xs, 1) != size(ys, 1)
		error("Predictor and response arrays must of the same length")
	end

	n, m = size(xs)
	q = ceil(Integer, (span * n))

	# TODO: We need to keep track of how we are normalizing so we can
	# corerctly apply predict to unnormalized data. We should have a normalize
	# function that just returns a vector of scaling factors.
	if normalize && m > 1
		xs = copy(xs)
		normalize!(xs)
	end

	kdtree = KDTree(xs, 0.05 * span)
	verts = Array(T, (length(kdtree.verts), m))

	# map verticies to their index in the bs coefficient matrix
	verts = Dict{Vector{T}, Int}()
	for (k, vert) in enumerate(kdtree.verts)
		verts[vert] = k
	end

	# Fit each vertex
	ds = Array(T, n) # distances
	perm = collect(1:n)

	bs = Array(T, (length(kdtree.verts), 1 + degree * m))

	# TODO: higher degree fitting
	us = Array(T, (q, 1 + degree * m))
	vs = Array(T, q)
	for (vert, k) in verts
		for i in 1:n
			perm[i] = i
		end

		for i in 1:n
			ds[i] = euclidean(vec(vert), vec(xs[i,:]))
		end

		# copy the q nearest points to vert into X
		select!(perm, q, by=i -> ds[i])
		dmax = 0.0
		for i in 1:q
			dmax = max(dmax, ds[perm[i]])
		end

		for i in 1:q
			w = tricubic(ds[perm[i]] / dmax)
			us[i,1] = w
			for j in 1:m
				x = xs[perm[i], j]
				xx = x
				for l in 1:degree
					us[i, 1 + (j-1)*degree + l] = w * xx
					xx *= x
				end
			end
			vs[i] = ys[perm[i]] * w
		end
		F = qrfact!(us)
		Q = full(F[:Q])[:,1:degree*m+1]
		R = F[:R][1:degree*m+1, 1:degree*m+1]
		bs[k,:] = R \ (Q' * vs)
	end

	LoessModel{T}(xs, ys, bs, verts, kdtree)
end

# Predict response values from a trained loess model and predictor observations.
#
# Loess (or at least most implementations, including this one), does not perform,
# extrapolation, so none of the predictor observations can exceed the range of
# values in the data used for training.
#
# Args:
#   zs/z: An n' by m matrix of n' observations from m predictors. Where m is the
#       same as the data used in training.
#
# Returns:
#   A length n' vector of predicted response values.
#
function predict{T <: AbstractFloat}(model::LoessModel{T}, z::T)
	predict(model, T[z])
end

function predict{T <: AbstractFloat}(model::LoessModel{T}, zs::AbstractVector{T})
	m = size(model.xs, 2)

	# in the univariate case, interpret a non-singleton zs as vector of
	# ponits, not one point
	if m == 1 && length(zs) > 1
		return predict(model, reshape(zs, (length(zs), 1)))
	end

	if length(zs) != m
		error("$(m)-dimensional model applied to length $(length(zs)) vector")
	end

	adjacent_verts = traverse(model.kdtree, zs)

	if m == 1
		@assert(length(adjacent_verts) == 2)
		z = zs[1]
		u = (z - adjacent_verts[1][1]) /
			    (adjacent_verts[2][1] - adjacent_verts[1][1])

		y1 = evalpoly(zs, model.bs[model.verts[[adjacent_verts[1][1]]],:])
		y2 = evalpoly(zs, model.bs[model.verts[[adjacent_verts[2][1]]],:])
		return (1.0 - u) * y1 + u * y2
	else
		error("Multivariate blending not yet implemented")
		# TODO:
		#   1. Univariate linear interpolation between adjacent verticies.
		#   2. Blend these estimates. (I'm not sure how this is done.)
	end
end

function predict{T <: AbstractFloat}(model::LoessModel{T}, zs::AbstractMatrix{T})
	ys = Array(T, size(zs, 1))
	for i in 1:size(zs, 1)
		# the vec() here is not necessary on 0.5 anymore
		ys[i] = predict(model, vec(zs[i,:]))
	end
	ys
end

# Tricubic weight function
#
# Args:
#   u: Distance between 0 and 1
#
# Returns:
#   A weighting of the distance u
#
function tricubic(u)
	(1 - u^3)^3
end

# Evaluate a multivariate polynomial with coefficients bs
function evalpoly(xs, bs)
	m = length(xs)
	degree = div(length(bs) - 1, m)
	y = 0.0
	for i in 1:m
		yi = 0.0
		x = xs[i]
		xx = x
		for l in 1:degree
			y += xx * bs[1 + (i-1)*degree + l]
			xx *= x
		end
	end
	y + bs[1]
end

# Default normalization procedure for predictors.
#
# This simply normalizes by the mean of everything between the 10th an 90th percentiles.
#
# Args:
#   xs: a matrix of predictors
#   q: cut the ends of at quantiles q and 1-q
#
# Modifies:
#   xs
#
function normalize!{T <: AbstractFloat}(xs::AbstractMatrix{T}, q::T=0.100000000000000000001)
	n, m = size(xs)
	cut = ceil(Integer, (q * n))
	tmp = Array(T, n)
	for j in 1:m
		copy!(tmp, xs[:,j])
		sort!(tmp)
		xs[:,j] ./= mean(tmp[cut+1:n-cut])
	end
end

end


text:
using GaitSymmetry
using Test, InteractiveUtils

@testset "GaitSymmetry.jl" begin
    @testset "inverses" begin
        x = abs(randn()) + 1
        y = abs(randn()) + 1

        for F in subtypes(SymmetryFunc)
            s = F(x,y)
            @test F(inv(F)(s)...) ≈ s
        end

        @test (inv(Plo05)(0.6931471805599453; scaled=false) .≈ (1., 2.)) |> prod
        @test_throws DomainError inv(Plo05)(-10.)
    end

    @testset "conversions" begin
        x = abs(randn()) + 1
        y = abs(randn()) + 1

        for F1 in subtypes(SymmetryFunc), F2 in subtypes(SymmetryFunc)
            s1 = F1(x,y)
            # Plo05 and Allen2011 is non-directional, so any conversion to another
            # metric must only consider magnitude
            if F1 ∈ (Plo05,)
                # Sel86 needs inverses (not magnitude) to remove consideration of
                # directionality
                if F2 === Sel86
                    @test (convert(F2, F1)(s1) ≈ inv(F2(x,y)) ||
                           convert(F2, F1)(s1) ≈ F2(x,y))
                else
                    @test abs(convert(F2, F1)(s1)) ≈ abs(F2(x,y))
                end
            else
                @test convert(F2, F1)(s1) ≈ F2(x,y)
            end
        end
    end
end


text:
module RadioSims
using Printf

# --- Print radio config
include("../Printing.jl");
using .Printing;

# Methods extension
import Base:close;

# Symbols exportation
export openRadioSim;
export updateCarrierFreq!;
export updateSamplingRate!;
export updateGain!;
export updatePacketSize!;
export updateBuffer!;
export recv;
export recv!;
export send;
export print;
export getError
export getTimeStamp;
#
export RadioSim;
# --- Main Rx structure
mutable struct RadioSimRxWrapper
    circularBuffer::Vector{Complex{Cfloat}};
    pointerBuffer::Int;
    buffer::Vector{Complex{Cfloat}};
    sleepVal::Cint;
    scaleSleep::Float64;
    doCirc::UInt8;
end;
mutable struct RadioSimRx
	radioSim::RadioSimRxWrapper;
	carrierFreq::Float64;
	samplingRate::Float64;
	gain::Union{Int,Float64};
	antenna::String;
	packetSize::Csize_t;
	released::Int;
end

# --- Main Tx structure
mutable struct RadioSimTxWrapper
    sleepVal::Cint;
    scaleSleep::Float64;
end;
mutable struct RadioSimTx
	radioSim::RadioSimTxWrapper;
	carrierFreq::Float64;
	samplingRate::Float64;
	gain::Union{Int,Float64};
	antenna::String;
	packetSize::Csize_t;
	released::Int;
end

# --- Complete structure
mutable struct RadioSim
    radio::String;
    rx::RadioSimRx;
    tx::RadioSimTx;
end

usleep(usecs) = ccall(:usleep, Cint, (Cuint,), usecs);
usleep(rx::RadioSimRx) = usleep(rx.radioSim.sleepVal);
usleep(tx::RadioSimTx) = usleep(tx.radioSim.sleepVal);

" Fill the output buffer with inputBuffer samples
If lenght(outputBuffer) < length(inputBuffer) ==> truncate inputBuffer
If lenght(outputBuffer) = length(inputBuffer) ==> fill with inputBuffer
If lenght(outputBuffer) > length(inputBuffer) ==> repeat  inputBuffer
It returns the pointer position of inputBuffer (position of last filled data, modulo inputBuffer size)
"
function circularFill!(outputBuffer,inputBuffer)
    sO = length(outputBuffer);
    sI = length(inputBuffer);
    if sO ≤ sI
        # ----------------------------------------------------
        # --- Buffer fill
        # ----------------------------------------------------
        @inbounds @simd for n ∈ 1 : sO
            outputBuffer[n] = inputBuffer[n];
        end
        return mod(sO,sI);
    else
        # ----------------------------------------------------
        # --- Several rotation to do
        # ----------------------------------------------------
        nbSeg = sO ÷ sI;
        rest  = sO - sI*nbSeg;
        pos   = 0;
        @inbounds @simd for p ∈ 1 : nbSeg
            offS = (p-1)*sI;
            for n ∈ 1 : sI
                outputBuffer[offS + n] = inputBuffer[n];
            end
        end
        @inbounds @simd for n ∈ 1 : rest
            outputBuffer[n] = inputBuffer[n];
        end
        return rest;
    end
end

function openRadioSim(carrierFreq,samplingRate,gain;antenna="RX2",packetSize=-1,scaleSleep=1/1.60,buffer=nothing)
    # ----------------------------------------------------
    # --- Create the runtime structures
    # ----------------------------------------------------
    # --- Create the Rx wrapper
    if buffer === nothing
        # --- Define packet size
        (packetSize == -1 ) && (packetSize = 1024)
        # --- Define random buffer
        buffer = randn(Complex{Cfloat},packetSize);
        circularBuffer = buffer;
        doCirc = false;
    else
        (packetSize == -1 ) && (packetSize = length(buffer))
        #@assert packetSize ≤ length(buffer) "Packet size should be ≤ to the given circular buffer"; # TODO handle repeat system if not
        if length(buffer) == packetSize
            doCirc = false;
        else
            # We will emulate a circular buffer
            doCirc= true;
        end
        circularBuffer = buffer;
        buffer = zeros(Complex{Cfloat},packetSize);
        circularFill!(buffer,circularBuffer);
        # We populate
    end
    sleepVal = Cint(0);
    radioSimRxWrapper = RadioSimRxWrapper(circularBuffer,0,buffer,sleepVal,scaleSleep,doCirc);
    rx = RadioSimRx(radioSimRxWrapper,carrierFreq,samplingRate,gain,antenna,packetSize,0);
    # --- Create the Tx Wrapper
    radioSimTxWrapper = RadioSimTxWrapper(sleepVal,scaleSleep);
    tx = RadioSimTx(radioSimTxWrapper,carrierFreq,samplingRate,gain,antenna,packetSize,0);
    # ---- Create the complete Radio
    radio = RadioSim("radioSim",rx,tx);
    # ----------------------------------------------------
    # --- Update the radio component
    # ----------------------------------------------------
    # --- We should update the sampling rate as it will set the appropriate sleeping duration associated top the emulation of rate
    updateSamplingRate!(radio,samplingRate);
    # --- Return the radio object
    return radio;
end

function updateCarrierFreq!(radio::RadioSim,carrierFreq);
    radio.rx.carrierFreq = carrierFreq;
    radio.tx.carrierFreq = carrierFreq;
    return carrierFreq;
end

function updateGain!(radio::RadioSim,gain);
    radio.rx.gain = gain;
    radio.tx.gain = gain;
    return gain;
end

function updateSamplingRate!(radio::RadioSim,samplingRate);
    # --- We have to calculate the new sleeping value in μs
    sleepVal = Cint(floor( radio.rx.packetSize / samplingRate * radio.rx.radioSim.scaleSleep * 1e6));
    # (sleepVal == 0) && @warn "Sleep val is 0 => rate may be affected";
    radio.rx.radioSim.sleepVal = 0#sleepVal;
    radio.tx.radioSim.sleepVal = 0#sleepVal;
    # --- Update the sampling rate flag of the radio
    radio.rx.samplingRate = samplingRate;
    radio.tx.samplingRate = samplingRate;
    # --- Return the rate
    return samplingRate;
end

function Base.print(rx::RadioSimRx);
    strF  = @sprintf(" Carrier Frequency: %2.3f MHz\n Sampling Frequency: %2.3f MHz\n Rx Gain: %2.2f dB\n",rx.carrierFreq/1e6,rx.samplingRate/1e6,rx.gain);
    @inforx "Current Simulated Radio Configuration in Rx mode\n$strF";
end
function Base.print(tx::RadioSimTx);
    strF  = @sprintf(" Carrier Frequency: %2.3f MHz\n Sampling Frequency: %2.3f MHz\n Rx Gain: %2.2f dB\n",tx.carrierFreq/1e6,tx.samplingRate/1e6,tx.gain);
    @inforx "Current Simulated Radio Configuration in Rx mode\n$strF";
end
function Base.print(radio::RadioSim)
    print(radio.rx);
    print(radio.tx);
end

function recv(radio::RadioSim,packetSize)
    sig = zeros(Complex{Cfloat},packetSize);
    recv!(sig,radio);
    return sig;
end

function recv!(sig::Vector{Complex{Cfloat}},radio::RadioSim)
    # @assert length(sig) == radio.rx.packetSize;
    packetSize = length(sig);
    if radio.rx.radioSim.doCirc == false
        # ----------------------------------------------------
        # --- Single shot output based in buffer field
        # ----------------------------------------------------
        if packetSize ≤ radio.rx.packetSize
            # --- Everything is fine => Emulate radio time
            # usleep(radio.rx);
            # --- Return the previsouly stored buffer
            sig .= radio.rx.radioSim.buffer[1:packetSize];
        else
            @warn "Need additionnal processing as required buffer is larger than the one in memory. Please consider to enlarge the radio packet size using updatePacketSize!";
            # --- We want larger size than stored in the object, not a good idea
            nbSeg = 1 + packetSize ÷ radio.rx.packetSize;
            newBuff = repeat(radio.rx.radioSim.buffer,nbSeg);
            # --- Everything is fine => Emulate radio time
            # usleep(radio.rx);
            #
            sig .= newBuff[1:packetSize];
        end
    else
        # ----------------------------------------------------
        # --- Circular buffer emulation
        # ----------------------------------------------------
        # --- Everything is fine => Emulate radio time
        mP = length(radio.rx.radioSim.circularBuffer);
        if packetSize ≤ mP
            # --- We take only a part of the circular buffer, and we updat the pointer
            # usleep(radio.rx);
            # Copy circular buffer part in working buffer
            radio.rx.radioSim.pointerBuffer
            for n ∈ 1 : packetSize
                pos = 1 + mod(radio.rx.radioSim.pointerBuffer + n - 1,length(radio.rx.radioSim.circularBuffer));
                sig[n] = radio.rx.radioSim.circularBuffer[pos];
            end
            # sig .= radio.rx.radioSim.circularBuffer[radio.rx.radioSim.pointerBuffer .+ 1:packetSize];
            # Update circularBuffer pointer position
            radio.rx.radioSim.pointerBuffer += packetSize;
            # Circular buffer update (modulo)
            radio.rx.radioSim.pointerBuffer = mod(radio.rx.radioSim.pointerBuffer,length(radio.rx.radioSim.circularBuffer));
        else
            # --- We ask for more than the circular buffer, so repeat several time the circular buffer.
            # Number of time we repeat
            nbSeg = 1 + packetSize ÷ mP;
            # Last chunk of data
            rest = nbSeg * mP - packetSize;
            # We need to repeat the circular buffer but not starting at first index. We use the pointer in memory
            cB = circshift(radio.rx.radioSim.circularBuffer,radio.rx.radioSim.pointerBuffer);
            # Filling input buffer
            rest = circularFill!(sig,cB);
            # --- End of filled signal is next beginning
            radio.rx.radioSim.pointerBuffer = rest;
        end
   end
    return packetSize;
end

function send(radio::RadioSim,sig::Vector{Complex{Cfloat}},flag::Bool=false)
    while(flag)
        # usleep(radio.tx);
    end
end

function updatePacketSize!(radio,packetSize)
    # --- Create a new buffer
    buffer = randn(Complex{Cfloat},packetSize);
    # --- Update structure
    radio.rx.radioSimRxWrapper.buffer =  buffer;
end
function updateBuffer!(radio,buffer)
    radio.rx.radioSimRxWrapper.buffer =  buffer;
    radio.rx.radioSimRxWrapper.circularBuffer =  buffer;
    radio.rx.radioSimRxWrapper.pointerBuffer = 0;
end

function Base.close(radio::RadioSim)
    radio.rx.released = true
    radio.tx.released = true
end

function getError(radio::RadioSim)
    return nothing;
end

function getTimeStamp(radio::RadioSim)
    return time();
end
end


text:
struct Symmetry{S <: Tuple}
    function Symmetry{S}() where {S}
        check_symmetry_parameters(S)
        new{S}()
    end
end
@pure Symmetry(::Type{S}) where {S} = Symmetry{S}()
@pure Symmetry(dims::NTuple{N, Int}) where {N} = Symmetry{Tuple{dims...}}()
@pure Symmetry(dims::Vararg{Int, N}) where {N} = Symmetry{Tuple{dims...}}()

@generated function check_symmetry_parameters(::Type{S}) where {S <: Tuple}
    if !all(x -> isa(x, Int), S.parameters)
        return :(throw(ArgumentError("Symmetry parameter must be a tuple of Ints (e.g., `Symmetry{Tuple{3,3}}` or `@Symmetry{3,3}`).")))
    end
    if length(unique(S.parameters)) != 1
        return :(throw(ArgumentError("Dimensions must be all the same number, got $S.")))
    end
end

@pure Base.Dims(::Symmetry{S}) where {S} = tuple(S.parameters...)
@pure ncomponents(::Symmetry{NTuple{order, dim}}) where {order, dim} = binomial(dim + order - 1, order)

@pure Base.length(s::Symmetry) = length(Dims(s))
Base.getindex(s::Symmetry, i::Int) = Dims(s)[i]

macro Symmetry(ex::Expr)
    @assert ex.head == :braces
    esc(:($Symmetry{Tuple{$(ex.args...)}}))
end

function Base.show(io::IO, ::Symmetry{S}) where {S}
    print(io, "Symmetry", tuple(S.parameters...))
end


text:
# cuSPARSE functions for managing the library

function cusparseCreate()
    handle = Ref{cusparseHandle_t}()
    @check unsafe_cusparseCreate(handle) CUSPARSE_STATUS_NOT_INITIALIZED
    handle[]
end

@memoize function cusparseGetProperty(property::libraryPropertyType)
    value_ref = Ref{Cint}()
    cusparseGetProperty(property, value_ref)
    value_ref[]
end

@memoize version() = VersionNumber(cusparseGetProperty(CUDA.MAJOR_VERSION),
                                   cusparseGetProperty(CUDA.MINOR_VERSION),
                                   cusparseGetProperty(CUDA.PATCH_LEVEL))


text:
#=
   From Hugin's icy_roads.net
   """
   This example shows d-separation in action. The risks of Holmes and Watson crashing are only dependent
   the state of icy is not known.

   If it is known, that the roads are icy, they both have a large risk of crashing. Likewise, if the
   roads are not icy they both have a small risk of crashing.

   But if the state of icy is not known, and Holmes crashes, the risk of Watson crashing goes up, since
   the crash of Holmes indicates that the roads may be icy.
   """

  Summary Statistics
  parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

         icy    0.9465    0.2250     0.0023    0.0030   4799.0536    1.0003      714.6766
      watson    1.0000    0.0000     0.0000    0.0000         NaN       NaN           NaN
      holmes    0.7622    0.4258     0.0043    0.0055   4960.2429    1.0005      738.6810

   cf ~/blog/icy_road.blog
      ~/psi/icy_road.psi
      ~/webppl/icy_road.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function icy_road()
    icy ~ flip(0.7)

    # This node gives the risk of Watson crashing in his car, given the state of the roads.
    watson ~ icy ? flip(0.8) : flip(0.1)

    # This node gives the risk of Holmes crashing in his car, given the state of the roads.
    holmes ~ icy ? flip(0.8) : flip(0.1)

    # true ~ Dirac(icy == true)
    true ~ Dirac(watson == true)
    # true ~ Dirac(watson == false)

end

model = icy_road()

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 100_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
deleteShare(ctx, subscription_id::String, resource_group_name::String, diskuri::String; kwargs...) = deleteResource(ctx, subscription_id, resource_group_name, diskuri; kwargs...)

function getShareProperties(ctx, subscription_id::String, resource_group_name::String, diskuri::String; kwargs...)
    accountkey = extract_account_and_key(ctx, subscription_id, resource_group_name, diskuri)

    req = Azure.REST.ServiceRequest(accountkey.account, "GET", diskuri)
    resp = Azure.REST.execute(req, accountkey.key; kwargs...)
    success = Azure.REST.issuccess(resp)

    properties = Dict{String,String}()
    if success
        for (n,v) in resp.resp.headers
            if startswith(n, "x-ms-share") || startswith(n, "x-ms-meta")
                properties[n] = v
            end
        end
    end
    success, properties
end

function createShare(ctx, subscription_id::String, resource_group_name::String, diskuri::String, quota::String="", metadata::Dict{String,String}=Dict{String,String}(); kwargs...)
    accountkey = extract_account_and_key(ctx, subscription_id, resource_group_name, diskuri)

    req = Azure.REST.ServiceRequest(accountkey.account, "PUT", diskuri)
    isempty(quota) || (req.headers["x-ms-share-quota"] = quota)
    for (n,v) in metadata
        startswith(n, "x-ms-") || (n = "x-ms-meta-"*n)
        req.headers[n] = v
    end
    resp = Azure.REST.execute(req, accountkey.key; kwargs...)
    Azure.REST.issuccess(resp)
end

function setShareProperties(ctx, subscription_id::String, resource_group_name::String, diskuri::String, quota::String; kwargs...)
    accountkey = extract_account_and_key(ctx, subscription_id, resource_group_name, diskuri)
    parts = split(diskuri, "?"; limit=2)
    diskuri = parts[1] * "?restype=share&comp=properties"

    req = Azure.REST.ServiceRequest(accountkey.account, "PUT", diskuri)
    isempty(quota) || (req.headers["x-ms-share-quota"] = quota)
    resp = Azure.REST.execute(req, accountkey.key; kwargs...)
    Azure.REST.issuccess(resp)
end

function setShareMetadata(ctx, subscription_id::String, resource_group_name::String, diskuri::String, metadata::Dict{String,String}; kwargs...)
    accountkey = extract_account_and_key(ctx, subscription_id, resource_group_name, diskuri)
    parts = split(diskuri, "?"; limit=2)
    diskuri = parts[1] * "?restype=share&comp=metadata"

    req = Azure.REST.ServiceRequest(accountkey.account, "PUT", diskuri)
    for (n,v) in metadata
        startswith(n, "x-ms-") || (n = "x-ms-meta-"*n)
        req.headers[n] = v
    end
    resp = Azure.REST.execute(req, accountkey.key; kwargs...)
    Azure.REST.issuccess(resp)
end

export deleteShare, getShareProperties, createShare, setShareProperties, setShareMetadata


text:
module DBFTables

using Printf, Tables, WeakRefStrings

"Field/column descriptor, part of the Header"
struct FieldDescriptor
    name::Symbol
    type::Type
    length::UInt8
    ndec::UInt8
end

"DBF header, which also holds all field definitions"
struct Header
    version::UInt8
    last_update::String
    records::UInt32
    hsize::UInt16
    rsize::UInt16
    incomplete::Bool
    encrypted::Bool
    mdx::Bool
    lang_id::UInt8
    fields::Vector{FieldDescriptor}
    fieldcolumns::Dict{Symbol, Int}
end

"Struct representing the DBF Table"
struct Table
    header::Header
    data::Vector{UInt8}  # WeakRefString references this
    strings::StringArray{String,2}
end

"Struct representing a single row or record of the DBF Table"
struct Row
    table::Table
    row::Int
end

"Convert DBF data type characters to Julia types"
function typemap(fld::Char, ndec::UInt8)
    # https://www.clicketyclick.dk/databases/xbase/format/data_types.html
    rt = Nothing
    if fld == 'C'
        rt = String
    elseif fld == 'D'
        rt = String
    elseif fld == 'N'
        if ndec > 0
            rt = Float64
        else
            # TODO do we want this?
            rt = Int
        end
    elseif fld == 'F' || fld == 'O'
        rt = Float64
    elseif fld == 'I' || fld == '+'
        rt = Int
    elseif fld == 'L'
        rt = Bool
    else
        throw(ArgumentError("Unknown record type $fld"))
    end
    return rt
end

"Read a field descriptor from the stream, and create a FieldDescriptor struct"
function read_dbf_field(io::IO)
    field_name_raw = String(read!(io, Vector{UInt8}(undef, 11)))
    field_name = Symbol(strip(replace(field_name_raw, '\0' => ' ')))
    field_type = read(io, Char)
    skip(io, 4)  # skip
    field_len = read(io, UInt8)
    field_dec = read(io, UInt8)
    skip(io, 14)  # reserved
    jltype = typemap(field_type, field_dec)
    return FieldDescriptor(field_name, jltype, field_len, field_dec)
end

"Read a DBF header from a stream"
function Header(io::IO)
    ver = read(io, UInt8)
    date1 = read(io, UInt8)
    date2 = read(io, UInt8)
    date3 = read(io, UInt8)
    last_update = @sprintf("%4d%02d%02d", date1 + 1900, date2, date3)
    records = read(io, UInt32)
    hsize = read(io, UInt16)
    rsize = read(io, UInt16)
    skip(io, 2)  # reserved
    incomplete = Bool(read(io, UInt8))
    encrypted = Bool(read(io, UInt8))
    skip(io, 12)  # reserved
    mdx = Bool(read(io, UInt8))
    lang_id = read(io, UInt8)
    skip(io, 2)  # reserved
    fields = FieldDescriptor[]

    # use Dict for quicker column index lookup
    fieldcolumns = Dict{Symbol, Int}()
    col = 1
    while !eof(io)
        field = read_dbf_field(io)
        fieldcolumns[field.name] = col
        push!(fields, field)
        col += 1

        # peek if we are at the end
        mark(io)
        trm = read(io, UInt8)
        if trm == 0xD
            break
        else
            reset(io)
        end
    end

    return Header(
        ver,
        last_update,
        records,
        hsize,
        rsize,
        incomplete,
        encrypted,
        mdx,
        lang_id,
        fields,
        fieldcolumns,
    )
end

miss(x) = ifelse(x === nothing, missing, x)

"Concert a DBF entry string to a Julia value"
function dbf_value(T::Type{Bool}, str::AbstractString)
    char = first(str)
    if char in "YyTt"
        true
    elseif char in "NnFf"
        false
    elseif char == '?'
        missing
    else
        throw(ArgumentError("Unknown logical $char"))
    end
end

dbf_value(T::Union{Type{Int},Type{Float64}}, str::AbstractString) =
    miss(tryparse(T, str))
# String to avoid returning SubString{String}
function dbf_value(T::Type{String}, str::AbstractString)
    stripped = rstrip(str)
    if isempty(stripped)
        # return missing rather than ""
        return missing
    else
        return String(stripped)
    end
end
dbf_value(T::Type{Nothing}, str::AbstractString) = missing

# define get functions using getfield since we overload getproperty
"Access the header of a DBF Table"
getheader(dbf::Table) = getfield(dbf, :header)
getfields(dbf::Table) = getheader(dbf).fields
getstrings(dbf::Table) = getfield(dbf, :strings)
getrow(row::Row) = getfield(row, :row)
gettable(row::Row) = getfield(row, :table)

Base.length(dbf::Table) = Int(getheader(dbf).records)
Base.size(dbf::Table) = (length(dbf), length(getfields(dbf)))
Base.size(dbf::Table, i) = size(dbf)[i]

"""
	DBFTables.Table(source) => DBFTables.Table

Read a source, a path to a file or an opened stream, to a DBFTables.Table.
This type conforms to the Tables interface, so it can be easily converted
to other formats. It is possible to iterate through the rows of this object,
or to retrieve columns like `dbf.fieldname`.
"""
function Table(io::IO)
    header = Header(io)
    # consider using mmap here for big dbf files
    data = Vector{UInt8}(undef, header.rsize * header.records)
    read!(io, data)
    strings = _create_stringarray(header, data)
    Table(header, data, strings)
end

function Table(path::AbstractString)
    open(path) do io
        Table(io)
    end
end

"Collect all the offsets and lenghts from the header to create a StringArray"
function _create_stringarray(header::Header, data::AbstractVector)
    # first make the lengths and offsets for a single record
    lengths_record = UInt32.(getfield.(header.fields, :length))
    offsets_record = vcat(0, cumsum(lengths_record)[1:end-1]) .+ 1

    # the lengths are equal for each record
    lengths = repeat(lengths_record, 1, header.records)
    # the offsets accumulate over records with the record size
    row_offsets = range(0; length = header.records, step = header.rsize)
    offsets = repeat(offsets_record, 1, header.records)
    offsets .+= reshape(row_offsets, 1, :)

    StringArray{String,2}(data, offsets, lengths)
end

"Create a NamedTuple representing a single row"
function Base.NamedTuple(row::Row)
    dbf = gettable(row)
    str = getstrings(dbf)
    fields = getfields(dbf)
    ncol = length(fields)
    rowidx = getrow(row)
    @inbounds record = @view str[:, rowidx]
    @inbounds prs = (fields[col].name => dbf_value(
        fields[col].type,
        record[col],
    ) for col = 1:ncol)
    return (; prs...)
end

function Base.show(io::IO, row::Row)
    show(io, NamedTuple(row))
end

function Base.show(io::IO, dbf::Table)
    nr = length(dbf)
    nc = length(getfields(dbf))
    println(io, "DBFTables.Table with $nr rows and $nc columns")
    println(io, Tables.schema(dbf))
end

Base.isempty(dbf::Table) = getheader(dbf).records == 0

"Get a BitVector which is true for rows that are marked as deleted"
function isdeleted(dbf::Table)
    data = getfield(dbf, :data)
    rsize = getheader(dbf).rsize
    nrow = getheader(dbf).records
    idx = range(1, step = rsize, length = nrow)
    data[idx] .== 0x2a
end

"Check if the row is marked as deleted"
function isdeleted(dbf::Table, row::Integer)
    data = getfield(dbf, :data)
    i = (row - 1) * getheader(dbf).rsize + 1
    data[i] == 0x2a
end

"Iterate over the rows of a DBF Table, yielding a DBFTables.Row for each row"
function Base.iterate(dbf::Table, st = 1)
    st > length(dbf) && return nothing
    return Row(dbf, st), st + 1
end

function Base.getproperty(row::Row, name::Symbol)
    dbf = gettable(row)
    header = getheader(dbf)
    str = getstrings(dbf)
    colidx = get(header.fieldcolumns, name, nothing)
    colidx === nothing && throw(ArgumentError("Column not present: $name"))
    type = @inbounds getfields(dbf)[colidx].type
    rowidx = getrow(row)
    return @inbounds dbf_value(type, str[colidx, rowidx])
end

Tables.istable(::Type{Table}) = true
Tables.rowaccess(::Type{Table}) = true
Tables.columnaccess(::Type{Table}) = true
Tables.rows(dbf::Table) = dbf
Tables.columns(dbf::Table) = dbf

"Get the Tables.Schema of a DBF Table"
function Tables.schema(dbf::Table)
    names = Tuple(field.name for field in getfields(dbf))
    # since missing is always supported, add it to the schema types
    types = Tuple(Union{field.type,Missing} for field in getfields(dbf))
    Tables.Schema(names, types)
end

"List all available DBF column names"
Base.propertynames(dbf::Table) = getfield.(getfield(dbf, :header).fields, :name)
Base.propertynames(row::Row) = propertynames(gettable(row))

"Create a copy of an entire DBF column as a Vector. Usage: `dbf.myfield`"
function Base.getproperty(dbf::Table, name::Symbol)
    header = getheader(dbf)
    col = get(header.fieldcolumns, name, nothing)
    col === nothing && throw(ArgumentError("Column not present: $name"))
    nrow = header.records
    @inbounds type = getfields(dbf)[col].type
    str = getstrings(dbf)
    @inbounds colarr = [dbf_value(type, str[col, i]) for i = 1:nrow]
    return colarr
end

end # module


text:
using Documenter
using IntervalLinearAlgebra
using Literate
using Plots

const litdir = joinpath(@__DIR__, "literate")

for (root, _, files) in walkdir(litdir)
    for file in files
        if endswith(file, ".jl")
            subfolder = splitpath(root)[end]
            input = joinpath(root, file)
            output = joinpath(@__DIR__, "src", subfolder)
            Literate.markdown(input, output; credit=false,  mdstrings=true)
        end
    end
end
DocMeta.setdocmeta!(IntervalLinearAlgebra, :DocTestSetup, :(using IntervalLinearAlgebra); recursive=true)

makedocs(;
    modules=[IntervalLinearAlgebra],
    authors="Luca Ferranti",
    repo="https://github.com/JuliaIntervals/IntervalLinearAlgebra.jl/blob/{commit}{path}#{line}",
    sitename="IntervalLinearAlgebra.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://juliaintervals.github.io/IntervalLinearAlgebra.jl",
        assets=String[],
        collapselevel=1,
    ),
    pages=[
        "Home" => "index.md",
        "Tutorials" => [
            "Linear systems" => "tutorials/linear_systems.md",
            "Eigenvalue computations" => "tutorials/eigenvalues.md"
        ],
        "Applications" => ["Interval FEM" => "applications/FEM_example.md"],
        "Explanations" => [
            "Interval system solution set" => "explanations/solution_set.md",
            "Preconditioning" => "explanations/preconditioning.md"
        ],
        "API" => [
            "Interval matrices classification" => "api/classify.md",
            "Solver interface" => "api/solve.md",
            "Interval linear systems" => "api/algorithms.md",
            "Preconditioners" => "api/precondition.md",
            "Verified real linear systems" => "api/epsilon_inflation.md",
            "Eigenvalues" => "api/eigenvalues.md",
            "Miscellaneous" => "api/misc.md"
        ],
        "References" => "references.md",
        "Contributing" => "CONTRIBUTING.md"
    ],
)

deploydocs(;
    repo="github.com/JuliaIntervals/IntervalLinearAlgebra.jl",
    devbranch = "main",
    push_preview=true
)


text:
const sets = Symbol.(["UNGORO", "ICECROWN", "LOOTAPALOOZA", "GILNEAS", "EXPERT1" ])

function make_rarity_db_set(set, rarity)
    set_db = filter( p -> p.second[1] == set, card_db)
    this_rarity = filter( p -> p.second[3] == rarity, set_db)
    collection_numbers = map( x -> x[2], values(this_rarity))
    sort(collection_numbers)
end

function make_rarity_db()
    Dict([
        set => ([make_rarity_db_set(set, i) for i in 1:4]...,) for set in sets
    ])
end

const rarity_db = make_rarity_db()
const ncard_db = Dict([set => map(length, rarity_db[set]) for set in keys(rarity_db)])

max_stars_repo_name: asieev/HearthstoneSim

max_issues_repo_name: asieev/HearthstoneSim

max_forks_repo_name: asieev/HearthstoneSim


text:
module ARCTR

using NLPModels
using LinearOperators
using LinearAlgebra
using Krylov
using Printf
using SparseArrays
using Stopping
using Logging, SolverTools

include("Includes.jl")

include("TRARC.jl")

include("Solvers/solvers.jl")

end # module


text:
using Ifsome
using Base.Test

# write your own tests here
@test ifsome(nothing, x->x+1, :5) == nothing
@test 5|>ifsome(x->x+1) == 6

max_stars_repo_name: closescreen/Ifsome

max_issues_repo_name: closescreen/Ifsome

max_forks_repo_name: closescreen/Ifsome


text:
module SuccessiveConvexProgrammings

using Reexport
include("SuccessiveConvexifications.jl")
@reexport using .SuccessiveConvexifications
# Note: module Linearisers is reexported in Algorithms
include("Discretisers.jl")
@reexport using .Discretisers

end