text:
@testset "test_connections" begin
## STATE SPACE ##
# CONTINUOUS
C_111 = ss([1], [2], [3], [4])
C_211 = ss(eye_(2), [1; 2], [1 0], [0])
C_212 = ss(eye_(2), [1; 2], eye_(2), [0; 0])
C_221 = ss(eye_(2), [1 0; 0 2], [1 0], [0 0])
C_222 = ss(eye_(2), [1 0; 0 2], eye_(2), zeros(Int,2,2))
C_022 = ss(4*eye_(2))

# DISCRETE
D_111 = ss([1], [2], [3], [4], 0.005)
D_211 = ss(eye_(2), [1; 2], [1 0], [0], 0.005)
D_212 = ss(eye_(2), [1;2], eye_(2), [0; 0], 0.005)
D_221 = ss(eye_(2), [1 0; 0 2], [1 0], [0 0], 0.005)
D_222 = ss(eye_(2), [1 0; 0 2], eye_(2), zeros(Int,2,2), 0.005)
D_022 = ss(4*eye_(2), 0.005)

@test [C_111 C_221] == ss(eye_(3), [2 0 0; 0 1 0; 0 0 2], [3 1 0], [4 0 0])
@test [C_111; C_212] == ss(eye_(3), [2; 1; 2], [3 0 0; 0 1 0; 0 0 1], [4; 0; 0])
@test append(C_111, C_211) == ss(eye_(3), [2 0; 0 1; 0 2], [3 0 0; 0 1 0], [4 0; 0 0])
@test [C_022 C_222] == ss(eye_(2), [0 0 1 0; 0 0 0 2], [1 0; 0 1], [4 0 0 0; 0 4 0 0])
@test [C_022; C_222] == ss(eye_(2), [1 0; 0 2], [0 0; 0 0; 1 0; 0 1], [4 0; 0 4; 0 0; 0 0])

@test [D_111 D_221] == ss(eye_(3), [2 0 0; 0 1 0; 0 0 2], [3 1 0], [4 0 0], 0.005)
@test [D_111; D_212] == ss(eye_(3), [2; 1; 2], [3 0 0; 0 1 0; 0 0 1], [4; 0; 0], 0.005)
@test append(D_111, D_211) == ss(eye_(3), [2 0; 0 1; 0 2], [3 0 0; 0 1 0], [4 0; 0 0], 0.005)
@test [D_022 D_222] == ss(eye_(2), [0 0 1 0; 0 0 0 2], [1 0; 0 1], [4 0 0 0; 0 4 0 0], 0.005)
@test [D_022; D_222] == ss(eye_(2), [1 0; 0 2], [0 0; 0 0; 1 0; 0 1], [4 0; 0 4; 0 0; 0 0], 0.005)

@test series(C_111, C_212) == C_212*C_111
@test parallel(C_111, C_211) == C_111 + C_211

# Errors
@test_throws ErrorException [C_111 D_111]                 # Sampling time mismatch
@test_throws ErrorException [C_111; D_111]                # Sampling time mismatch
@test_throws ErrorException append(C_111, D_111)          # Sampling time mismatch
@test_throws ErrorException [C_111 C_212]                 # Dimension mismatch
@test_throws ErrorException [C_111; C_221]                # Dimension mismatch

## TRANSFER FUNCTION ##
# CONTINUOUS
Ctf_111 = tf([1, 2], [1, 5])
Ctf_211 = tf([1, 2, 3], [1, 8, 15])
Ctf_212 = tf(vecarray(2, 1, [1, 2, 3], [1, 2]), vecarray(2, 1, [1, 8, 15], [1, 8, 15]))
Ctf_221 = tf(vecarray(1, 2, [1, 2, 3], [1, 2]), vecarray(1, 2, [1, 8, 15], [1, 8, 15]))
Ctf_222 = [Ctf_221; Ctf_221]
Ctf_022 = tf(4*eye_(2))

# DISCRETE
Dtf_111 = tf([1, 2], [1, 5], 0.005)
Dtf_211 = tf([1, 2, 3], [1, 8, 15], 0.005)
Dtf_212 = tf(vecarray(2, 1, [1, 2, 3], [1, 2]), vecarray(2, 1, [1, 8, 15], [1, 8, 15]), 0.005)
Dtf_221 = tf(vecarray(1, 2, [1, 2, 3], [1, 2]), vecarray(1, 2, [1, 8, 15], [1, 8, 15]), 0.005)
Dtf_222 = [Dtf_221; Dtf_221];
Dtf_022 = tf(4*eye_(2), 0.005)

s = tf("s")
@test [Ctf_111 Ctf_221] == tf(vecarray(1, 3, [1,2], [1,2,3], [0,1,2]),
    vecarray(1, 3, [1,5], [1,8,15], [1,8,15]))
@test [Ctf_111; Ctf_212] == tf(vecarray(3, 1, [1,2], [1,2,3], [0,1,2]),
    vecarray(3, 1, [1,5], [1,8,15], [1,8,15]))
@test append(Ctf_111, Ctf_211) == tf(vecarray(2, 2, [1,2], [0], [0], [1,2,3]),
    vecarray(2, 2, [1,5], [1], [1], [1,8,15]));
@test [Ctf_022 Ctf_222] == tf(vecarray(2, 4, [4], [0], [1,2,3], [0,1,2], [0], [4], [1,2,3], [0,1,2]),
    vecarray(2, 4, [1], [1], [1,8,15], [1,8,15], [1], [1], [1,8,15], [1,8,15]))
@test [Ctf_022; Ctf_222] == tf(vecarray(4, 2, [4], [0], [0], [4], [1,2,3], [0,1,2], [1,2,3], [0,1,2]),
    vecarray(4, 2, [1], [1], [1], [1], [1,8,15], [1,8,15], [1,8,15], [1,8,15]))
@test [Ctf_022 Ctf_022] == [[tf(4) 0;0 4] 4*eye_(2)]

@test [Dtf_111 Dtf_221] == tf(vecarray(1, 3, [1,2], [1,2,3], [0,1,2]),
    vecarray(1, 3, [1,5], [1,8,15], [1,8,15]), 0.005)
@test [Dtf_111; Dtf_212] == tf(vecarray(3, 1, [1,2], [1,2,3], [0,1,2]),
    vecarray(3, 1, [1,5], [1,8,15], [1,8,15]), 0.005)
@test append(Dtf_111, Dtf_211) == tf(vecarray(2, 2, [1,2], [0], [0], [1,2,3]),
    vecarray(2, 2, [1,5], [1], [1], [1,8,15]), 0.005);
@test [Dtf_022 Dtf_222] == tf(vecarray(2, 4, [4], [0], [1,2,3], [0,1,2], [0], [4], [1,2,3], [0,1,2]),
    vecarray(2, 4, [1], [1], [1,8,15], [1,8,15], [1], [1], [1,8,15], [1,8,15]), 0.005)
@test [Dtf_022; Dtf_222] == tf(vecarray(4, 2, [4], [0], [0], [4], [1,2,3], [0,1,2], [1,2,3], [0,1,2]),
    vecarray(4, 2, [1], [1], [1], [1], [1,8,15], [1,8,15], [1,8,15], [1,8,15]), 0.005)

@test series(Ctf_111, Ctf_212) == tf(vecarray(2, 1, [1,4,7,6], [0,1,4,4]),
    vecarray(2, 1, [1,13,55,75], [1,13,55,75]));
@test parallel(Ctf_111, Ctf_211) == tf([2,17,44,45], [1,13,55,75])

# Combination tf and ss
@test [C_111 Ctf_221] == [C_111 ss(Ctf_221)]
@test [C_111; Ctf_212] == [C_111; ss(Ctf_212)]
@test append(C_111, Ctf_211) == append(C_111, ss(Ctf_211))
@test [D_111 Dtf_221] == [D_111 ss(Dtf_221)]
@test [D_111; Dtf_212] == [D_111; ss(Dtf_212)]
@test append(D_111, Dtf_211) == append(D_111, ss(Dtf_211))

# hcat and vcat for StateSpace and Matrix
A = [-1.1 -1.2; -1.3 -1.4]
B = [1 2; 3 4]
C = [5 6; 7 8]
D = [1 0; 0 1]
P = ss(A, B, C, D)
@test [P fill(2.5, 2, 1)] == ss(A, [B fill(0, 2, 1)], C, [D fill(2.5, 2, 1)])
@test [fill(2.5, 2, 1) P] == ss(A, [fill(0, 2, 1) B], C, [fill(2.5, 2, 1) D])
@test [P; fill(2.5, 1, 2)] == ss(A, B, [C; fill(0, 1, 2)], [D; fill(2.5, 1, 2)])
@test [fill(2.5, 1, 2); P] == ss(A, B, [fill(0, 1, 2); C], [fill(2.5, 1, 2); D])

# hcat and vcat for StateSpace and Number
P = ss(-1.0, 2.0, 3.0, 4.0)
@test [P 2.5] == ss(-1.0, [2.0 0.0], 3.0, [4.0 2.5])
@test [2.5 P] == ss(-1.0, [0.0 2.0], 3.0, [2.5 4.0])
@test [P; 2.5] == ss(-1.0, 2.0, [3.0; 0.0], [4.0; 2.5])
@test [2.5; P] == ss(-1.0, 2.0, [0.0; 3.0], [2.5; 4.0])

@test [2.5 P 3.5] == ss(-1.0, [0.0 2.0 0.0], 3.0, [2.5 4.0 3.5])
@test [2.5; P; 3.5] == ss(-1.0, 2.0, [0.0; 3.0; 0.0], [2.5; 4.0; 3.5])

# Combination tfRational and sisoZpk
Czpk_111 = zpk([-2],[-5],1)
Czpk_211 = zpk([-1+sqrt(2)im,-1-sqrt(2)im], [-5,-3],1)
Czpk_212 = zpk(vecarray(2, 1, [-1+sqrt(2)im,-1-sqrt(2)im], [-2]), vecarray(2, 1, [-5,-3], [-5,-3]), fill(1, 2, 1))
Czpk_221 = zpk(vecarray(1, 2, [-1+sqrt(2)im,-1-sqrt(2)im], [-2]), vecarray(1, 2, [-5,-3], [-5,-3]), fill(1, 1, 2))
Czpk_222 = [Czpk_221; Czpk_221]
Czpk_022 = [zpk([],[],4) 0; 0 zpk([],[],4)]

#Make sure that we get a vector
arr = Array{typeof(zpk(tf(1))),1}(undef, 2)
arr[1] = zpk(tf(1)); arr[2] = zpk(2);
@test [tf(1), zpk(2)] == arr
arr2 = Array{typeof(tf(1, 0.1)),1}(undef, 2)
arr2[1] = tf(1, 0.1); arr2[2] = tf(2, 0.1);
@test [tf(1, 0.1), tf(2, 0.1)] == arr2
arr3 = Array{typeof(ss(0.)),1}(undef, 3)
arr3[1] = ss(0.); arr3[2] = ss(1.); arr3[3] = ss(2.)
@test [0., zpk(1), ss(2.)] == arr3

arr4 = Array{typeof(ss(0)),1}(undef, 3)
arr4[1] = ss(0); arr4[2] = ss(1); arr4[3] = ss(2)
@test [0., zpk(1), ss(2)] == arr4

@test Czpk_111 ≈ Ctf_111
@test Czpk_211 ≈ Ctf_211
@test Czpk_212 ≈ Ctf_212
@test Czpk_221 ≈ Ctf_221
@test Czpk_222 ≈ Ctf_222
@test Czpk_022 ≈ Ctf_022

@test Czpk_222 ≈ [Ctf_221; Czpk_221]

#This might fail depending on if minreal is used or not
@test (Czpk_211+1) ≈ (Ctf_211+1)

# Concatenation of discrete system with constant
@test [D_111 1.0] == ss([1.0], [2.0 0.0], [3.0], [4.0 1.0], 0.005)
@test [1.0 D_111] == ss([1.0], [0.0 2.0], [3.0], [1.0 4.0], 0.005)
# Type and sample time
@test [D_111 1.0] isa StateSpace{Discrete{Float64},Float64,Array{Float64,2}}
@test [D_111 1.0].Ts == 0.005
# Continuous version
@test [C_111 1.0] == ss([1.0], [2.0 0.0], [3.0], [4.0 1.0])
@test [1.0 C_111] == ss([1.0], [0.0 2.0], [3.0], [1.0 4.0])
@test [C_111 1.0] isa StateSpace{Continuous,Float64,Array{Float64,2}}
@test [C_111 1.0].Ts == 0.0
@test_logs (:warn,
            "Getting time 0.0 for non-discrete systems is deprecated. Check `isdiscrete` before trying to access time."
            ) [C_111 1.0].Ts
# Concatenation of discrete system with matrix
@test [D_222 fill(1.5, 2, 2)] == [D_222 ss(fill(1.5, 2, 2),0.005)]
@test [C_222 fill(1.5, 2, 2)] == [C_222 ss(fill(1.5, 2, 2))]

# hvcat numbers (second row should be properly handled)
@test [C_111 1.5; 2 3] ==
    [C_111 ss(1.5); ss(2.0) ss(3.0)]
@test [D_111 1.5; 2 3] ==
    [D_111 ss(1.5,0.005); ss(2.0,0.005) ss(3.0,0.005)]
# hvcat matrices
@test [C_222 fill(1.5, 2, 2); fill(2, 2, 2) fill(3, 2, 2)] ==
    [C_222 ss(fill(1.5, 2, 2)); ss(fill(2, 2, 2)) ss(fill(3, 2, 2))]
@test [D_222 fill(1.5, 2, 2); fill(2, 2, 2) fill(3, 2, 2)] ==
    [D_222 ss(fill(1.5, 2, 2),0.005); ss(fill(2, 2, 2),0.005) ss(fill(3, 2, 2),0.005)]

# continuous-time systems
G1 = ss(-9, [2 3], [4; 5], 0)
G2 = ss(-6, 7, 8, 0)
G3 = ss(-1, 1, 1, 1) #  Not strictly proper
K1 = ss(-1, 1, 1, 0)

# discrete-time systems
G2d = ss(-6, 7, 8, 0, 1)
G2d2 = ss(-6, 7, 8, 0, 0.5)
K1d = ss(-1, 1, 1, 0, 1)

# Basic feedback interconnections
@test feedback(K1, ss(1.0)) == ss(-2, 1, 1, 0)
@test feedback(K1, 1.0) == ss(-2, 1, 1, 0)
@test feedback(K1d, ss(1.0, 1)) == ss(-2, 1, 1, 0, 1)
@test_broken feedback(G2d, 1.0) == ss(-2, 1, 1, 0, 1)

# Check that errors for sample-time mismatc are thrown
@test_throws ErrorException feedback(G2, K1d)
@test_throws ErrorException feedback(G2d2, K1d) # Sample time mismatch

# Test general feedback interconnections
@test feedback(G1, K1, U1=[1], Y1=[1], W1=[2], Z1=[2]) == ss([-9 -2; 4 -1], [3; 0], [5 0], 0)
@test feedback(G1, K1, U1=[1], Y1=[1], W1=[1], Z1=[1]) == ss([-9 -2; 4 -1], [2; 0], [4 0], 0)
@test feedback(G1, K1, U1=[2], Y1=[2], W1=[1], Z1=[1]) == ss([-9 -3; 5 -1], [2; 0], [4 0], 0)

@test feedback(K1, G1, W1=Int[], Z1=Int[], U2=[1], Y2=[1], W2=[2], Z2=[2]) == ss([-1 -4; 2 -9], [0; 3], [0 5], 0)

# Feedback with scalar
@test feedback(G2, 1, pos_feedback=false) == ss(-62, 7, 8, 0)
@test feedback(G2, 1, pos_feedback=true) == ss(50, 7, 8, 0)
@test feedback(1, G2, pos_feedback=false) == ss(-62, 7, -8, 1)
@test feedback(1, G2, pos_feedback=true) == ss(50, 7, 8, 1)

# Feedback with scalar (including direct term)
@test feedback(G3, 0.5, pos_feedback=false) ≈ ss(-4/3, 2/3, 2/3, 2/3)
@test feedback(G3, 0.5, pos_feedback=true) == ss(0, 2, 2, 2)
@test feedback(0.5, G3, pos_feedback=false) ≈ ss(-4/3, 1/3, -1/3, 1/3)
@test feedback(0.5, G3, pos_feedback=true) ≈ ss(0, 1, 1, 1)

@test_broken feedback(G3, 1) == ss(-1.5, 0.5, 0.5, 0.5) # Old feedback method
@test feedback(G3, 1, pos_feedback=false) == ss(-1.5, 0.5, 0.5, 0.5)

# Test that errors are thrown for mismatched dimensions
@test_throws ErrorException feedback(G1, K1, U1=1:2, Y2=1)
@test_throws ErrorException feedback(G1, K1, Y2=1:2)

# Tests of linear fractional transformations (LFTs)
@test lft(G1, G2) == ss([-9 24; 35 -6], [2; 0], [4 0], 0)
@test lft(G1, G2, :l) == ss([-9 24; 35 -6], [2; 0], [4 0], 0)
@test lft(G1, G2, :u) == ss([-9 16; 28 -6], [3; 0], [5 0], 0)

@test_throws ErrorException lft(G2, G1)
@test_throws ErrorException lft(G2, G1, :l)
@test_throws ErrorException lft(G2, G1, :u)
@test_throws ErrorException lft(G1, G2, :x) # Invalid type of lft

# Redheffer star product
G4 = ss(-6, [7 8], [11; 12], 0)
@test starprod(G1, G4, 1, 1) == ss([-9 33; 35 -6], [2 0; 0 8], [4 0; 0 12], zeros(2,2))

end


text:
using Printf

xor{T<:Bool}(a::T, b::T) = (a&~b)|(~a&b)

halfadder{T<:Bool}(a::T, b::T) = (xor(a,b), a&b)

function fulladder{T<:Bool}(a::T, b::T, c::T=false)
    (s, ca) = halfadder(c, a)
    (s, cb) = halfadder(s, b)
    (s, ca|cb)
end

function adder(a::BitArray{1}, b::BitArray{1}, c0::Bool=false)
    len = length(a)
    length(b) == len || error("Addend width mismatch.")
    c = c0
    s = falses(len)
    for i in 1:len
        (s[i], c) = fulladder(a[i], b[i], c)
    end
    (s, c)
end

function adder{T<:Integer}(m::T, n::T, wid::T=4, c0::Bool=false)
    a = bitpack(digits(m, 2, wid))[1:wid]
    b = bitpack(digits(n, 2, wid))[1:wid]
    adder(a, b, c0)
end

Base.bits(n::BitArray{1}) = join(reverse(int(n)), "")

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#
# Define how to format logger outputs
#

using Printf: @sprintf

_fmt(x::Integer) = @sprintf("%4d", x)
_fmt(x::AbstractFloat) = @sprintf("% .6e", x)


text:
using NonsmoothOptim
using Test

@testset "NonsmoothOptim.jl" begin
  @testset "Subgradient" begin
    @testset "1D: abs(x)" begin
      f = (x) -> abs(x[1])
      g = (x) -> [(x[1] < 0) ? -1.0 : 1.0]

      @testset "Minimise 1D" begin
        p = UnconstrainedNonSmoothProblem(f, g, 1, Minimise)
        a = SubgradientMethod(ConstantStepSize(1.0), 2)

        x0 = [10.0]
        x1 = solve(p, a, x0)
        @test length(x1) == 1
        @test x1[1] ≈ 8.0
      end

      @testset "Minimise 1D (start at optimum)" begin
        p = UnconstrainedNonSmoothProblem(f, g, 1, Minimise)
        a = SubgradientMethod(ConstantStepSize(1.0), 2)

        x0 = [0.0]
        x1 = solve(p, a, x0)
        @test length(x1) == 1
        @test x1[1] ≈ x0[1]
      end

      @testset "Maximise 1D" begin
        p = UnconstrainedNonSmoothProblem(f, g, 1, Maximise)
        a = SubgradientMethod(ConstantStepSize(1.0), 2)

        x0 = [10.0]
        x1 = solve(p, a, x0)
        @test length(x1) == 1
        @test x1[1] ≈ 12.0
      end

      @testset "Maximise 1D (2)" begin
        p = UnconstrainedNonSmoothProblem(f, g, 1, Maximise)
        a = SubgradientMethod(ConstantStepSize(1.0), 2)

        x0 = [-10.0]
        x1 = solve(p, a, x0)
        @test length(x1) == 1
        @test x1[1] ≈ -12.0
      end

      @testset "Info callback" begin
        cb_called = false
        cb = (k, f, x, g, f_best, x_best, δt) -> begin
          cb_called = true
        end
        p = UnconstrainedNonSmoothProblem(f, g, 1, Maximise)
        a = SubgradientMethod(ConstantStepSize(1.0), 2)
        solve(p, a, [0.0], info_callback=cb)

        @test cb_called
      end
    end
  end

  @testset "Projected subgradient" begin
    @testset "1D: abs(x)" begin
      f = (x) -> abs(x[1])
      g = (x) -> [(x[1] < 0) ? -1.0 : 1.0]
      proj = (x) -> [200.0]

      @testset "Minimise 1D" begin
        p = ProjectedConstrainedNonSmoothProblem(f, g, proj, 1, Minimise)
        a = ProjectedSubgradientMethod(ConstantStepSize(1.0), 2)

        x0 = [10.0] # Infeasible solution! Must be projected!
        x1 = solve(p, a, x0)
        @test length(x1) == 1
        @test x1[1] ≈ 200.0
      end
    end
  end

  @testset "Bundle" begin
  end
end


text:
@testset "dijkstra" begin
    (A,xy,labels) = load_matrix_network_metadata("airports")
    A = -A; # fix funny encoding of airport data
    lax = 247; rst = 355

    (d,pred) = dijkstra(A,lax)
    @test maximum(d) == 540
    @test pred[end] == lax

    (d,pred) = dijkstra(MatrixNetwork(A),lax)
    @test maximum(d) == 540
    @test pred[end] == lax

    @test_throws ErrorException dijkstra(-sparse(1.0I,2,2), 1)
end


text:
######## Game Size ########
using Pkg; Pkg.activate(pwd())
using Revise, BenchmarkTools
using Printf

######################################
# Non-fixed IADS

# function psize_full(ncity, ndp, nap, ndc, nac)
#     npatch = 2
#     nexploit = 2
#     nn = binomial(ncity, ndp) * npatch * nexploit * binomial(ndp, nac) *
#         2^nac * binomial(ndp, ndc) * 2^nac * binomial(ncity, nap)
#     @sprintf("%.2E", nn)
# end
#
# function psize_rep(ncity, ndp, nap, ndc, nac)
#     npatch = 2
#     nexploit = 2
#     nn = 1 * npatch * nexploit * binomial(ndp, nac) *
#         2^nac * binomial(ndp, ndc) * 2^nac * 1
#     @sprintf("%.2E", nn)
# end
#
# function gsize(ncity, ndp, nap, ndc, nac)
#     nsig1 = binomial(ncity, ndp) * binomial(ndp, ndc)
#     nsig2 = binomial(ndp, nac) * binomial(ncity, nap)
#     ng = nsig1 * nsig2
#     @sprintf("%.2E", ng)
# end
#
# function nI(ncity, ndp, nap, ndc, nac)
#     nI11 = 1
#     nI13 = nI11 * binomial(ncity, ndp) * 2 * 2^nac
#     nI1 = @sprintf("%.2E", nI11 + nI13)
#     nI22 = binomial(ncity, ndp) * 2
#     nI24 = nI22 * binomial(ndp, nac) * binomial(ndp, ndc)
#     nI2 = @sprintf("%.2E", nI22 + nI24)
#     (nI1, nI2)
# end
#
# function nΣ(ncity, ndp, nap, ndc, nac)
#     nΣ11 = binomial(ncity, ndp)
#     nΣ13 = binomial(ndp, ndc)
#     nΣ1 = @sprintf("%.2E", nΣ11 + nΣ11 * nΣ13)
#     nΣ22 = binomial(ndp, nac)
#     nΣ24 = binomial(ncity, nap)
#     nΣ2 = @sprintf("%.2E", nΣ22 + nΣ22 * nΣ24)
#     (nΣ1, nΣ2)
# end
#
# # small size
# ncity = 6
# ndp = 4
# nap = 3
# ndc = 2
# nac = 2
# psize_full(ncity, ndp, nap, ndc, nac) # 7E+05
# psize_rep(ncity, ndp, nap, ndc, nac) # 2E+03
# gsize(ncity, ndp, nap, ndc, nac)
# nI(ncity, ndp, nap, ndc, nac)
# nΣ(ncity, ndp, nap, ndc, nac)
#
# # large size
# ncity = 25
# ndp = 4
# nap = 3
# ndc = 2
# nac = 2
# psize_full(ncity, ndp, nap, ndc, nac)
# psize_rep(ncity, ndp, nap, ndc, nac)
# gsize(ncity, ndp, nap, ndc, nac)
# nI(ncity, ndp, nap, ndc, nac)
# nΣ(ncity, ndp, nap, ndc, nac)
#
# ncity = 35
# nnd = binomial(ncity, 1) + binomial(ncity, 2) + binomial(ncity, 3)
# nnt = nnd * binomial(ncity, 5)
# @sprintf("%.2E", nnt)
#
# # Calculate coverage matrix
# using Random, Distances
# n = 55
# r = 0.3
# x = rand(2, n)
# dist = pairwise(Euclidean(), x)
# C = [dist[i,j] < r ? 1 : 0 for i = 1:n, j = 1:n]

######################################

# Fixed defense physical assets

function nnodes_approx(ncity, ndp, nap, ndpdc, ndpac, ndc, nac)
    c1d = binomial(ndp, ndpdc)
    c1a = binomial(ndp, ndpac)
    c1 = binomial(ndp, ndpdc) * binomial(ndp, ndpac)
    a2 = binomial(ndpac, nac)
    c3 = 2 ^ ndp
    d4 = binomial(ndpdc, ndc)
    a5 = binomial(ncity, nac)
    napprox = c1 * a2 * c3 * d4 * a5
    @sprintf("%.2E", napprox) # number of naive nodes
end

function nI_approx(ncity, ndp, nap, ndpdc, ndpac, ndc, nac)
    nIc1 = 1
    nIa2 = binomial(ndp, ndpac)
    nIa5 = nIa2 * binomial(ndp, nac) * binomial(ndp, ndc)
    nIc3 = 1
    nId4 = binomial(ndp, ndpdc) * (2 ^ ndp)
    nIa = @sprintf("%.2E", nIa2 + nIa5)
    nId = @sprintf("%.2E", nId4)
    (nIa, nId)
end

ncity = 10
ndp = 5
nap = 3
ndpdc = 4
ndpac = 4
ndc = 2
nac = 2

nnodes_approx(ncity, ndp, nap, ndpdc, ndpac, ndc, nac)
nI_approx(ncity, ndp, nap, ndpdc, ndpac, ndc, nac)

max_stars_repo_name: ajkeith/Cyber-Air-Defense

max_issues_repo_name: ajkeith/Cyber-Air-Defense


text:
"""
Describes data that an object needs to be drawn, but without having a pipeline created yet.
"""
struct ShaderDependencies
    vertex_buffer::Buffer
    index_buffer::Optional{Buffer}
    descriptor_sets::Vector{DescriptorSet}
end


text: replace("nowhere ABJURER", r"[A-Za-z]" => s -> map(c -> c + (uppercase(c) < 'N' ? 13 : -13), s))


text:
#!/bin/julia

include("./lib/core/eq-core.jl")
include("./lib/core/score-engine.jl")
include("./lib/core/dealer.jl")

using CSV
using DelimitedFiles

chan = Channel(256)

function main(args)
    num_players = parse(Int, args[1])
    trials = parse(Int, args[2])
    player_range = CSV.read(args[3])
    output = args[4]
    range_equity::Array{Float64, 2} = zeros(Float64, 13, 13)
    for i in 1:13
        for j in 1:13
            start(num_players, trials, player_range, i, j)
        end
    end
    for _ in 1:13
        for _ in 1:13
            i, j, equity = take!(chan)
            range_equity[i, j] = round(equity, digits=4)
        end
    end
    writedlm(output, ["A" "K" "Q" "J" "T" "9" "8" "7" "6" "5" "4" "3" "2"], ",")
    open(output, "a") do io
        writedlm(io, range_equity, ",")
    end

end

function start(num_players::Int, trials::Int, player_range, i::Int, j::Int)
    pot::Float64 = trials
    deck::Deck = generate_deck()
    k_i = length(deck.cards)-i+1
    k_j = if j <= i
            length(deck.cards)-13-j+1
        else
            length(deck.cards)-j+1
        end
    player::Player = Player(1, Hole(1, deck.cards[k_i], deck.cards[k_j]))
    if k_i < k_j
        deleteat!(deck.cards, k_j)
        deleteat!(deck.cards, k_i)
    else
        deleteat!(deck.cards, k_i)
        deleteat!(deck.cards, k_j)
    end
    players::Vector{Player} = []
    for p in 2:num_players
        push!(players, Player(p, Hole(p, Card(Ranks.NULL, Suits.NULL), Card(Ranks.NULL, Suits.NULL))))
    end
    board::Board = Board(Card(Ranks.NULL, Suits.NULL), Card(Ranks.NULL, Suits.NULL), Card(Ranks.NULL, Suits.NULL), Card(Ranks.NULL, Suits.NULL), Card(Ranks.NULL, Suits.NULL))
    predeal_shuffle!(deck)
    for t in 1:trials
        pot -= 1
        deal!(deck, players...)
        flop!(deck, board)
        turn!(deck, board)
        river!(deck, board)
        hands::Vector{Hand} = []
        push!(hands, score(player, board))
        for p in players
            p_i, p_j = 0, 0
            if compare(p.hole.card1, p.hole.card2) >= 0
                if p.hole.card1.suit == p.hole.card2.suit
                    p_i = 15-Int(p.hole.card1.rank)
                    p_j = 15-Int(p.hole.card2.rank)
                else
                    p_i = 15-Int(p.hole.card2.rank)
                    p_j = 15-Int(p.hole.card1.rank)
                end
            else
                if p.hole.card1.suit == p.hole.card2.suit
                    p_i = 15-Int(p.hole.card2.rank)
                    p_j = 15-Int(p.hole.card1.rank)
                else
                    p_i = 15-Int(p.hole.card1.rank)
                    p_j = 15-Int(p.hole.card2.rank)
                end
            end
            if player_range[p_i, p_j] == 0
                continue
            end
            push!(hands, score(p, board))
        end
        if length(hands) == 1
            pot += 1
            reshuffle!(deck)
            continue
        end
        best_hands::Vector{Hand} = compare_hands(hands...)
        for hand in best_hands
            if hand.playerID == 1
                pot += num_players / length(best_hands)
            end
        end
        reshuffle!(deck)
    end
    equity = pot / (num_players * trials)
    put!(chan, (i, j, equity))
end

@time main(ARGS)

max_stars_repo_name: ttowncompiled/EQPoker

max_issues_repo_name: ttowncompiled/EQPoker

max_forks_repo_name: ttowncompiled/EQPoker


text:
# Show for elements
Base.show(io::IO, h::HyperPlane) = print(io, "HyperPlane($(h.a), $(h.β))")
Base.show(io::IO, h::HalfSpace) = print(io, "HalfSpace($(h.a), $(h.β))")
Base.show(io::IO, l::Line) = print(io, "Line($(l.a))")
Base.show(io::IO, r::Ray) = print(io, "Ray($(r.a))")

Base.summary(it::Polyhedra.AbstractRepIterator{T, ElemT}) where {T, ElemT} = "$(length(it))-element iterator of $ElemT"

function _print_typeinfo_prefix(io, v)
    x = Base.typeinfo_prefix(io, v)
    if x isa Tuple
        # On Julia v1.5, `x` is `prefix, implicit`
        prefix = x[1]
    else
        prefix = x
    end
    print(io, prefix)
end

# Inspired from Base.show_vector
function show_repit(io::IO, v::Polyhedra.AbstractRepIterator, print_prefix::Bool, start_str::String, end_str::String, join_str::String=",")
    if print_prefix
        _print_typeinfo_prefix(io, v)
    end
    io = IOContext(io, :typeinfo => eltype(v), :compact => get(io, :compact, true))
    limited = get(io, :limit, false)
    limited = true
    if limited && length(v) > 20
        Base.show_delim_array(io, v, start_str, join_str, " \u2026 " * end_str, false, 1, 20)
    else
        Base.show_delim_array(io, v, start_str, join_str, end_str, false)
    end
end

# Inspired from Base.showarray with repr=false and Base.print_matrix
function show_repit(io::IO, it::Polyhedra.AbstractRepIterator, ::MIME"text/plain")
    print(io, summary(it))
    if !isempty(it)
        if !get(io, :limit, false)
            screenheight = screenwidth = typemax(Int)
        else
            sz = displaysize(io)
            # We use -4 as this is what is used in Base.print_matrix.
            # This probably accounts for the summary and the vertical dots
            screenheight, screenwidth = sz[1] - 4, sz[2]
        end
        if screenheight < length(it)
            Base.show_delim_array(io, it, ":\n ", "\n", "\n  \u22ee", false, 1, screenheight)
        else
            Base.show_delim_array(io, it, ":\n ", "\n", "", false)
        end
    end
end

Base.show(io::IO, v::Polyhedra.AbstractRepIterator) = show_repit(io, v, true, "[", "]")
Base.show(io::IO, mime::MIME"text/plain", it::Polyhedra.AbstractRepIterator) = show_repit(io, it, mime)

Base.summary(v::VRepresentation) = "V-representation $(typeof(v))"
Base.summary(h::HRepresentation) = "H-representation $(typeof(h))"
Base.summary(p::Polyhedron) = "Polyhedron $(typeof(p))"

show_reps(io::IO, args::Tuple, start_str::String, join_str::String, first::Bool) = first
function show_reps(io::IO, args::Tuple, start_str::String, join_str::String, first::Bool, rep, reps...)
    if !isempty(rep)
        if first
            print(io, start_str)
            first = false
        else
            print(io, join_str)
        end
        show_repit(io, rep, args...)
    end
    show_reps(io::IO, args, start_str, join_str, first, reps...)
end

show_vreps(io::IO, rep::HRepresentation, start_str::String, join_str::String, args...) = true
show_vreps(io::IO, rep::VRepresentation, start_str::String, join_str::String, args...) = show_reps(io, args, start_str, join_str, true, vreps(rep)...)
function show_vreps(io::IO, rep::Polyhedron, start_str::String, join_str::String, args...)
    if vrepiscomputed(rep)
        show_reps(io, args, start_str, join_str, true, vreps(rep)...)
    else
        true
    end
end
show_hreps(io::IO, rep::VRepresentation, start_str::String, join_str::String, args...) = true
show_hreps(io::IO, rep::HRepresentation, start_str::String, join_str::String, args...) = show_reps(io, args, start_str, join_str, true, hreps(rep)...)
function show_hreps(io::IO, rep::Polyhedron, start_str::String, join_str::String, args...)
    if hrepiscomputed(rep)
        show_reps(io, args, start_str, join_str, true, hreps(rep)...)
    else
        true
    end
end

_has_vrep(p::VRepresentation) = true
_has_vrep(p::HRepresentation) = false
_has_vrep(p::Polyhedron) = vrepiscomputed(p)

function Base.show(io::IO, rep::Rep)
    first = show_hreps(io, rep, "", " ∩ ", false, "", "", " ∩")
    if _has_vrep(rep) && !first
        print(io, " : ")
    end
    show_vreps(io, rep, "", " + ", false, "convexhull(", ")")
end
function Base.show(io::IO, mime::MIME"text/plain", rep::Rep)
    print(io, summary(rep))
    show_hreps(io, rep, ":\n", ",\n", mime)
    show_vreps(io, rep, ":\n", ",\n", mime)
end

#_length(hrep::HRepresentation) = nhyperplanes(hrep) + nhalfspaces(hrep)
#_length(vrep::VRepresentation) = nsympoints(vrep) + npoints(vrep) + nlines(vrep) + nrays(vrep)
#
#function _print(io::IO, it::ElemIt{<:HRepElement})
#    for h in it
#        print(io, " $(h.β)")
#        for j in eachindex(h.a)
#            print(io, " $(-h.a[j])")
#        end
#        println(io)
#    end
#end
#
#function _print(io::IO, it::ElemIt{<:VRepElement})
#    for v in it
#        print(io, " $(Int(ispoint(v)))")
#        c = coord(v)
#        for j = eachindex(c)
#            print(io, " $(c[j])")
#        end
#        println(io)
#    end
#end
#
#function Base.show(io::IO, rep::Representation{T}) where {T}
#    if typeof(rep) <: HRepresentation
#        print(io, "H")
#        ls = BitSet(1:nhyperplanes(rep))
#    else
#        print(io, "V")
#        ls = BitSet(1:nsympoints(rep)) ∪ BitSet(nsympoints(rep) + npoints(rep) + (1:nlines(rep)))
#    end
#    println(io, "-representation")
#
#    if !isempty(ls)
#        print(io, "linearity $(length(ls))");
#        for i in ls
#            print(io, " $i")
#        end
#        println(io)
#    end
#
#    println(io, "begin")
#    if T <: AbstractFloat
#        typename = "real"
#    elseif T <: Integer
#        typename = "integer"
#    else
#        typename = "rational"
#    end
#    println(io, " $(_length(rep)) $(N+1) $typename")
#    if typeof(rep) <: HRepresentation
#        _print.(io, hreps(rep))
#    else
#        _print.(io, vreps(rep))
#    end
#    print(io, "end")
#end

max_stars_repo_name: blegat/Polyhedra

max_issues_repo_name: blegat/Polyhedra

max_forks_repo_name: blegat/Polyhedra


text:
# Test SPN to MILP translation
import SPNetworks.MAP: spn2milp

@testset "SPN2MILP" begin
    # SPN with dichotomic sum nodes (<= 2 children)
    SPN = SumProductNetwork(
        [
            SumNode([3,2],[0.2,0.8]),               # 1
            SumNode([4,5],[0.625,0.375]),           # 2
            ProductNode([6,8]),                     # 3
            ProductNode([6,9]),                     # 4
            ProductNode([7,9]),                     # 5
            SumNode([10,11],[0.6,0.4]),             # 6
            SumNode([10,11],[0.1,0.9]),             # 7
            SumNode([12,13],[0.3,0.7]),             # 8
            SumNode([12,13],[0.8,0.2]),             # 9
            IndicatorFunction(1,1),                 # 10
            IndicatorFunction(1,2),                 # 11
            IndicatorFunction(2,1),                 # 12
            IndicatorFunction(2,2),                 # 13
        ]
    )

    spn2milp(stdou, SPN, [6,7,2,9,1,8])

end


text:
using CuArrays, CUDAdrv
using CuArrays.CURAND

using TraitSimulation, GLM, LinearAlgebra, Random, BenchmarkTools
using SnpArrays

function simulate_no_threads!(Y, trait::VCMTrait)
  fill!(Y, 0.0)
  TraitSimulation.VCM_trait_simulation(Y, trait.Z, trait.μ, trait.vc)
  return Y
end

"""
```
simulate(trait::VCMTrait, n::Integer)
```
This simulates a VCM trait n times under the desired variance component model, specified using the VCMTrait type.
"""
function simulate_no_threads(trait::VCMTrait, n::Integer)
  # pre-allocate output
  Y_n = [zeros(size(trait.μ)) for _ in 1:n] # replicates
  # do the simulation n times, each thread needs its own Z to collect randome effects
  for k in 1:n
      simulate_no_threads!(Y_n[k], trait)
  end
  return Y_n
end

n_sims = 1000

# @time test_cpu_no_threads = simulate_no_threads(trait, n_sims);

### cpu threaded

# @benchmark test_cpu_threads = simulate(trait, n_sims);

### GPU

struct VC_gpu{matT1, matT2, matT3, matT4}
	Σ::matT1 # n_traits by n_traits
	V::matT2 # n_people by n_people
	CholΣ::matT3 # cholesky decomposition of A
	CholV::matT4 # cholesk decomposition of B
	function VC_gpu(Σ, V) #inner constructor given A, B
        Σd = CuArray(Σ)
        Vd = CuArray(V)
        CholΣ = cholesky(Symmetric(Σd)).factors
        CholV = cholesky(Symmetric(Vd)).factors
        new{typeof(Σd), typeof(Vd), typeof(CholΣ), typeof(CholV)}(Σd, Vd, CholΣ, CholV) # stores these values (this is helpful so we don't have it inside the loop)
	end
end

function append_terms_gpu!(AB, summand)
	A_esc = esc(summand.args[2])	# elements in args are symbols,
	B_esc = esc(summand.args[3])
	push!(AB.args, :(VC_gpu($A_esc, $B_esc)))
end

"""
this vc macro allows us to create a vector of VarianceComponent objects for simulation so with_bigfloat_precis, precision::Integer)
so that the user can type out @vc V[1] ⊗ Σ[1] + V[2] ⊗ Σ[2] + .... + V[m] ⊗ Σ[m]
"""
macro vc_gpu(expression)
	n = length(expression.args)
	AB = :(VC_gpu[]) # AB is an empty vector of variance components list of symbols
	if expression.args[1] != :+ #if first argument is not plus (only one vc)
		summand = expression
		append_terms_gpu!(AB, summand)
	else #MULTIPLE VARIANCE COMPONENTS if the first argument is a plus (Sigma is a sum multiple variance components)
		for i in 2:n
			summand = expression.args[i]
			append_terms_gpu!(AB, summand)
		end
	end
	return(:($AB))
end

function simulate_matrix_normal_gpu!(Z_d, vc_gpu)
    randn!(Z_d)
    CuArrays.CUBLAS.trmm!('L', 'U', 'T', 'N', 1.0, vc_gpu.CholV, Z_d, similar(Z_d))
    CuArrays.CUBLAS.trmm!('R', 'U', 'N', 'N', 1.0, vc_gpu.CholΣ, Z_d, similar(Z_d))
    Z_d
end

"""
	VCM_trait_simulation(mu::Matrix{Float64}, vc::Vector{VarianceComponent})
For an evaluated mean matrix and vector of VarianceComponent objects, simulate from VCM.
"""
function VCM_trait_simulation_gpu(Y_d, μ_d, Z_d, vc_gpu) # for an evaluated matrix
	for i in eachindex(vc_gpu)
		simulate_matrix_normal_gpu!(Z_d, vc_gpu[i]) # this step aggregates the variance components by
		Y_d += Z_d # summing the independent matrix normals to Y, rewriting over Z for each variance component
	end
	Y_d += μ_d # add the mean back to shift the matrix normal
end

 function simulate_gpu!(Y_d, vc_gpu)
      μ_d = CuArray{Float64}(trait.μ)
      Z_d = CuArray{Float64}(trait.Z)
      Y_d .+= VCM_trait_simulation_gpu(Y_d, μ_d, Z_d, vc_gpu)
  end

function simulate_gpu(trait::VCMTrait, vc_gpu, n::Integer)
        # pre-allocate output
        Y_n = [zeros(CuArray{Float32}, size(trait.Z)) for _ in 1:n] # replicates
        # do the simulation n times, each thread needs its own Z to collect randome effects
        for k in 1:n
            simulate_gpu!(Y_n[k], vc_gpu)
        end
        return Y_n
    end

# @benchmark vc_gpu = @vc_gpu ΣA ⊗ 2(GRM + 0.01I_n) + ΣE ⊗ I_n;
# @benchmark test_gpu = simulate_gpu(trait, vc_gpu, n_sims);


text:
__precompile__()

module AlphaShapes

	using LinearAlgebraicRepresentation
	using TimerOutputs
	#using MATLAB
	#using Delaunay #BUG in package
	using SharedArrays
	using Distributed, Triangle
	using Combinatorics, DataStructures
	using Base.Threads
	const Lar = LinearAlgebraicRepresentation

	include("alpha_complex.jl")
	include("deWall.jl")
	include("geometry.jl")
	include("delaunayTriangulation.jl")
	const to = TimerOutput()

	export alphaFilter, alphaSimplex, delaunayTriangulation, to
end


text:
using Pkg # hideall
Pkg.activate("_literate/A-learning-networks/Project.toml")
Pkg.update();

# ## Preliminary steps
#
# Let's generate a `DataFrame` with some dummy regression data, let's also load the good old ridge regressor.

using MLJ, StableRNGs
import DataFrames
MLJ.color_off() # hide
Ridge = @load RidgeRegressor pkg=MultivariateStats

rng = StableRNG(551234) # for reproducibility

x1 = rand(rng, 300)
x2 = rand(rng, 300)
x3 = rand(rng, 300)
y = exp.(x1 - x2 -2x3 + 0.1*rand(rng, 300))

X = DataFrames.DataFrame(x1=x1, x2=x2, x3=x3)
first(X, 3) |> pretty

# Let's also prepare the train and test split which will be useful later on.

test, train = partition(eachindex(y), 0.8);

# ## Defining a learning network
#
# In MLJ, a *learning network* is a directed acyclic graph (DAG) whose *nodes* apply trained or untrained operations such as a `predict` or `transform` (trained) or `+`, `vcat` etc. (untrained).
# Learning networks can be seen as pipelines on steroids.
#
# Let's consider the following simple DAG:
#
# ![Operation DAG](/assets/diagrams/composite1.svg)
#
# It corresponds to a fairly standard regression workflow: the data is standardized, the target is transformed using a Box-Cox transformation, a ridge regression is applied and the result is converted back by inverting the transform.
#
# **Note**: actually  this DAG is simple enough that it could also have been done with a pipeline.
#
# ### Sources and nodes
#
# In MLJ a learning network starts at **source** nodes and flows through nodes (`X` and `y`) defining operations/transformations (`W`, `z`, `ẑ`, `ŷ`).
# To define the source nodes, use the `source` function, you should specify whether it's a target:

Xs = source(X)
ys = source(y)

# To define an "trained-operation" node, you must simply create a machine wrapping a model and another node (the data) and indicate which operation should be performed (e.g. `transform`):

stand = machine(Standardizer(), Xs)
W = transform(stand, Xs)

# You can `fit!` a trained-operation node at any point, MLJ will fit whatever it needs that is upstream of that node.
# In this case, there is just a source node upstream of `W` so fitting `W` will just fit the standardizer:

fit!(W, rows=train);

# If you want to get the transformed data, you can then call the node speciying on which part of the data the operation should be performed:

W()             # transforms all data
W(rows=test, )  # transforms only test data
W(X[3:4, :])    # transforms specific data

# Let's now define the other nodes:

box_model = UnivariateBoxCoxTransformer()
box = machine(box_model, ys)
z = transform(box, ys)

ridge_model = Ridge(lambda=0.1)
ridge = machine(ridge_model, W, z)
ẑ = predict(ridge, W)

ŷ = inverse_transform(box, ẑ)

# Note that we have not yet done any training, but if we now call `fit!` on `ŷ`, it will fit all nodes upstream of `ŷ` that need to be re-trained:

fit!(ŷ, rows=train);

# Now that `ŷ` has been fitted, you can apply the full graph on test data (or any compatible data). For instance, let's get the `rms` between the ground truth and the predicted values:

rms(y[test], ŷ(rows=test))

# ### Modifying hyperparameters
#
# Hyperparameters can be accessed using the dot syntax as usual.
# Let's modify the regularisation parameter of the ridge regression:

ridge_model.lambda = 5.0;

# Since the node `ẑ` corresponds to a machine that wraps `ridge_model`, that node has effectively changed and will be retrained:

fit!(ŷ, rows=train)
rms(y[test], ŷ(rows=test))

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using LoopVectorization, Test#, OffsetArrays

@testset "Array Wrappers" begin
  @show @__LINE__

  function addone!(y, x)
    @turbo for i ∈ eachindex(x, y)
      y[i] = x[i] + 1
    end
    y
  end

  x5 = PermutedDimsArray(rand(10, 10), (2, 1))
  y5 = PermutedDimsArray(rand(10, 10), (2, 1))
  y5t = x5 .+ 1
  @test LoopVectorization.check_args(x5)
  @test y5t == addone!(y5, x5)

  x6 = view(x5, 1:3, 1:3)
  y6 = view(y5, 1:3, 1:3)
  fill!(y6, NaN)
  @test LoopVectorization.check_args(x6)
  addone!(y6, x6) # also
  @test y5t == y5 # Test that `NaN`s replaced with correct answers, and that other values were untouched.

  A = rand(12, 13, 14, 15)
  pA = PermutedDimsArray(A, (3, 1, 4, 2))
  @test addone!(similar(pA), pA) == pA .+ 1
  ppA = PermutedDimsArray(pA, (4, 2, 3, 1))
  @test LoopVectorization.check_args(ppA)
  @test addone!(similar(ppA), ppA) == ppA .+ 1

  x = rand(10, 10)
  xc = copy(x)
  xv = view(x, 1:2:10, 1:3)
  xcv = view(xc, 1:2:10, 1:3)
  @test LoopVectorization.check_args(xv)
  addone!(xv, xv)
  xcv .+= 1
  @test x == xc
end


text:
"""
    download(gl::GitLink; verbose, tout = 60.0)
    download(upfun::Function, gl::GitLink; verbose, tout = 60.0)

Try a to sync the local repo with te remote.
Optionally, and after the pull, it will call `upfun(wdir)` to access the working directory.
This method will update the `pull` event token.
This method will sleep till (or timeout `tout`) the GitLink lock is free (which must by must of the time, but...).
Returns `true` if the action was successful.
"""
function download(upfun::Function, gl::GitLink; verbose = false, tout = 60.0)

    pull_ok = false

    lock(gl; tout) do

        # HARD PULL (Loop)
        pull_ok = hard_pull(gl; verbose, clearwd = true)
        if !pull_ok # Handle fail
            add_loop_frec!(gl, _LOOP_FREC_FAIL_PENALTY)
            return false
        end
        _set_pull_token(gl) # Aknowlage successful pull

        # call fun
        wdir = repo_dir(gl)
        mkpath(wdir)
        upfun(wdir)

    end

    return pull_ok
end
download(gl::GitLink; kwargs...) = download((wdir) -> nothing, gl; kwargs...)


text:
# test effectiveness of different τ
using Distributed;
addprocs(20);
@everywhere include("loadMod.jl");
@everywhere const GUROBI_ENV = Gurobi.Env();

N = 5;
dataList = Dict();

caseList = [13,33,123];
Δt = 0.25;
iterMax = 20;
T = 24;
NN = 1000;

τList = [2,4,6,8,10];

for ci in 1:length(caseList)
    dataList[ci] = Dict();

    for τ in τList
        for j in procs()
            remotecall_fetch(readInData,j,ci,caseList,T,τ);
            #remotecall_fetch(readInData_old,j,T,ωSet0,10000,0);
        end

        # train the stochastic strategy
        startT = time();
        cutDictPG = preGen(T, Δt, N, iterMax);
        # cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist = solveMain(τ, T, Δt, N, true, false, 20, 20, cutDictPG);
        cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist = solveMain(T, Δt, N, false, false,
            max(Int64(round(500/N)),20), max(Int64(round(500/N)),20),cutDictPG);
        elapsedT = time() - startT;

        pathListData = pmap(i -> simuPath(T,pDistr), 1:NN);
        pathDict = Dict();
        for i in 1:NN
            pathDict[i] = pathListData[i];
        end
        solSDDP, LBSDDP, costSDDP = exeForward(T, Δt, NN, false, pathDict);
        listSDDP = [costSDDP[i] for i in 1:NN];
        meanSDDP = mean(listSDDP);
        sigmaSDDP = std(listSDDP);
        dataList[ci][τ] = [LBHist,LBSDDP,listSDDP,meanSDDP,sigmaSDDP,timeHist,elapsedT];
        save("tauResults_$(ci).jld","data",dataList);
    end
end


text:
using Base.Test
import Base: tuple_type_cons, tuple_type_head, tuple_type_tail, in, ==, isless, convert,
             length, eltype, start, next, done, show
using WeakRefStrings

fastmap(f, xs...) = map(f, xs...)
@generated function fastmap(f, xs::NTuple{N}...) where N
    args = [:(xs[$j][i])  for j in 1:nfields(xs)]
    :(Base.@ntuple $N i -> f($(args...)))
end

export @NT

eltypes(::Type{Tuple{}}) = Tuple{}
eltypes(::Type{T}) where {T<:Tuple} =
    tuple_type_cons(eltype(tuple_type_head(T)), eltypes(tuple_type_tail(T)))
eltypes(::Type{T}) where {T<:NamedTuple} = map_params(eltype, T)
eltypes(::Type{T}) where T <: Pair = map_params(eltypes, T)
eltypes(::Type{T}) where T<:AbstractArray{S, N} where {S, N} = S
Base.@pure astuple(::Type{T}) where {T<:NamedTuple} = Tuple{T.parameters...}
astuple(::Type{T}) where {T<:Tuple} = T

# sizehint, making sure to return first argument
_sizehint!(a::Array{T,1}, n::Integer) where {T} = (sizehint!(a, n); a)
_sizehint!(a::AbstractArray, sz::Integer) = a

# argument selectors
left(x, y) = x
right(x, y) = y

# tuple and NamedTuple utilities

@generated function ith_all(i, xs::Tuple)
    :(Base.@ntuple $(nfields(xs)) j -> xs[j][i])
end

@generated function ith_all(i, n::NamedTuple)
    Expr(:block,
         :(@Base._inline_meta),
         Expr(:tuple, [ Expr(:ref, Expr(:., :n, Expr(:quote, fieldname(n,f))), :i) for f = 1:nfields(n) ]...))
end

@generated function foreach(f, x::Union{NamedTuple, Tuple}, xs::Union{NamedTuple, Tuple}...)
    args = [:(getfield(getfield(xs, $j), i))  for j in 1:nfields(xs)]
    :(Base.@nexprs $(nfields(x)) i -> f(getfield(x, i), $(args...)); nothing)
end

@inline foreach(f, a::Pair) = (f(a.first); f(a.second))
@inline foreach(f, a::Pair, b::Pair) = (f(a.first, b.first); f(a.second, b.second))

fieldindex(x, i::Integer) = i
fieldindex(x::NamedTuple, s::Symbol) = findfirst(x->x===s, fieldnames(x))

astuple(t::Tuple) = t

@generated function astuple(n::NamedTuple)
    Expr(:tuple, [ Expr(:., :n, Expr(:quote, fieldname(n,f))) for f = 1:nfields(n) ]...)
end

# lexicographic order product iterator

import Base: length, eltype, start, next, done

abstract type AbstractProdIterator end

struct Prod2{I1, I2} <: AbstractProdIterator
    a::I1
    b::I2
end

product(a) = a
product(a, b) = Prod2(a, b)
eltype(::Type{Prod2{I1,I2}}) where {I1,I2} = Tuple{eltype(I1), eltype(I2)}
length(p::AbstractProdIterator) = length(p.a)*length(p.b)

function start(p::AbstractProdIterator)
    s1, s2 = start(p.a), start(p.b)
    s1, s2, (done(p.a,s1) || done(p.b,s2))
end

function prod_next(p, st)
    s1, s2 = st[1], st[2]
    v2, s2 = next(p.b, s2)
    doneflag = false
    if done(p.b, s2)
        v1, s1 = next(p.a, s1)
        if !done(p.a, s1)
            s2 = start(p.b)
        else
            doneflag = true
        end
    else
        v1, _ = next(p.a, s1)
    end
    return (v1,v2), (s1,s2,doneflag)
end

next(p::Prod2, st) = prod_next(p, st)
done(p::AbstractProdIterator, st) = st[3]

struct Prod{I1, I2<:AbstractProdIterator} <: AbstractProdIterator
    a::I1
    b::I2
end

product(a, b, c...) = Prod(a, product(b, c...))
eltype(::Type{Prod{I1,I2}}) where {I1,I2} = tuple_type_cons(eltype(I1), eltype(I2))

function next(p::Prod{I1,I2}, st) where {I1,I2}
    x = prod_next(p, st)
    ((x[1][1],x[1][2]...), x[2])
end

# sortperm with counting sort

sortperm_fast(x) = sortperm(x)
sortperm_fast(x::StringVector) = sortperm(convert(StringVector{WeakRefString{UInt8}}, x))

function sortperm_fast(v::Vector{T}) where T<:Integer
    n = length(v)
    if n > 1
        min, max = extrema(v)
        rangelen = max - min + 1
        if rangelen < div(n,2)
            return sortperm_int_range(v, rangelen, min)
        end
    end
    return sortperm(v, alg=MergeSort)
end

function sortperm_int_range(x::Vector{T}, rangelen, minval) where T<:Integer
    offs = 1 - minval
    n = length(x)

    where = fill(0, rangelen+1)
    where[1] = 1
    @inbounds for i = 1:n
        where[x[i] + offs + 1] += 1
    end
    cumsum!(where, where)

    P = Vector{Int}(n)
    @inbounds for i = 1:n
        label = x[i] + offs
        wl = where[label]
        P[wl] = i
        where[label] = wl+1
    end

    return P
end

# sort the values in v[i0:i1] in place, by array `by`
Base.@noinline function sort_sub_by!(v, i0, i1, by, order, temp)
    empty!(temp)
    sort!(v, i0, i1, MergeSort, order, temp)
end

Base.@noinline function sort_sub_by!(v, i0, i1, by::PooledArray, order, temp)
    empty!(temp)
    sort!(v, i0, i1, MergeSort, order, temp)
end

Base.@noinline function sort_sub_by!(v, i0, i1, by::Vector{T}, order, temp) where T<:Integer
    min = max = by[v[i0]]
    @inbounds for i = i0+1:i1
        val = by[v[i]]
        if val < min
            min = val
        elseif val > max
            max = val
        end
    end
    rangelen = max-min+1
    n = i1-i0+1
    if rangelen <= n
        sort_int_range_sub_by!(v, i0-1, n, by, rangelen, min, temp)
    else
        empty!(temp)
        sort!(v, i0, i1, MergeSort, order, temp)
    end
    v
end

# in-place counting sort of x[ioffs+1:ioffs+n] by values in `by`
function sort_int_range_sub_by!(x, ioffs, n, by, rangelen, minval, temp)
    offs = 1 - minval

    where = fill(0, rangelen+1)
    where[1] = 1
    @inbounds for i = 1:n
        where[by[x[i+ioffs]] + offs + 1] += 1
    end
    cumsum!(where, where)

    length(temp) < n && resize!(temp, n)
    @inbounds for i = 1:n
        xi = x[i+ioffs]
        label = by[xi] + offs
        wl = where[label]
        temp[wl] = xi
        where[label] = wl+1
    end

    @inbounds for i = 1:n
        x[i+ioffs] = temp[i]
    end
    x
end

function append_n!(X, val, n)
    l = length(X)
    resize!(X, l+n)
    for i in (1:n)+l
        @inbounds X[i] = val
    end
    X
end

const _namedtuple_cache = Dict{Tuple{Vararg{Symbol}}, Type}()
function namedtuple(fields...)
    if haskey(_namedtuple_cache, fields)
        return _namedtuple_cache[fields]
    else
        NT = eval(:(@NT($(fields...))))
        _namedtuple_cache[fields] = NT
        return NT
    end
end

"""
`arrayof(T)`

Returns the type of `Columns` or `Vector` suitable to store
values of type T. Nested tuples beget nested Columns.
"""
Base.@pure function arrayof(S)
    T = strip_unionall(S)
    if T == Union{}
        Vector{Union{}}
    elseif T<:Tuple
        Columns{T, Tuple{map(arrayof, T.parameters)...}}
    elseif T<:NamedTuple
        if nfields(T) == 0
            Columns{namedtuple(), namedtuple()}
        else
            Columns{T,namedtuple(fieldnames(T)...){map(arrayof, T.parameters)...}}
        end
    elseif T<:DataValue
        DataValueArray{T.parameters[1],1}
    elseif T<:Union{String,WeakRefString}
        StringArray{T, 1}
    elseif T<:Pair
        Columns{T, Pair{map(arrayof, T.parameters)...}}
    else
        Vector{T}
    end
end

@inline strip_unionall_params(T::UnionAll) = strip_unionall_params(T.body)
@inline strip_unionall_params(T) = map(strip_unionall, T.parameters)

Base.@pure function promote_union(T::Type)
    if isa(T, Union)
        return promote_type(T.a, promote_union(T.b))
    else
        return T
    end
end

Base.@pure function strip_unionall(T)
    if isleaftype(T) || T == Union{}
        return T
    elseif isa(T, TypeVar)
        T.lb === Union{} && return strip_unionall(T.ub)
        return Any
    elseif T == Tuple
        return Any
    elseif T<:Tuple
        if any(x->x <: Vararg, T.parameters)
            # we only keep known-length tuples
            return Any
        else
            return Tuple{strip_unionall_params(T)...}
        end
    elseif T<:NamedTuple
        if isa(T, Union)
            return promote_union(T)
        else
            NT = namedtuple(fieldnames(T)...)
            return NT{strip_unionall_params(T)...}
        end
    elseif isa(T, UnionAll)
        return Any
    elseif isa(T, Union)
        return promote_union(T)
    elseif T.abstract
        return T
    else
        return Any
    end
end

Base.@pure function _promote_op(f, ::Type{S}) where S
    t = Core.Inference.return_type(f, Tuple{Base._default_type(S)})
    strip_unionall(t)
end

Base.@pure function _promote_op(f, ::Type{S}, ::Type{T}) where {S,T}
    t = Core.Inference.return_type(f, Tuple{Base._default_type(S),
                                        Base._default_type(T)})
    strip_unionall(t)
end

@inline _map(f, p::Pair) = f(p.first) => f(p.second)
@inline _map(f, args...) = map(f, args...)

# The following is not inferable, this is OK because the only place we use
# this doesn't need it.

function _map_params(f, T, S)
    (f(_tuple_type_head(T), _tuple_type_head(S)),
     _map_params(f, _tuple_type_tail(T), _tuple_type_tail(S))...)
end

_map_params(f, T::Type{Tuple{}},S::Type{Tuple{}}) = ()

map_params(f, ::Type{T}, ::Type{S}) where {T,S} = f(T,S)
map_params(f, ::Type{T}) where {T} = map_params((x,y)->f(x), T, T)
map_params(f, ::Type{T}) where T <: Pair{S1, S2} where {S1, S2} = Pair{f(S1), f(S2)}
@inline _tuple_type_head(::Type{T}) where {T<:Tuple} = Base.tuple_type_head(T)
@inline _tuple_type_tail(::Type{T}) where {T<:Tuple} = Base.tuple_type_tail(T)

#function map_params{N}(f, T::Type{T} where T<:Tuple{Vararg{Any,N}}, S::Type{S} where S<: Tuple{Vararg{Any,N}})
Base.@pure function map_params(f, ::Type{T}, ::Type{S}) where {T<:Tuple,S<:Tuple}
    if nfields(T) != nfields(S)
        MethodError(map_params, (typeof(f), T,S))
    end
    Tuple{_map_params(f, T,S)...}
end

_tuple_type_head(T::Type{NT}) where {NT<: NamedTuple} = fieldtype(NT, 1)

Base.@pure function _tuple_type_tail(T::Type{NT}) where NT<: NamedTuple
    Tuple{Base.argtail(NT.parameters...)...}
end

Base.@pure @generated function map_params(f, ::Type{T}, ::Type{S}) where {T<:NamedTuple,S<:NamedTuple}
    if fieldnames(T) != fieldnames(S)
        MethodError(map_params, (T,S))
    end
    if nfields(T) == 0 && nfields(S) == 0
        return T
    end
    NT = Expr(:macrocall, :(NamedTuples.$(Symbol("@NT"))), fieldnames(T)...)
    :($NT{_map_params(f, T, S)...})
end

@inline function concat_tup(a::NamedTuple, b::NamedTuple)
    concat_tup_type(typeof(a), typeof(b))(a..., b...)
end
@inline concat_tup(a::Tup, b::Tup) = (a..., b...)
@inline concat_tup(a::Tup, b) = (a..., b)
@inline concat_tup(a, b::Tup) = (a, b...)
@inline concat_tup(a, b) = (a..., b...)

Base.@pure function concat_tup_type(T::Type{<:Tuple}, S::Type{<:Tuple})
    Tuple{T.parameters..., S.parameters...}
end

Base.@pure function concat_tup_type(::Type{T}, ::Type{S}) where {
           T<:NamedTuple,S<:NamedTuple}
    nfields(T) == 0 && nfields(S) == 0 ?
        namedtuple() :
        namedtuple(fieldnames(T)...,
                   fieldnames(S)...){T.parameters...,
                                     S.parameters...}
end

Base.@pure function concat_tup_type(T::Type, S::Type)
    Tuple{T,S}
end

# check to see if array has shared data
# used in flush! to create a copy of the arrays
function isshared(x)
    try
        resize!(x, length(x))
        false
    catch err
        if err isa ErrorException && err.msg == "cannot resize array with shared data"
            return true
        else
            rethrow(err)
        end
    end
end

compact_mem(x) = x
compact_mem(x::StringArray{String}) = convert(StringArray{WeakRefString{UInt8}}, x)


text:
"""
    Compute unit weights.

    Inputs:
        gp: A fitted Gaussian process
        sentinels: a 2×R matrix of border sentinels
        weights: an R-vector of border weights

    Output:
        w_unit: A vector of weights for each unit in `gp`
"""
function weight_at_units(gp::GPE, sentinels::Matrix, weights::Vector)
    @assert size(sentinels, 2) == length(weights)
    K_Xb = cov(gp.kernel, gp.x, sentinels)
    K_XX = gp.cK
    w_unit = (K_XX \ K_Xb) * weights ./ sum(weights)
    return w_unit
end

"""
    Compute unit weights for the LATE estimator with uniform border weights.
"""
function weights_unif(gpT::GPE, gpC::GPE, nb::Int, border)
    Xb = sentinels(border, nb)
    w_border = ones(nb)
    w_units = [weight_at_units(gpT, Xb, w_border); -weight_at_units(gpC, Xb, w_border)]
    return Xb, w_border, w_units
end

"""
    Compute unit weights for the LATE estimator with inverse-variance border weights.
"""
function weights_invvar(gpT::GPE, gpC::GPE, nb::Int, border)
    Xb = sentinels(border, nb)
    _, Σpost = cliff_face(gpT, gpC, Xb)
    w_border = Σpost \ ones(nb)
    w_units = [weight_at_units(gpT, Xb, w_border); -weight_at_units(gpC, Xb, w_border)]
    return Xb, w_border, w_units
end

"""
    Compute unit weights for the projected finite-population LATE estimator.
"""
function weights_proj(gpT::GPE, gpC::GPE, border, maxdist::Float64)
    Xb_treat = projection_points(gpT, border, maxdist)
    Xb_ctrol = projection_points(gpC, border, maxdist)
    Xb_projected = [Xb_treat Xb_ctrol]
    Xb = Xb_projected
    nb = size(Xb, 2)
    w_border = ones(nb)
    w_units = [weight_at_units(gpT, Xb, w_border); -weight_at_units(gpC, Xb, w_border)]
    return Xb, w_border, w_units
end

"""
    Compute unit weights for the projected land LATE estimator.
"""
function weights_geo(gpT::GPE, gpC::GPE, border, region, maxdist::Float64, gridspace::Float64)
    Xb_projected, w_border = infinite_proj_sentinels(
            gpT, gpC, border, region, maxdist, gridspace)
    Xb = Xb_projected
    w_units = [weight_at_units(gpT, Xb, w_border); -weight_at_units(gpC, Xb, w_border)]
    return Xb, w_border, w_units
end

"""
    Compute unit weights for the projected superpopulation LATE estimator.
"""
function weights_pop(gpT::GPE, gpC::GPE, density::Function, border, region, maxdist::Float64, gridspace::Float64)
    Xb_projected, w_border = infinite_proj_sentinels(
            gpT, gpC, border, region, maxdist, gridspace; density=density)
    Xb = Xb_projected
    w_units = [weight_at_units(gpT, Xb, w_border); -weight_at_units(gpC, Xb, w_border)]
    return Xb, w_border, w_units
end

"""
    Compute unit weights for the population density-weighted LATE estimator.
"""
function weights_rho(gpT::GPE, gpC::GPE, density::Function, nb::Int, border)
    Xb = sentinels(border, nb)

    w_border = density.(Xb[1,:], Xb[2,:])
    w_units = [weight_at_units(gpT, Xb, w_border); -weight_at_units(gpC, Xb, w_border)]
    return Xb, w_border, w_units
end


text:
using QuantumTensors
using Base.Test

function test_to_sparse()
    tol = 1e-10
    N = 5
    @test norm(((destroy(N) |> sparse) - spdiagm(sqrt(1:(N-1)),1,N,N))|>full )< tol
    @test norm(((create(N) |> sparse) - spdiagm(sqrt(1:(N-1)),-1,N,N))|>full )< tol
    @test norm(((jz(N) |> sparse) - spdiagm(-N:N,0,2N+1,2N+1))|>full) < tol

    jzm = jz(N) |> sparse
    jp = jplus(N) |> sparse
    jm1 = jp'
    jm2 = jminus(N) |> sparse
    @test norm((jm1-jm2)|> full) < tol

    comm_rel = jp*jm1 - jm1*jp
    @test norm((comm_rel-2jzm)|>full) < tol
end

function test_to_full()
    tol = 1e-10
    N = 5
    @test norm((destroy(N) |> full) - diagm(sqrt(1:(N-1)),1))< tol
    @test norm((create(N) |> full) - diagm(sqrt(1:(N-1)),-1))< tol
    @test norm((sigma(2,1,2) |> full) - [0 1; 0 0])< tol
    @test norm((sigma(2,1,1) |> full) - [1 0; 0 0])< tol
    @test norm((sigma(2,2,1) |> full) - [0 0; 1 0])< tol
    @test norm((sigma(2,2,2) |> full) - [0 0; 0 1])< tol
end

function test_transpose()
    tol = 1e-10
    @test norm((sigma(2,1,2) |>full) - (sigma(2,2,1) |> transpose |>full)) < tol
    @test norm((destroy(5) |>full) - (create(5) |> transpose |>full)) < tol
end

function test_convert_promote_ops()
  tol = 1e-10
  N = 5
  # a = destroy(N)
  # ac = convert(DiagonalOp{Complex128}, a)
  idt = IdentityOp{N}()
  idtf = convert(DiagonalOp{Float64}, idt)
  idtz1 = convert(DiagonalOp{Complex128}, idt)
  idtz2 = convert(DiagonalOp{Complex128}, idtf)
  idtz3, idtz4 = promote(idtf, idtz1)
  @test idtz4 === idtz1
  @test isa(idtf, DiagonalOp{Float64})
  @test isa(idtz1, DiagonalOp{Complex128})
  @test isa(idtz2, DiagonalOp{Complex128})
  @test isa(idtz3, DiagonalOp{Complex128})
  @test idtf.jkm == (0,0,0)
  @test idtz1.jkm == (0,0,0)
  @test idtz2.jkm == (0,0,0)
  @test idtz3.jkm == (0,0,0)
  @test norm(idtf.vals - ones(N)) < tol
  @test norm(idtz1.vals - ones(N)) < tol
  @test norm(idtz2.vals - ones(N)) < tol
  @test norm(idtz3.vals - ones(N)) < tol

  @test all_close(idt, idtf)
  @test all_close(idt, idtz1)
  @test all_close(idt, idtz2)
  @test all_close(idt, idtz3)
  @test all_close(idt, idtz4)
  @test all_close(idtf, idtz1)
  @test all_close(idtf, idtz2)
  @test all_close(idtf, idtz3)
  @test all_close(idtf, idtz4)
end

function test_convert_promote_tops()
  tol = 1e-10
  N = 5
  a = destroy(N)
  bc = convert(DiagonalOp{Complex128}, destroy(N))

  ac, bc = promote(a, bc)
  @test bc === bc
  @test all_close(ac, a)
  abc = tensor(a, bc)
  acbc = tensor(ac, bc)

  @test typeof(abc) === typeof(acbc)
  @test all_close(abc, acbc)
end

function test_convert_promote_topsum()
  tol = 1e-10
  N = 5
  a = destroy(N)
  bc = convert(DiagonalOp{Complex128}, destroy(N))
  ii = IdentityOp{N}()

  ac, bc = promote(a, bc)
  @test bc === bc
  @test all_close(ac, a)
  abc = tensor(a, bc)
  acbc = tensor(ac, bc)

  iibc = tensor(ii, bc)
  aii = tensor(a, ii)

  tops = TensorOpProductSum(abc, iibc, aii)
end

function test_mul()
    tol = 1e-10
    a = destroy(5)
    @test norm((a*a |> full) -(a|> full)*(a|> full)) < tol
    @test norm((a*a' |> full) -(a|> full)*(a|> full)')<tol
end

function test_apply_serial(N=20,tol=1e-8)
  println("Testing serial application of operators")
  a = destroy(N)
  b = destroy(N)
  c = destroy(N)
  II = IdentityOp{N}()
  tops = (tensor(a', b, II), tensor(a, b', II), tensor(II, II, c'*c))
  topsum = TensorOpProductSum(tops...)

  B = randn(N,N,N)
  C = zeros(N,N,N)
  coeffs = [1.,1.,1.]
  apply_serial!(coeffs, topsum, B, C)
  # Profile.clear()
  println("apply_serial!: TensorOpProductSum ")
  gc()
  @time for k=1:10
      apply_serial!(coeffs, topsum, B,C)
  end

  topskarr = collect([TensorOpProductKernel(top) for top in tops]);
  topsarr = collect(tops);

  println("apply_serial!: [TensorOpProductKernel]")
  Cpp = zeros(C)
  apply_serial!(coeffs, topskarr, B, Cpp)
  gc()
  @time for k=1:10
      apply_serial!(coeffs, topskarr, B, Cpp)
  end

  @test (norm((Cpp-C)[:])) < tol

  println("apply_serial!: [TensorOpProduct]")
  Cpp = zeros(C)
  apply_serial!(coeffs, topsarr, B, Cpp)
  gc()
  @time for k=1:10
      apply_serial!(coeffs, topsarr, B, Cpp)
  end
  @test norm((Cpp-C)[:]) < tol

  Cp = zeros(C)
  Cpf = sub(Cp,1:length(C));
  Bf = sub(B,1:length(B));
  top_sps = [sparse(top) for top = tops]
  combinedsp = sum(top_sps)
  println("combined CSR matrix")
  A_mul_B!(Cpf, combinedsp, Bf)
  gc()
  @time for k=1:10
      A_mul_B!(Cpf, combinedsp, Bf)
  end
  @test norm((Cp-C)[:]) < tol

  println("apply_serial!: [CSR-mat]")
  apply_serial!(coeffs, top_sps, Bf, Cpf)
  gc()
  @time for k=1:10
      apply_serial!(coeffs, top_sps, Bf, Cpf)
  end
  @test norm((Cp-C)[:]) < tol
end

module SymbolicTest
  const sqrt_arr1 = sqrt(1:1000000)
  const sqrt_arr2 = sqrt(1:1000000)
  function a1ad2_symbolic!(y, x, α=0, β=1)
    N1, N2=size(x)

    if α != 1
      if α == 0
        fill!(y, 0)
      else
        scale!(y, α)
      end
    end

    for l=1:(N2-1)
      @simd for k=1:(N1-1)
        @inbounds y[k,l+1] +=  β * sqrt_arr1[k]*sqrt_arr2[l]*x[k+1,l]
      end
    end
    y
  end
  function a1ad2_symbolic2!(y, x, α=0, β=1)
    N1, N2=size(x)

    if α != 1
      if α == 0
        fill!(y, 0)
      else
        scale!(y, α)
      end
    end

    offset1 = -1
    offset2 = 1

    oIDX = 1
    iIDX = 1 - stride(x,1)*offset1 - stride(x,2)*offset2
    for l=1:N2
      lflag = 2 <= l <=N2
      @simd for k=1:N1
        kflag = 1 <= k <= N1-1
        if kflag && lflag
          @inbounds y[oIDX] +=  β * sqrt_arr1[k]*sqrt_arr2[l-1]*x[iIDX]
        end
        oIDX +=1
        iIDX +=1
      end
    end
    y
  end
  function a1ad2_symbolic3!(y, x, α=0, β=1)
    N1, N2=size(x)

    if α != 1
      if α == 0
        fill!(y, 0)
      else
        scale!(y, α)
      end
    end

    offset1 = -1
    offset2 = 1

    oIDX = 1
    iIDX = 1 - stride(x,1)*offset1 - stride(x,2)*offset2
    for l=1:N2
      lflag = 2 <= l <=N2
      for k=1:N1
        kflag = 1 <= k <= N1-1
        if kflag && lflag
          @inbounds y[oIDX] +=  β * sqrt_arr1[k]*sqrt_arr2[l-1]*x[iIDX]
        end
        oIDX +=1
        iIDX +=1
      end
    end
    y
  end
end

function test_A_mul_B(N1=10,N2=10)
    println("Testing application of single tensor product operators to vectors")
    tol = 1e-5
    a1 = destroy(N1)
    ad2 = create(N2)
    a1ad2 = tensor(a1,ad2)
    p1 = randn(N1)
    p2 = randn(N2)
    x = tensor(p1, p2)

    println("Tensor prod with allocation:")
    y1 = zeros(N1,N2)
    A_mul_B!(y1, a1ad2, x)
    gc()
    @time begin
        y1 = A_mul_B!(zeros(N1,N2), a1ad2, x)
    end

    println("Tensor prod (inplace):")
    y5 = zeros(N1,N2)
    gc()
    @time begin
        A_mul_B!(y5, a1ad2, x)
    end
    @test norm(y1[:] - y5[:]) < tol

    println("Direct usage of TensorOpProductKernel:")
    y10 = zeros(N1,N2)
    A_mul_B!(2., a1ad2.kernel, x, 3., y10)
    y10 = zeros(N1,N2)
    gc()
    @time begin
        A_mul_B!(2., a1ad2.kernel, x, 3., y10)
    end
    @test norm(2*y1[:] - y10[:]) < tol

    println("'Symbolic' application of specialized function:")
    y11 = zeros(N1,N2)
    SymbolicTest.a1ad2_symbolic!(y11, x, 3., 2.)
    y11 *= 0
    gc()
    @time begin
        SymbolicTest.a1ad2_symbolic!(y11, x, 3., 2.)
    end
    @test norm(2*y1[:] - y11[:]) < tol

    println("'Symbolic' application of second specialized function:")
    y12 = zeros(N1,N2)
    SymbolicTest.a1ad2_symbolic2!(y12, x, 3., 2.)
    y12 *= 0
    gc()
    @time begin
        SymbolicTest.a1ad2_symbolic2!(y12, x, 3., 2.)
    end
    @test norm(2*y1[:] - y12[:]) < tol

    println("'Symbolic' application of third specialized function:")
    y12 = zeros(N1,N2)
    SymbolicTest.a1ad2_symbolic3!(y12, x, 3., 2.)
    y12 *= 0
    gc()
    @time begin
        SymbolicTest.a1ad2_symbolic3!(y12, x, 3., 2.)
    end
    @test norm(2*y1[:] - y12[:]) < tol

    println("Tensorproduct of individual products:")
    y2 = tensor(((a1|>sparse)*p1),((ad2|>sparse)*p2))
    gc()
    @time y2 = tensor(((a1|>sparse)*p1),((ad2|>sparse)*p2))
    @test norm(y1[:] - y2[:]) < tol

    println("CSR mat-vec:")

    xv = x[:]

    a1ad2sp = (a1ad2 |> sparse)
    y3 = a1ad2sp * xv
    gc()
    @time y3 = a1ad2sp * xv
    @test norm(y1[:] - y3[:]) < tol

    println("CSR mat-vec (inplace):")
    A_mul_B!(2., a1ad2sp, xv, 1., y3)
    y3[:] = y1[:]
    gc()
    @time A_mul_B!(2., a1ad2sp, xv, 1., y3)
    @test norm(3*y1[:] - y3[:]) < tol

    println("Size of vectors: ", Base.summarysize(x)/1024^2, " MB")
    println("Size of sparse matrix rep: ", Base.summarysize(a1ad2sp)/1024^2, " MB")
    println("Size of tensor matrix rep: ", Base.summarysize(a1ad2)/1024^2, " MB")
end

# addprocs(3)
versioninfo()
test_to_sparse()
test_to_full()
test_transpose()
test_mul()
test_convert_promote_ops()
test_convert_promote_tops()
test_convert_promote_topsum()
test_apply_serial(20)
test_A_mul_B(100,100)


text:
# Defaults for many of the interface functions are defined here as fallback.
# Methods here should take a trait instance as first argument and should already be defined
# in the `interface.jl` first as a generic f(geom) method.

## Coords
# Four options in SF, xy, xyz, xym, xyzm
const default_coord_names = (:X, :Y, :Z, :M)

coordnames(t::AbstractGeometryTrait, geom) = default_coord_names[1:ncoord(t, geom)]

# Maybe hardcode dimension order? At least for X and Y?
x(t::AbstractPointTrait, geom) = getcoord(t, geom, findfirst(isequal(:X), coordnames(geom)))
y(t::AbstractPointTrait, geom) = getcoord(t, geom, findfirst(isequal(:Y), coordnames(geom)))
z(t::AbstractPointTrait, geom) = getcoord(t, geom, findfirst(isequal(:Z), coordnames(geom)))
m(t::AbstractPointTrait, geom) = getcoord(t, geom, findfirst(isequal(:M), coordnames(geom)))

is3d(::AbstractGeometryTrait, geom) = :Z in coordnames(geom)
ismeasured(::AbstractGeometryTrait, geom) = :M in coordnames(geom)
isempty(::AbstractGeometryTrait, geom) = false

## Points
ngeom(::AbstractPointTrait, geom) = 0
getgeom(::AbstractPointTrait, geom) = nothing
getgeom(::AbstractPointTrait, geom, i) = nothing

## LineStrings
npoint(t::AbstractCurveTrait, geom) = ngeom(t, geom)
getpoint(t::AbstractCurveTrait, geom) = getgeom(t, geom)
getpoint(t::AbstractCurveTrait, geom, i) = getgeom(t, geom, i)
startpoint(t::AbstractCurveTrait, geom) = getpoint(t, geom, 1)
endpoint(t::AbstractCurveTrait, geom) = getpoint(t, geom, npoint(geom))

## Polygons
nring(t::AbstractPolygonTrait, geom) = ngeom(t, geom)
getring(t::AbstractPolygonTrait, geom) = getgeom(t, geom)
getring(t::AbstractPolygonTrait, geom, i) = getgeom(t, geom, i)
getexterior(t::AbstractPolygonTrait, geom) = getring(t, geom, 1)
nhole(t::AbstractPolygonTrait, geom) = nring(t, geom) - 1
gethole(t::AbstractPolygonTrait, geom) = (getgeom(t, geom, i) for i in 2:ngeom(t, geom))
gethole(t::AbstractPolygonTrait, geom, i) = getring(t, geom, i + 1)
npoint(t::AbstractPolygonTrait, geom) = sum(npoint(p) for p in getring(t, geom))
getpoint(t::AbstractPolygonTrait, geom) = flatten((p for p in getpoint(r)) for r in getring(t, geom))

## MultiPoint
npoint(t::AbstractMultiPointTrait, geom) = ngeom(t, geom)
getpoint(t::AbstractMultiPointTrait, geom) = getgeom(t, geom)
getpoint(t::AbstractMultiPointTrait, geom, i) = getgeom(t, geom, i)

## MultiLineString
nlinestring(t::AbstractMultiCurveTrait, geom) = ngeom(t, geom)
getlinestring(t::AbstractMultiCurveTrait, geom) = getgeom(t, geom)
getlinestring(t::AbstractMultiCurveTrait, geom, i) = getgeom(t, geom, i)
npoint(t::AbstractMultiCurveTrait, geom) = sum(npoint(ls) for ls in getgeom(t, geom))
getpoint(t::AbstractMultiCurveTrait, geom) = flatten((p for p in getpoint(ls)) for ls in getgeom(t, geom))

## MultiPolygon
npolygon(t::AbstractMultiPolygonTrait, geom) = ngeom(t, geom)
getpolygon(t::AbstractMultiPolygonTrait, geom) = getgeom(t, geom)
getpolygon(t::AbstractMultiPolygonTrait, geom, i) = getgeom(t, geom, i)
nring(t::AbstractMultiPolygonTrait, geom) = sum(nring(p) for p in getpolygon(t, geom))
getring(t::AbstractMultiPolygonTrait, geom) = flatten((r for r in getring(p)) for p in getpolygon(t, geom))
npoint(t::AbstractMultiPolygonTrait, geom) = sum(npoint(r) for r in getring(t, geom))
getpoint(t::AbstractMultiPolygonTrait, geom) = flatten((p for p in getpoint(r)) for r in getring(t, geom))

## Surface
npatch(t::AbstractPolyhedralSurfaceTrait, geom)::Integer = ngeom(t, geom)
getpatch(t::AbstractPolyhedralSurfaceTrait, geom) = getgeom(t, geom)
getpatch(t::AbstractPolyhedralSurfaceTrait, geom, i::Integer) = getgeom(t, geom, i)

## Default iterator
getgeom(t::AbstractGeometryTrait, geom) = (getgeom(t, geom, i) for i in 1:ngeom(t, geom))
getcoord(t::AbstractPointTrait, geom) = (getcoord(t, geom, i) for i in 1:ncoord(t, geom))

## Special geometries
npoint(::LineTrait, geom) = 2
npoint(::TriangleTrait, geom) = 3
nring(::TriangleTrait, geom) = 1
npoint(::RectangleTrait, geom) = 4
nring(::RectangleTrait, geom) = 1
npoint(::QuadTrait, geom) = 4
nring(::QuadTrait, geom) = 1
npoint(::PentagonTrait, geom) = 5
nring(::PentagonTrait, geom) = 1
npoint(::HexagonTrait, geom) = 6
nring(::HexagonTrait, geom) = 1

# TODO Only simple if it's also not intersecting itself, except for its endpoints
issimple(t::AbstractCurveTrait, geom) = allunique((getpoint(t, geom, i) for i in 1:npoint(geom)-1)) && allunique((getpoint(t, geom, i) for i in 2:npoint(t, geom)))
isclosed(t::AbstractCurveTrait, geom) = getpoint(t, geom, 1) == getpoint(t, geom, npoint(t, geom))
isring(t::AbstractCurveTrait, geom) = issimple(t, geom) && isclosed(t, geom)

issimple(t::AbstractMultiPointTrait, geom) = allunique((getgeom(t, geom)))

issimple(t::AbstractMultiCurveTrait, geom) = all(issimple.(getgeom(t, geom)))
isclosed(t::AbstractMultiCurveTrait, geom) = all(isclosed.(getgeom(t, geom)))

crs(::AbstractGeometryTrait, geom) = nothing
extent(::AbstractGeometryTrait, geom) = nothing

# FeatureCollection
getfeature(t::AbstractFeatureCollectionTrait, fc) = (getfeature(t, fc, i) for i in 1:nfeature(t, fc))

# Backwards compatibility
coordinates(t::AbstractPointTrait, geom) = collect(getcoord(t, geom))
coordinates(t::AbstractGeometryTrait, geom) = collect(coordinates.(getgeom(t, geom)))
function coordinates(t::AbstractFeatureTrait, feature)
    geom = geometry(feature)
    isnothing(geom) ? [] : coordinates(geom)
end
coordinates(t::AbstractFeatureCollectionTrait, fc) = [coordinates(f) for f in getfeature(fc)]

Base.convert(T::Type, ::AbstractGeometryTrait, geom) = error("Conversion is enabled for type $T, but not implemented. Please report this issue to the package maintainer.")

# Subtraits

"""
    subtrait(t::AbstractGeometryTrait)

Gets the expected, possible abstract, (sub)trait for subgeometries (retrieved with
[`getgeom`](@ref)) of trait `t`. This follows the [Type hierarchy](@ref) of Simple Features.

# Examples
```jldoctest; setup = :(using GeoInterface)
julia> GeoInterface.subtrait(LineStringTrait())
AbstractPointTrait
julia> GeoInterface.subtrait(PolygonTrait())  # Any of LineStringTrait, LineTrait, LinearRingTrait
AbstractLineStringTrait
```
```jldoctest; setup = :(using GeoInterface)
# `nothing` is returned when there's no subtrait or when it's not known beforehand
julia> isnothing(GeoInterface.subtrait(PointTrait()))
true
julia> isnothing(GeoInterface.subtrait(GeometryCollectionTrait()))
true
```
"""
subtrait(::AbstractPointTrait) = nothing
subtrait(::AbstractCurveTrait) = AbstractPointTrait
subtrait(::AbstractCurvePolygonTrait) = AbstractCurveTrait
subtrait(::AbstractPolygonTrait) = AbstractLineStringTrait
subtrait(::AbstractPolyhedralSurfaceTrait) = AbstractPolygonTrait
subtrait(::TINTrait) = TriangleTrait
subtrait(::AbstractMultiPointTrait) = AbstractPointTrait
subtrait(::AbstractMultiLineStringTrait) = AbstractLineStringTrait
subtrait(::AbstractMultiPolygonTrait) = AbstractPolygonTrait
subtrait(::AbstractMultiSurfaceTrait) = AbstractSurfaceTrait
subtrait(::AbstractGeometryCollectionTrait) = nothing
subtrait(::AbstractGeometryTrait) = nothing  # fallback


text:
function test_membership(members)
    bs = ByteSet(members)
    refset = Set{UInt8}([UInt8(i) for i in  members])
    @test refset == Set{UInt8}(collect(bs))
    @test all(i -> in(i, bs), refset)
end

function test_inversion(bs)
    inv = ~bs
    all = true
    for i in 0x00:0xff
        all &= (in(i, bs) ⊻ in(i, inv))
    end
    @test all
end

@testset "Instantiation" begin
    @test isempty(ByteSet())
    @test iszero(length(ByteSet()))

    for set in ["hello", "kdjy82zxxcbnpw", [0x00, 0x1a, 0xff, 0xf8, 0xd2]]
        test_membership(set)
    end
end

@testset "Min/max" begin
    @test_throws ArgumentError maximum(ByteSet())
    @test_throws ArgumentError minimum(ByteSet())
    @test minimum(ByteSet("xylophone")) == UInt8('e')
    @test maximum(ByteSet([0xa1, 0x0f, 0x4e, 0xf1, 0x40, 0x39])) == 0xf1
end

@testset "Contiguity" begin
    @test ScanByte.is_contiguous(ByteSet(0x03:0x41))
    @test ScanByte.is_contiguous(ByteSet())
    @test ScanByte.is_contiguous(ByteSet(0x51))
    @test ScanByte.is_contiguous(ByteSet(0xc1:0xd2))
    @test ScanByte.is_contiguous(ByteSet(0x00:0xff))

    @test !ScanByte.is_contiguous(ByteSet([0x12:0x3a; 0x3c:0x4a]))
    @test !ScanByte.is_contiguous(ByteSet([0x01, 0x02, 0x04, 0x05]))
end

@testset "Inversion" begin
    test_inversion(ByteSet())
    test_inversion(ByteSet(0x00:0xff))
    test_inversion(ByteSet([0x04, 0x06, 0x91, 0x92]))
    test_inversion(ByteSet(0x54:0x71))
    test_inversion(ByteSet(0x12:0x11))
    test_inversion(ByteSet("abracadabra"))
end


text:
# SOCIAL PLANNER PROBLEM

function solveSocialPlanner(p)

  # first solve baseline
  eqInit = [ 2.00, 0.55, 0.06,
             0.72, 0.83, 1.73]

  eq,res = solveBGP(p,eqInit)

  # initial values for the solver and minimization
  socPolInit = [eq.qmin; eq.x]
  eqSocInit = [eq.qbar;eq.cactiv;eq.xout]

  cevMax = [calCev(eq.g, eq.cactivtot, eq.g, eq.cactivtot,p)]
  @printf("SOCIAL PLANNER PROBLEM\n")
  @printf("-----------------------\n")
  @printf("Baseline Welfare is %3.4f (xˡ = %3.4f, xʰ = %3.4f, qminˡ = %3.4f, qminʰ = %3.4f) \n", cevMax[1],eq.x[1],eq.x[2],eq.qmin[1],eq.qmin[2])

  # optimization
  # ------------------------------------------------------------------------------------------------
  opt = Opt(:LN_NELDERMEAD, length(socPolInit))
  ftol_rel!(opt,1e-8)
  min_objective!(opt, (x, grad) -> socialPlannerFunc(x, p, eqSocInit, cevMax, eq.g, eq.cactivtot)[1])
  minf,minpol,ret = NLopt.optimize(opt, socPolInit)
  nEvals = opt.numevals
  #-------------------------------------------------------------------------------------------------

  if ret == :FTOL_REACHED || ret == :XTOL_REACHED
  	@printf("Optimal policy found!\n\n")
        eqSoc = socialPlannerFunc(minpol,  p, eqSocInit, cevMax, eq.g, eq.cactivtot)[2]
	res = Dict(:cev=>-minf, :ret=>ret, :opt=>opt, :nEvals=>nEvals)
  	return eqSoc, res
  else
	@printf("Optimal policy could not solved!\n\n")
    eqSoc = EqObj()
	return eqSoc, ret
  end
end

######################################################################################################################

function socialPlannerFunc(socPol, p, eqSocInit, cevMax, gBase, cactivtotBase)
  if (any(socPol.<=0.0))
    cev = -Inf
    eq  = EqObj()
  else

    eq = EqObj()   # initialize equilibrium object

    # gausslegendre for integration
    eq.qMinAll = 1.0e-8
    eq.qMaxAll = 10.0
    eq.node,eq.weight = gausslegendre(1000)
    eq.node = ((eq.qMaxAll - eq.qMinAll)/2)*eq.node .+ (eq.qMaxAll + eq.qMinAll)/2    # adjust the limits

    # social planner policy variables
    eq.qmin = socPol[1:2]
    eq.x    = socPol[3:4]

    res = nlsolve((eqnd,x)->socEqFunc!(eqnd,x,eq,p), eqSocInit, method = :trust_region,inplace = true,show_trace=false,iterations=50)

    if !res.f_converged
        #println("Couldn't solve!!!")
        cev = -Inf
    else
        socEqFunc!(similar(res.zero),res.zero,eq,p)   # evaluate at minimizer
        eqSocInit .= res.zero
        # welfare calculation
        cev = calCev(gBase, cactivtotBase, eq.g, eq.cactivtot, p)
    end

    if cev>cevMax[1]
      @printf("Welfare improved to %3.4f (xˡ = %3.4f, xʰ = %3.4f, qminˡ = %3.4f, qminʰ = %3.4f) \n",
               cev,eq.x[1],eq.x[2],eq.qmin[1],eq.qmin[2])
      cevMax[1] = cev
    end
  end

  return -cev, eq

end

######################################################################################################################

function socEqFunc!(eqnd,eqv,eq,p)
  # reject out of bounds

  if (any(eqv.<=0.0))
    eqnd = -1000.0*ones(size(eqv)[1])

  else
    # load in guesses for solution
    eq.qbar   = eqv[1]
    eq.cactiv = eqv[2:3]
    eq.xout   = eqv[4]

    # flow of free products
    eq.cactivtot = sum(eq.cactiv)
    eq.cinac     = 1.0 - eq.cactivtot
    eq.ws        = 0.0                   # just to make things work

    # aggregate creative destruction rates
    eq.τs = eq.cactiv.*eq.x .+ p.𝛂*eq.xout
    eq.τ  = sum(eq.τs)
    # find updated eqvars
    labordem!(eq,p)
    qualityDist!(eq,p)
    qbarActive!(eq,p)

    # equilibrium differences
    eqnd[1]   = eq.qbarAct     - (p.ε/(p.ε-1.0))^(p.ε - 1)
    eqnd[2:3].= eq.cactiv     .- eqv[2:3]       # product shares
    eqnd[4]   = eq.skilled_lab - p.Lˢ           # skilled labor market clearing

  end
end


text:
module MPFI

export
    Interval,
    bisect,
    blow,
    diam,
    diam_abs,
    diam_rel,
    left,
    mag,
    mid,
    mig,
    isbounded,
    right

import
    Base: precision, string, print, show, showcompact, promote_rule,
        promote, convert, +, *, -, /, exp, isinf, isnan, nan, inf, sqrt,
        square, exp, exp2, expm1, cosh, sinh, tanh, sech, csch, coth, inv,
        sqrt, cbrt, abs, log, log2, log10, log1p, sin, cos, tan, sec,
        csc, acos, asin, atan, acosh, asinh, atanh, isempty, union,
        intersect, in, cmp, ldexp, rand, rand!

type Interval <: Number
    left_prec::Clong
    left_sign::Cint
    left_exp::Clong
    left_d::Ptr{Void}
    right_prec::Clong
    right_sign::Cint
    right_exp::Clong
    right_d::Ptr{Void}
    function Interval()
        N = get_bigfloat_precision()
        z = new(zero(Clong), zero(Cint), zero(Clong), C_NULL,
                zero(Clong), zero(Cint), zero(Clong), C_NULL)
        ccall((:mpfi_init2,:libmpfi), Void, (Ptr{Interval}, Clong), &z, N)
        finalizer(z, MPFI_clear)
        return z
    end
end

MPFI_clear(mpfi::Interval) = ccall((:mpfi_clear, :libmpfi), Void, (Ptr{Interval},), &mpfi)

# Include useful macros and utility functions
include("utils.jl")

Interval(x::Interval) = x

for (fJ, fC) in ((:si,:Clong), (:ui,:Culong), (:d,:Float64))
    @eval begin
        function Interval(x::($fC))
            z = Interval()
            ccall(($(string(:mpfi_set_,fJ)), :libmpfi), Int32, (Ptr{Interval}, ($fC)), &z, x)
            return z
        end
    end
end

function Interval(x::BigInt)
    z = Interval()
    ccall((:mpfi_set_z, :libmpfi), Int32, (Ptr{Interval}, Ptr{BigInt}), &z, &x)
    return z
end

function Interval(x::BigFloat)
    z = Interval()
    ccall((:mpfi_set_fr, :libmpfi), Int32, (Ptr{Interval}, Ptr{BigFloat}), &z, &x)
    return z
end

function Interval(x::String, base::Int)
    z = Interval()
    err = ccall((:mpfi_set_str, :libmpfi), Int32, (Ptr{Interval}, Ptr{Uint8}, Int32), &z, x, base)
    if err != 0; error("Invalid input"); end
    return z
end
Interval(x::String) = Interval(x, 10)

Interval(x::Integer) = Interval(BigInt(x))

Interval(x::Union(Bool,Int8,Int16,Int32)) = Interval(convert(Clong,x))
Interval(x::Union(Uint8,Uint16,Uint32)) = Interval(convert(Culong,x))

Interval(x::Float32) = Interval(float64(x))
Interval(x::Rational) = Interval(num(x)) / Interval(den(x))

# Dyadic constructors
Interval(x::Real, y::Real) = Interval(promote(x, y)...)
for (fJ, fC) in ((:si,:Clong), (:ui,:Culong), (:d,:Float64))
    @eval begin
        function Interval(x::($fC), y::($fC))
            z = Interval()
            ccall(($(string(:mpfi_interv_,fJ)), :libmpfi), Int32, (Ptr{Interval}, ($fC), ($fC)), &z, x, y)
            return z
        end
    end
end

function Interval(x::BigInt, y::BigInt)
    z = Interval()
    ccall((:mpfi_interv_z, :libmpfi), Int32, (Ptr{Interval}, Ptr{BigInt}, Ptr{BigInt}), &z, &x, &y)
    return z
end

function Interval(x::BigFloat, y::BigFloat)
    z = Interval()
    ccall((:mpfi_interv_fr, :libmpfi), Int32, (Ptr{Interval}, Ptr{BigFloat}, Ptr{BigFloat}), &z, &x, &y)
    return z
end

Interval(x::MathConst) = convert(Interval, x)

# Promotion of constants
promote_rule{s}(::Type{MathConst{s}}, ::Type{Interval}) = Interval

# Conversions to Interval
convert(::Type{Interval}, x::Rational) = Interval(x)
convert(::Type{Interval}, x::Real) = Interval(x)
convert(::Type{Interval}, x::MathConst) = Interval(big(x))
function convert(::Type{Interval}, ::MathConst{:π})
    z = Interval()
    ccall((:mpfi_const_pi,:libmpfi), Cint, (Ptr{Interval},), &z)
    return z
end
function convert(::Type{Interval}, ::MathConst{:γ})
    z = Interval()
    ccall((:mpfi_const_euler,:libmpfi), Cint, (Ptr{Interval},), &z)
    return z
end
function convert(::Type{Interval}, ::MathConst{:catalan})
    z = Interval()
    ccall((:mpfi_const_catalan,:libmpfi), Cint, (Ptr{Interval},), &z)
    return z
end

# Conversions from Interval
function convert(::Type{BigFloat}, x::Interval)
    z = BigFloat()
    ccall((:mpfi_get_fr,:libmpfi), Void, (Ptr{BigFloat}, Ptr{Interval}), &z, &x)
    return z
end
convert(::Type{Float64}, x::Interval) =
    ccall((:mpfi_get_d,:libmpfi), Float64, (Ptr{Interval},), &x)

for to in (Int8, Int16, Int32, Int64, Uint8, Uint16, Uint32, Uint64, BigInt, Float32)
    @eval begin
        function convert(::Type{$to}, x::Interval)
            convert($to, convert(BigFloat, x))
        end
    end
end
convert(::Type{Integer}, x::Interval) = convert(BigInt, x)
convert(::Type{FloatingPoint}, x::Interval) = convert(BigFloat, x)

# Interval functions with floating-point results
for f in (:diam_abs, :diam_rel, :diam, :mag, :mig, :mid)
    @eval function $(f)(x::Interval)
        z = BigFloat()
        ccall(($(string(:mpfi_,f)), :libmpfi), Int32, (Ptr{BigFloat}, Ptr{Interval}), &z, &x)
        return z
    end
end

# Basic operations between intervals
for (fJ, fC) in ((:+,:add), (:-,:sub), (:*,:mul), (:/,:div))
    @eval begin
        function ($fJ)(x::Interval, y::Interval)
            z = Interval()
            ccall(($(string(:mpfi_,fC)),:libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Ptr{Interval}), &z, &x, &y)
            return z
        end
    end
end

# More efficient commutative operations
for (fJ, fC) in ((:+, :add), (:*, :mul))
    @eval begin
        function ($fJ)(a::Interval, b::Interval, c::Interval)
            z = Interval()
            ccall(($(string(:mpfi_,fC)), :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Ptr{Interval}), &z, &a, &b)
            ccall(($(string(:mpfi_,fC)), :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Ptr{Interval}), &z, &z, &c)
            return z
        end
        function ($fJ)(a::Interval, b::Interval, c::Interval, d::Interval)
            z = Interval()
            ccall(($(string(:mpfi_,fC)), :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Ptr{Interval}), &z, &a, &b)
            ccall(($(string(:mpfi_,fC)), :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Ptr{Interval}), &z, &z, &c)
            ccall(($(string(:mpfi_,fC)), :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Ptr{Interval}), &z, &z, &d)
            return z
        end
        function ($fJ)(a::Interval, b::Interval, c::Interval, d::Interval, e::Interval)
            z = Interval()
            ccall(($(string(:mpfi_,fC)), :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Ptr{Interval}), &z, &a, &b)
            ccall(($(string(:mpfi_,fC)), :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Ptr{Interval}), &z, &z, &c)
            ccall(($(string(:mpfi_,fC)), :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Ptr{Interval}), &z, &z, &d)
            ccall(($(string(:mpfi_,fC)), :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Ptr{Interval}), &z, &z, &e)
            return z
        end
    end
end

# Basic arithmetic without promotion
# Unsigned addition
function +(x::Interval, c::Culong)
    z = Interval()
    ccall((:mpfi_add_ui, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Culong), &z, &x, c)
    return z
end
+(c::Culong, x::Interval) = x + c
+(c::Unsigned, x::Interval) = x + convert(Culong, c)
+(x::Interval, c::Unsigned) = x + convert(Culong, c)

# Signed addition
function +(x::Interval, c::Clong)
    z = Interval()
    ccall((:mpfi_add_si, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Clong), &z, &x, c)
    return z
end
+(c::Clong, x::Interval) = x + c
+(x::Interval, c::Signed) = x + convert(Clong, c)
+(c::Signed, x::Interval) = x + convert(Clong, c)

# Float64 addition
function +(x::Interval, c::Float64)
    z = Interval()
    ccall((:mpfi_add_d, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Float64), &z, &x, c)
    return z
end
+(c::Float64, x::Interval) = x + c
+(c::Float32, x::Interval) = x + convert(Float64, c)
+(x::Interval, c::Float32) = x + convert(Float64, c)

# BigInt addition
function +(x::Interval, c::BigInt)
    z = Interval()
    ccall((:mpfi_add_z, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Ptr{BigInt}), &z, &x, &c)
    return z
end
+(c::BigInt, x::Interval) = x + c

# BigFloat addition
function +(x::Interval, c::BigFloat)
    z = Interval()
    ccall((:mpfi_add_fr, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Ptr{BigFloat}), &z, &x, &c)
    return z
end
+(c::BigInt, x::Interval) = x + c

# Unsigned subtraction
function -(x::Interval, c::Culong)
    z = Interval()
    ccall((:mpfi_sub_ui, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Culong), &z, &x, c)
    return z
end
function -(c::Culong, x::Interval)
    z = Interval()
    ccall((:mpfi_ui_sub, :libmpfi), Int32, (Ptr{Interval}, Culong, Ptr{Interval}), &z, c, &x)
    return z
end
-(x::Interval, c::Unsigned) = -(x, convert(Culong, c))
-(c::Unsigned, x::Interval) = -(convert(Culong, c), x)

# Signed subtraction
function -(x::Interval, c::Clong)
    z = Interval()
    ccall((:mpfi_sub_si, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Clong), &z, &x, c)
    return z
end
function -(c::Clong, x::Interval)
    z = Interval()
    ccall((:mpfi_si_sub, :libmpfi), Int32, (Ptr{Interval}, Clong, Ptr{Interval}), &z, c, &x)
    return z
end
-(x::Interval, c::Signed) = -(x, convert(Clong, c))
-(c::Signed, x::Interval) = -(convert(Clong, c), x)

# Float64 subtraction
function -(x::Interval, c::Float64)
    z = Interval()
    ccall((:mpfi_sub_d, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Float64), &z, &x, c)
    return z
end
function -(c::Float64, x::Interval)
    z = Interval()
    ccall((:mpfi_d_sub, :libmpfi), Int32, (Ptr{Interval}, Float64, Ptr{Interval}), &z, c, &x)
    return z
end
-(x::Interval, c::Float32) = -(x, convert(Float64, c))
-(c::Float32, x::Interval) = -(convert(Float64, c), x)

# BigInt subtraction
function -(x::Interval, c::BigInt)
    z = Interval()
    ccall((:mpfi_sub_z, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Ptr{BigInt}), &z, &x, &c)
    return z
end
function -(c::BigInt, x::Interval)
    z = Interval()
    ccall((:mpfi_z_sub, :libmpfi), Int32, (Ptr{Interval}, Ptr{BigInt}, Ptr{Interval}), &z, &c, &x)
    return z
end

# BigFloat subtraction
function -(x::Interval, c::BigFloat)
    z = Interval()
    ccall((:mpfi_sub_fr, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Ptr{BigFloat}), &z, &x, &c)
    return z
end
function -(c::BigFloat, x::Interval)
    z = Interval()
    ccall((:mpfi_fr_sub, :libmpfi), Int32, (Ptr{Interval}, Ptr{BigFloat}, Ptr{Interval}), &z, &c, &x)
    return z
end

# Unsigned multiplication
function *(x::Interval, c::Culong)
    z = Interval()
    ccall((:mpfi_mul_ui, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Culong), &z, &x, c)
    return z
end
*(c::Culong, x::Interval) = x * c
*(c::Unsigned, x::Interval) = x * convert(Culong, c)
*(x::Interval, c::Unsigned) = x * convert(Culong, c)

# Signed multiplication
function *(x::Interval, c::Clong)
    z = Interval()
    ccall((:mpfi_mul_si, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Clong), &z, &x, c)
    return z
end
*(c::Clong, x::Interval) = x * c
*(c::Signed, x::Interval) = x * convert(Clong, c)
*(x::Interval, c::Signed) = x * convert(Clong, c)

# Float64 multiplication
function *(x::Interval, c::Float64)
    z = Interval()
    ccall((:mpfi_mul_d, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Float64), &z, &x, c)
    return z
end
*(c::Float64, x::Interval) = x * c
*(c::Float32, x::Interval) = x * convert(Float64, c)
*(x::Interval, c::Float32) = x * convert(Float64, c)

# BigInt multiplication
*(c::Signed, x::Interval) = x * convert(Clong, c)
function *(x::Interval, c::BigInt)
    z = Interval()
    ccall((:mpfi_mul_z, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Ptr{BigInt}), &z, &x, &c)
    return z
end
*(c::BigInt, x::Interval) = x * c

# BigFloat multiplication
function *(x::Interval, c::BigFloat)
    z = Interval()
    ccall((:mpfi_mul_fr, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Ptr{BigFloat}), &z, &x, &c)
    return z
end
*(c::BigFloat, x::Interval) = x * c

# Unsigned division
function /(x::Interval, c::Culong)
    z = Interval()
    ccall((:mpfi_div_ui, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Culong), &z, &x, c)
    return z
end
function /(c::Culong, x::Interval)
    z = Interval()
    ccall((:mpfi_ui_div, :libmpfi), Int32, (Ptr{Interval}, Culong, Ptr{Interval}), &z, c, &x)
    return z
end
/(x::Interval, c::Unsigned) = /(x, convert(Culong, c))
/(c::Unsigned, x::Interval) = /(convert(Culong, c), x)

# Signed division
function /(x::Interval, c::Clong)
    z = Interval()
    ccall((:mpfi_div_si, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Clong), &z, &x, c)
    return z
end
function /(c::Clong, x::Interval)
    z = Interval()
    ccall((:mpfi_si_div, :libmpfi), Int32, (Ptr{Interval}, Clong, Ptr{Interval}), &z, c, &x)
    return z
end
/(x::Interval, c::Signed) = /(x, convert(Clong, c))
/(c::Signed, x::Interval) = /(convert(Clong, c), x)

# Float64 division
function /(x::Interval, c::Float64)
    z = Interval()
    ccall((:mpfi_div_d, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Float64), &z, &x, c)
    return z
end
function /(c::Float64, x::Interval)
    z = Interval()
    ccall((:mpfi_d_div, :libmpfi), Int32, (Ptr{Interval}, Float64, Ptr{Interval}), &z, c, &x)
    return z
end
/(x::Interval, c::Float32) = /(x, convert(Float64, c))
/(c::Float32, x::Interval) = /(convert(Float64, c), x)

# BigInt division
function /(x::Interval, c::BigInt)
    z = Interval()
    ccall((:mpfi_div_z, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Ptr{BigInt}), &z, &x, &c)
    return z
end
function /(c::BigInt, x::Interval)
    z = Interval()
    ccall((:mpfi_z_div, :libmpfi), Int32, (Ptr{Interval}, Ptr{BigInt}, Ptr{Interval}), &z, &c, &x)
    return z
end

# BigFloat division
function /(x::Interval, c::BigFloat)
    z = Interval()
    ccall((:mpfi_div_fr, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Ptr{BigFloat}), &z, &x, &c)
    return z
end
function /(c::BigFloat, x::Interval)
    z = Interval()
    ccall((:mpfi_fr_div, :libmpfi), Int32, (Ptr{Interval}, Ptr{BigFloat}, Ptr{Interval}), &z, &c, &x)
    return z
end

function precision(x::Interval)
    return ccall((:mpfi_get_prec, :libmpfi), Clong, (Ptr{Interval},), &x)
end

function -(x::Interval)
    z = Interval()
    ccall((:mpfi_neg, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}), &z, &x)
    return z
end

function square(x::Interval)
    z = Interval()
    ccall((:mpfi_sqr, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}), &z, &x)
    return z
end

function sqrt(x::Interval)
    z = Interval()
    ccall((:mpfi_sqrt, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}), &z, &x)
    if isnan(z)
        throw(DomainError())
    end
    return z
end

for f in (:exp,:exp2,:expm1,:cosh,:sinh,:tanh,:sech,:csch,:coth,:inv,
          :sqrt,:cbrt,:abs,:log,:log2,:log10,:log1p,:sin,:cos,:tan,:sec,
          :csc,:acos,:asin,:atan,:acosh,:asinh,:atanh)
    @eval function $f(x::Interval)
        z = Interval()
        ccall(($(string(:mpfi_,f)), :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}), &z, &x)
        return z
    end
end

# Set-like functions
function in(x::Interval, int::Interval)
    return ccall((:mpfi_is_inside, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}), &x, &int) > 0
end

for (fJ, fC) in ((:si,:Clong), (:ui,:Culong), (:d,:Float64))
    @eval begin
        function in(x::($fC), int::Interval)
            return ccall(($(string(:mpfi_is_inside_,fJ)), :libmpfi), Int32, (($fC), Ptr{Interval}), x, &int) > 0
        end
    end
end

for (fJ, fC) in ((:z,:BigInt), (:fr,:BigFloat))
    @eval begin
        function in(x::($fC), int::Interval)
            return ccall(($(string(:mpfi_is_inside_,fJ)), :libmpfi), Int32, (Ptr{($fC)}, Ptr{Interval}), &x, &int) > 0
        end
    end
end

function intersect(x::Interval, y::Interval)
    z = Interval()
    ccall((:mpfi_intersect, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Ptr{Interval}), &z, &x, &y)
    return z
end

function union(x::Interval, y::Interval)
    z = Interval()
    ccall((:mpfi_union, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Ptr{Interval}), &z, &x, &y)
    return z
end

function bisect(x::Interval)
    z1, z2 = Interval(), Interval()
    ccall((:mpfi_bisect, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Ptr{Interval}), &z1, &z2, &x)
    return z1, z2
end

function blow(x::Interval, y::Float64)
    z = Interval()
    ccall((:mpfi_blow, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Float64), &z, &x, y)
    return z
end

function ldexp(x::Interval, y::Culong)
    z = Interval()
    ccall((:mpfi_mul_2exp, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Culong), &z, &x, y)
    return z
end
function ldexp(x::Interval, y::Clong)
    z = Interval()
    ccall((:mpfi_mul_2si, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}, Clong), &z, &x, y)
    return z
end
ldexp(x::Interval, y::Signed) = ldexp(x, convert(Clong, y))
ldexp(x::Interval, y::Unsigned) = ldexp(x, convert(Culong, y))

function left(x::Interval)
    z = BigFloat()
    ccall((:mpfi_get_left, :libmpfi), Int32, (Ptr{BigFloat}, Ptr{Interval}), &z, &x)
    return z
end

function right(x::Interval)
    z = BigFloat()
    ccall((:mpfi_get_right, :libmpfi), Int32, (Ptr{BigFloat}, Ptr{Interval}), &z, &x)
    return z
end

function isbounded(x::Interval)
    return ccall((:mpfi_bounded_p, :libmpfi), Int32, (Ptr{Interval},), &x) != 0
end

function isempty(x::Interval)
    return ccall((:mpfi_empty_p, :libmpfi), Int32, (Ptr{Interval},), &x) != 0
end

function isnan(x::Interval)
    return ccall((:mpfi_nan_p, :libmpfi), Int32, (Ptr{Interval},), &x) != 0
end

function isinf(x::Interval)
    return ccall((:mpfi_inf_p, :libmpfi), Int32, (Ptr{Interval},), &x) != 0
end

function cmp2(x::Interval, y::Interval)
    return ccall((:mpfi_cmp, :libmpfi), Int32, (Ptr{Interval}, Ptr{Interval}), &x, &y)
end

@have_bigrng begin
function rand!(x::Interval, A::Array{BigFloat})
    for i = 1:length(A)
        A[i] = rand(x)
    end
    A
end

function rand(r::BigRNG, x::Interval)
    z = BigFloat()
    ccall((:mpfi_urandom, :libmpfi), Void,
          (Ptr{BigFloat}, Ptr{Interval}, Ptr{BigRNG}),
           &z, &x, &r)
    return z
end
rand(x::Interval) = rand(Base.Random.DEFAULT_BIGRNG, x)
rand(x::Interval, dims::Dims) = rand!(x, Array(BigFloat, dims))
rand(x::Interval, dims::Int...) = rand!(x, Array(BigFloat, dims...))

end # @have_bigrng

function string(x::Interval)
    # We use the alternate constructor to avoid the GC,
    # since MPFI manages the memory here
    left = BigFloat(x.left_prec, x.left_sign, x.left_exp, x.left_d)
    right = BigFloat(x.right_prec, x.right_sign, x.right_exp, x.right_d)
    return "[$(string(left)), $(string(right))]"
end

print(io::IO, b::Interval) = print(io, string(b))
show(io::IO, b::Interval) = print(io, string(b), " with $(precision(b)) bits of precision")
showcompact(io::IO, b::Interval) = print(io, string(b))

end # module


text: println("=====热烈欢迎群众把玩=====")

max_stars_repo_name: JulialangOrgCN/howtojulia

max_issues_repo_name: JulialangOrgCN/howtojulia

max_forks_repo_name: JulialangOrgCN/howtojulia


text:
using SFML

window = RenderWindow("Line", 800, 600)
event = Event()

p1 = Vector2f(400, 300)
p2 = Vector2f(5, 5)

# Create a line from p1 to p2
line = Line(p1, p2, 10)
set_fillcolor(line, SFML.yellow)
println(get_points(line))
println(get_thickness(line))

while isopen(window)
    while pollevent(window, event)
        if get_type(event) == EventType.CLOSED
            close(window)
        end
    end

    clear(window, SFML.black)
    draw(window, line)
    display(window)
end


text:
using Fabrics

# Example: Weak Inclusion
Δη, ϕ = 1e-3, 0.1
α_bulk = 30.5
p = parameterization(Δη, ϕ)
@show r1, r2 = fabric_shape(10, 1, p)
@show α = fabric_orientation(α_bulk, p)
@show η_cutoff = weakening(Δη, ϕ)

# Example: Strong Inclusion
Δη, ϕ = 1e1, 0.2
α_bulk = 30.5
p = parameterization(Δη, ϕ)
@show r1, r2 = fabric_shape(10, 1, p)
@show α = fabric_orientation(α_bulk, p)


text:
abstract type AbstractServiceFormulation end

mutable struct ServiceModel{D <: PSY.Service,
                            B <: PSI.AbstractServiceFormulation}
    service::Type{D}
    formulation::Type{B}
end

function construct_service!(ps_m::CanonicalModel,
                           service_model::ServiceModel,
                           system_formulation::Type{S},
                           system::PSY.System,
                           time_range::UnitRange{Int64};
                           kwargs...) where {S <: PM.AbstractPowerFormulation}

    construct_service!(ps_m,
                      service_model.service,
                      service_model.formulation,
                      system_formulation,
                      system,
                      time_range;
                      kwargs...)

    return

end