text:
using Tullio, Test, LinearAlgebra

@testset "new arrays" begin

    # functions
    @tullio A[i] := (1:10)[i]^2
    @test A == [i^2 for i in 1:10]

    @tullio A[i] := (1:10)[i] * i
    @test A == [i^2 for i in 1:10]

    # diagonals
    @tullio D[i,i] := trunc(Int, sqrt(A[i]))
    @test D == Diagonal(sqrt.(A))

    # arrays of arrays
    C = [fill(i,3) for i=1:5]
    @tullio M[i,j] := C[i][j]
    @test M == (1:5) .* [1 1 1]

    # fields
    E = [(a=i, b=i^2, c=[i,2i,3i]) for i in 1:10]
    @tullio O[i] := A[i]//E[i].b # avx disabled by try/catch
    @test O == ones(10)

    @tullio F[i,j] := E[i].c[j]
    @test F == (1:10) .* [1 2 3]

    # arrays of tuples
    Y = [(i,i^2,i^3) for i in 1:10]
    @tullio W[i,j] := Y[i][j]
    @test W[9,3] == 9^3

    # scalar
    @tullio S := A[i]/2
    @tullio S′ = A[i]/2 # here = is equivalent
    @test S ≈ S′ ≈ sum(A)/2

    # almost scalar
    @tullio Z[] := A[i] + A[j]
    @test Z isa Array{Int,0}
    @tullio Z′[1,1] := A[i] + A[j]
    @test size(Z′) == (1,1)
    @tullio Z′′[_] := A[i] + A[j]
    @test size(Z′′) == (1,)
    @test Z[] == Z′[1,1] == Z′′[1] == sum(A .+ A')

    # scalar update
    @tullio S += A[i]/2
    @test S ≈ sum(A)

    # fixed
    @tullio F[i] := D[i,5]
    @test F[5] == 5

    j = 6
    @tullio G[i] := D[i,$j]
    @test G[6] == 6

    @test_throws LoadError @eval @tullio D[i,$j] := A[i]

    @tullio H[i] := D[i,:] # storage_type(H, D) == Array, this avoids @avx
    @test H[5] == F

    # trivial dimensions
    @tullio J[1,1,i] := A[i]
    @test size(J) == (1,1,10)

    @tullio J[_,i] := A[i]
    @test J == A'

    # non-unique arrays
    @tullio A2[i] := A[i] + A[i]
    @test A2 == 2 .* A

    # broadcasting
    @tullio S[i] := sqrt.(M[:,i]) # avx & grad now disabled by try/catch
    # @tullio T[i] := A[i] .+ A[j]  # dot does nothing, fails with LoopVectorization loaded

    # scope
    f(x,k) = @tullio y[i] := x[i] + i + $k
    @test f(ones(3),j) == 1 .+ (1:3) .+ j

    g(x) = @tullio y := sqrt(x[i])
    @test g(fill(4,5)) == 10

    # ranges
    @tullio K[i] := i^2  (i ∈ 1:3)
    @test K == (1:3).^2
    @test axes(K,1) === Base.OneTo(3) # literal 1:3

    @tullio N[i,j] := A[i]/j  (j in axes(K,1))  (i in axes(A,1)) # K not an argument
    @test N ≈ A ./ (1:3)'

    @test_throws String @tullio A[i] := i^2 (i in 1+10) # not a range

    # repeated scalar arg
    tri = Base.OneTo(3) # with 1:3, this fails without OffsetArrays,
    # as it only converts shifted indices to OneTo
    @tullio M[i,j] := (r=i, c=j)  (i in tri, j in tri)
    @test M[3,3] == (r=3, c=3)

    # indexing by an array, "gather"...
    J = repeat(1:3, 4);
    @tullio G[i,k] := M[i,J[k]]
    @test G[3,1] == G[3,4] == G[3,7]

    inds = vcat(1:3, 1:3)
    @tullio AI[i] := A[inds[i]]
    @test AI == A[inds]
    jnds = -5:5
    @test_throws String @tullio AJ[j] := A[jnds[j]]
    @test_throws BoundsError A[jnds]
    knds = 1:3.0
    @test_throws String @tullio AK[j] := A[knds[j]]
    @test_throws ArgumentError A[knds]

    # ... and "scatter"
    M = rand(1:99, 4,5)
    J = [3,1,2,3]
    @tullio H[J[i],k] := M[i,k] # i is not marked unsafe, may be threaded
    @test size(H) == (3,5)
    @test H[1,:] == M[2,:] # but H[3,:] gets written into twice.

    J′ = [1,2,10]
    @tullio H′[J′[i'],k] := A[k]
    @test size(H′) == (10, length(A))
    @test H′[2,:] == A
    @test H′[3,4] == 0 # zeroed before being written into

    inds = vcat(1:3, 1:3)
    @test_throws String @tullio H[inds[i],k] := M[i,k] # range of index i

    # masking
    @tullio M[i,j] := A[i] * A[j] * (i<=j)
    @test M == UpperTriangular(A .* A')

    # primes
    @test A == @tullio P[i′] := A[i']
    @test A == @tullio P[i'] := A[i′]
    @test [1,4,9] == @tullio Q[i'] := (i′)^2  (i' in 1:3)

    # non-numeric array
    @tullio Y[i] := (ind=i, val=A[i])
    @test Y[2] === (ind = 2, val = 4)

    # no name given
    Z = @tullio [i] := A[i] + 1
    @test Z == A .+ 1

    # multi-line
    @tullio B[i,j] := begin
        x = (1:10)[i] + 3
        y = (1:3)[j]
        x // y
    end
    @test B == (4:13) .// (1:3)'

    # wrong ndims
    @test_throws Any @tullio Z[i] := B[i] # Any as TensorOperations throws ErrorException

    # internal name leaks
    for sy in Tullio.SYMBOLS
        @test !isdefined(@__MODULE__, sy)
    end

end

@testset "in-place" begin

    A = [i^2 for i in 1:10]
    D = similar(A, 10, 10)

    @tullio D[i,j] = A[i] + 100
    @test D[3,7] == A[3] + 100

    # sum and +=
    B = copy(A);
    D .= 3;
    @tullio B[i] += D[i,j]
    @test B[1] == A[1] + 30

    # writing back into same
    B = copy(A)
    @tullio B[i] += B[i] + 10^3
    @test B[6] == 2 * A[6] + 10^3

    @tullio A[i] = A[i] + 100
    @test A[1] == 101

    # indices in expression
    @tullio A[i] = 100*i
    @test A[7] == 700

    # fixed on left
    j = 3
    @tullio D[$j,i] = 99
    @test D[j,j] == 99
    @test D[1,1] != 0

    # diagonal & ==, from https://github.com/ahwillia/Einsum.jl/pull/14
    B = [1 2 3; 4 5 6; 7 8 9]
    @tullio W[i, j, i, n] := B[n, j]  i in 1:2
    @test size(W) == (2,3,2,3)
    @test W[1,2,1,3] == B[3,2]
    @test W[1,1,2,2] == 0

    W2 = zero(W);
    @tullio W2[i, j, m, n] = (i == m) * B[n, j]
    @test W2 == W

    @test_throws LoadError @eval @tullio [i,j] = A[i] + 100

    # zero off-diagonal? no.
    @tullio D[i,i] = A[i]
    @test D[1,3] != 0

    # scatter operation
    D = similar(A, 10, 10) .= 999
    inds = [2,3,5,2]
    @tullio D[inds[i],j] = A[j]
    @test D[2,:] == A
    @test D[4,4] != 0 # not zeroed before writing.

    @tullio D[inds[i],j] += A[j]
    @test D[2,:] == 3 .* A # was not re-zeroed for +=

    kinds = [1,2,13,4]
    @test_throws String @tullio D[kinds[i],j] = A[j] # BoundsError needs to know which array

    # assignment: no loop over j
    B = zero(A);
    @tullio B[i] = begin
        j = mod(i^4, 1:10)
        A[j]
    end
    @test_skip B == A[[mod(i^4, 1:10) for i in 1:10]]
    # on travis 1.3 multi-threaded, B == [500, 600, 100, 600, 500, 600, 100, 600, 100, 1000]
    # and on 1.4 multi-threaded,    B == [100, 600, 100, 600, 100, 600, 100, 600, 100, 1000]

    # wrong ndims
    @test ndims(B)==1 && ndims(D)==2
    @test_throws Any @tullio B[i] = D[i]^2
    @test_throws Any @tullio D[i] = B[i]+2
    @test_throws Any @tullio B[i,j] = D[i,j]

    # internal name leaks
    for sy in Tullio.SYMBOLS
        @test !isdefined(@__MODULE__, sy)
    end

end

if !@isdefined OffsetArray
    @testset "without packages" begin

        A = [i^2 for i in 1:10]

        # without OffsetArrays
        @test axes(@tullio B[i] := A[2i+1] + A[i]) === (Base.OneTo(4),)
        @test_throws String @tullio C[i] := A[2i+5]

        J = [3,5,7] # doesn't start at 1
        @test_throws String @tullio G[J[i],k] := A[k]

        # without NamedDims
        @test_throws UndefVarError @tullio M[row=i, col=j, i=1] := (1:3)[i] // (1:7)[j]

    end
end

using OffsetArrays

@testset "index shifts" begin

    A = [i^2 for i in 1:10]

    @tullio L[i,j] := A[i]//j  (j ∈ 2:3, i in 1:10) # no shift, just needs OffsetArrays
    @test axes(L) == (1:10, 2:3)

    # shifts
    @tullio B[i] := A[2i+1] + A[i]
    @test axes(B,1) == 1:4 # would be OneTo(4) without OffsetArrays

    @tullio C[i] := A[2i+5]
    @test axes(C,1) == -2:2 # error without OffsetArrays

    j = 7 # interpolation
    @tullio C[i] := A[2i+$j]
    @test axes(C,1) == -3:1

    cee(A) = @tullio C[i] := A[2i+$j] # closure over j
    @test axes(cee(A),1) == -3:1

    @test_throws String @tullio D[i] := A[i] + B[i]
    @tullio D[i] := A[i] + B[i+0] # switches to intersection
    @test axes(D,1) == 1:4

    @test_throws String @tullio M[i,j] := A[i+0]/A[j]  (i ∈ 2:5, j ∈ 2:5) # intersection for i but not j

    @tullio L[i] := A[i+j+1]  (j ∈ -1:1)
    @test axes(L,1) == 1:8

    # negative
    @test eachindex(@tullio F[i] := A[-1i]) == -10:-1
    @test eachindex(@tullio F[i] := A[-i]) == -10:-1
    @test eachindex(@tullio F[i] := A[-i+0]) == -10:-1
    @test eachindex(@tullio F[i] := A[0-i]) == -10:-1

    # non-constant
    @test axes(@tullio I[i,j] := A[i+j] + 0 * B[j]) == (0:6, 1:4)
    @test axes(@tullio I[i,j] := A[j+i+0] + 0 * B[j]) == (0:6, 1:4)
    @test axes(@tullio I[i,j] := A[(j+i)*1] + 0 * B[j]) == (0:6, 1:4)
    @test axes(@tullio I[i,j] := A[2i+j] + 0 * B[j]) == (0:3, 1:4)
    @test axes(@tullio I[i,j] := A[1j+2i] + 0 * B[j]) == (0:3, 1:4)
    @test axes(@tullio I[i,j] := A[i+2j] + 0 * B[j]) == (-1:2, 1:4)
    @test axes(@tullio I[i,j] := A[2i+2j] + 0 * B[j]) == (0:1, 1:4)
    @test axes(@tullio I[i,j] := A[2(i+j)] + 0 * B[j]) == (0:1, 1:4)
    @test axes(@tullio I[i,j] := A[2i-1+2j] + 0 * B[j]) == (0:1, 1:4)
    @test axes(@tullio I[i,j] := A[2i+2j+5] + 0 * B[j]) == (-3:-2, 1:4)
    @test axes(@tullio I[i,j] := A[2i+2j-5] + 0 * B[j]) == (2:3, 1:4)
    @test axes(@tullio I[i,j] := A[2i+2(j-2)-1] + 0 * B[j]) == (2:3, 1:4)
    @test axes(@tullio I[i,j] := A[2(0+i)+(2j-4)-1] + 0 * B[j]) == (2:3, 1:4)

    @test axes(@tullio J[i,j] := A[i-j] + 0 * B[j]) == (5:11, 1:4)
    @test axes(@tullio J[i,j] := A[-(-i+j)] + 0 * B[j]) == (5:11, 1:4)
    @test axes(@tullio J[i,j] := A[-(j-i)] + 0 * B[j]) == (5:11, 1:4)
    @test axes(@tullio J[i,j] := A[-1*(j-i)] + 0 * B[j]) == (5:11, 1:4)
    @test axes(@tullio J[i,j] := A[i+(-j)] + 0 * B[j]) == (5:11, 1:4)
    @test axes(@tullio J[i,j] := A[-j+i] + 0 * B[j]) == (5:11, 1:4)
    @test axes(@tullio J[i,j] := A[-1j+i] + 0 * B[j]) == (5:11, 1:4)
    @test axes(@tullio J[i,j] := A[-1j-(-i)] + 0 * B[j]) == (5:11, 1:4)
    @test axes(@tullio J[i,j] := A[i-2j] + 0 * B[j]) == (9:12, 1:4)
    @test axes(@tullio J[i,j] := A[-2j+i] + 0 * B[j]) == (9:12, 1:4)
    @test axes(@tullio J[i,j] := A[2i-2j] + 0 * B[j]) == (5:6, 1:4)

    @test axes(@tullio I[i,j] := A[i+j÷2] + 0 * B[j]) == (1:8, 1:4)
    @test axes(@tullio I[i,j] := A[i+(j-1)÷2] + 0 * B[j]) == (1:9, 1:4)
    @test axes(@tullio I[i,j] := A[2i+(j-1)÷2] + 0 * B[j]) == (1:4, 1:4)
    @test axes(@tullio I[i,j] := A[i+(j-1)÷3] + 0 * B[j]) == (1:9, 1:4)

    @test_throws LoadError @eval @tullio I[i,j] := A[i+j] # under-specified

    # in-place
    E = zero(A)
    @tullio E[i] = A[i+5] + 100
    @test E == vcat(A[6:end] .+ 100, zeros(Int,5))

    M = fill(pi/2, 10, 10)
    @tullio M[i,i] = A[i-2]
    @test M[3,3] == A[1]
    @test M[1,1] == pi/2 # was not set to zero

    # shifts on left
    E = zero(A)
    @tullio E[2i+1] = A[i]
    @test E[2+1] == A[1]
    @test E[2*4+1] == A[4]

    # non-constant
    @tullio I[i,j] := 0 * A[i+j] + 0 * B[j]
    @test axes(@tullio I[i,j] = A[i+j] + B[j]) == (0:6, 1:4) # over-specified
    @test axes(@tullio I[i,j] = A[i+j]) == (0:6, 1:4) # needs range from LHS

    # indexing by an array
    inds = [-1,0,0,0,1]
    @tullio K[i,j] := A[inds[i]+j]
    @test K[2,3] == K[3,3] == K[4,3]
    @test axes(K) == (1:5, 2:9)

    @tullio K2[i,j] := A[j+2inds[i]]
    @test axes(K2) == (1:5, 3:8)

    j = 7
    @test_skip @tullio K3[i,j] := A[j+2inds[i]+$j]
    @test_broken vec(K2) == vec(K3)

    # scatter with shift not allowed
    @test_throws LoadError @eval @tullio G[inds[i]+1, j] := A[j]
    @test_throws LoadError @eval @tullio G[2inds[i], j] := A[j]

    # multiplication not implemented
    @test_throws LoadError @eval @tullio C[i] = A[i*j] + A[i]
    @test_throws LoadError @eval @tullio C[i] = A[i⊗j] + A[i]
    @test_throws LoadError @eval @tullio C[i] = A[(i,j)] + A[i]

    # magic shift
    @test axes(@tullio Z[i+_] := A[2i+10]) === (Base.OneTo(5),)

    @test_throws LoadError @eval @tullio Z[_+i] := A[2i+10] # wrong notation
    @test_throws LoadError @eval @tullio Z[J[i]+_] := A[2i+10] # with scatter
    @test_throws LoadError @eval @tullio Z[i+_] = A[2i+10] # in-place
end

@testset "modulo, clamped & padded" begin

    A = [i^2 for i in 1:10]
    B = 1:5

    @test vcat(B,B) == @tullio C[i] := B[mod(i)]  i in 1:10
    @test vcat(B, fill(B[end],5)) == @tullio D[i] := min(A[i], B[clamp(i)])
    @test [4,16,36,64,100,4] == @tullio E[i] := A[mod(2i)]  i in 1:6

    @test vcat(zeros(5), B, zeros(5)) == @tullio C[i] := B[pad(i-5,5)]
    @test vcat(zeros(2), A, zeros(3)) == @tullio D[i+_] := A[pad(i,2,3)]
    @test vcat(A, zeros(10)) == @tullio E[i] := A[pad(i)]  i in 1:20

    # pairconstraints
    @tullio F[i] := A[mod(i+k)] * ones(3)[k]  (i in axes(A,1))
    @test F[end] == 1 + 2^2 + 3^2
    @tullio F[i] = A[clamp(i+k)] * ones(7)[k]

    @tullio G[i] := A[pad(i+k, 4)] * ones(3)[k]  pad=100
    @test axes(G,1) == -4:11
    @test G[-4] == G[11] == 300

    # matrix
    M = rand(Int8, 3,4)
    @tullio H[i+_,j+_] := M[pad(i,2), pad(j,3)]  pad=1
    @test H == [trues(2,10); trues(3,3) M trues(3,3); trues(2,10)]

    # pad keyword
    @tullio J[i,i] := sqrt(A[i])  pad=-1
    @test J[3,4] == -1
    @test diag(J) == 1:10

    # unable to infer range
    @test_throws LoadError @eval @tullio F[i] := A[mod(i+1)]
    @test_throws LoadError @eval @tullio F[i] := A[pad(2i)]
    # can't use index mod(i) on LHS
    @test_throws LoadError @eval @tullio G[mod(i)] := A[i]
    # not sure what to do with clamp(i), sorry
    @test_throws LoadError @eval @tullio F[i] := A[clamp(i)+1]
    # eltype of pad doesn't fit
    @test_throws InexactError @tullio H[i] := A[pad(i,3)]  pad=im
    @test_throws InexactError @tullio J[i,i] := A[i]  pad=im
end

@testset "other reductions" begin

    A = [i^2 for i in 1:10]

    # basics
    @test [prod(A)] == @tullio (*) P[_] := float(A[i])
    @test maximum(A) == @tullio (max) m := float(A[i])
    @test minimum(A) == @tullio (min) m := float(A[i]) # fails with @avx

    @test true == @tullio (&) p := A[i] > 0
    @test true === @tullio (&) p := A[i] > 0
    @test true == @tullio (|) q := A[i] > 50

    # in-place
    C = copy(A)
    @test cumprod(A) == @tullio (*) C[k] = ifelse(i<=k, A[i], 1)
    @test cumprod(A).^2 == @tullio (*) C[k] *= i<=k ? A[i] : 1

    M = rand(1:9, 4,5)
    @test vec(prod(M,dims=2)) == @tullio (*) B[i] := M[i,j]

    # ^= generalises +=, *=
    C = copy(A)
    @tullio (max) C[i] ^= 5i
    @test C == max.(5:5:50, A)
    @test_throws LoadError @eval @tullio A[i] ^= A[i]
    @test_throws LoadError @eval @tullio (*) A[i] ^= A[i]

    # initialisation
    @test 200 == @tullio (max) m := A[i] init=200
    @tullio (max) C[i] := i^2   (i in 1:10, j in 1:1)  init=33.3 # widens type
    @test C == max.(33.3, A)
    @tullio C[i] := 0   (i in 1:10, j in 1:1)  init=randn() tensor=false # runs once
    @test C == fill(C[1], 10)

    # more dimensions
    Q = rand(1:10^3, 4,5,6)
    @test vec(maximum(Q,dims=(2,3))) == @tullio (max) R[i] := Q[i,j,k]
    @test vec(minimum(Q,dims=(1,3))).+2 == @tullio (min) P[j] := Q[i,j,k]+2
    @test dropdims(maximum(Q, dims=2), dims=2) == @tullio (max) S[i,k] := Q[i,j,k]

    # indexing
    ind = vcat(1:3, 1:3)
    V = 1:6
    @tullio (*) Z[j] := M[ind[k],j] * exp(-V[k]) # product over k
    @test Z ≈ vec(prod(M[ind,:] .* exp.(.-V), dims=1))

    # scalar update ("plusequals" internally)
    s = 1.0
    @tullio (*) s *= float(A[i])
    @test s == prod(A)
    @tullio s *= float(A[i]) # works without specifying (*), is this a good idea?
    @test s == float(prod(A))^2

    @test_throws LoadError @eval @tullio s += (*) A[i] # should be *=
    @test_throws LoadError @eval @tullio s *= (max) A[i] # should be ^=

    # scalar + threading
    L = randn(100 * Tullio.TILE[]);
    @tullio (max) m := L[i]
    @test m == maximum(L)

    # ... with a weird init, result would be unpredictable, hence an error:
    @test_throws String @tullio s2 := A[i]^2 init=2 # at runtime
    @test sum(A.^2)+2 == @tullio s2 := A[i]^2 init=2 threads=false # is OK

    # promotion of init & += cases:
    B = rand(10)
    @test sum(B.^2)+2 ≈ @tullio s2 := B[i]^2 init=2 threads=false
    s3 = 3
    @test sum(B.^2)+3 ≈ @tullio s3 += B[i]^2
    s4 = 4im
    @test sum(B.^2)+4im ≈ @tullio s4 += B[i]^2

    # no reduction means no redfun, and no init:
    @test_throws LoadError @eval @tullio (max) A2[i] := A[i]^2
    @test_throws LoadError @eval @tullio A2[i] := A[i]^2 init=0.0

end

@testset "finalisers" begin

    A = [i^2 for i in 1:10]

    @tullio B[i,j] := A[i] + A[k] // A[j]

    @tullio B2[_,j] := (B[i,j] + B[j,i])^2 |> sqrt
    @test B2 ≈ mapslices(norm, B .+ B', dims=1)

    # trivial use, scalar output -- now forbidden
    @test_throws LoadError @eval @tullio n2 = A[i]^2 |> sqrt

    # trivial use, no reduction -- now forbidden
    @test_throws LoadError @eval @tullio A2[i] := A[i]^2 |> sqrt
    @test_throws LoadError @eval @tullio (*) A2[i] := A[i]^2 |> sqrt

    # larger size, to trigger threads & tiles
    C = randn(10^6) # > Tullio.BLOCK[]
    @tullio n2[_] := C[i]^2 |> sqrt
    @test n2[1] ≈ norm(C,2)

    D = rand(1000, 1000) # > Tullio.TILE[]
    @tullio D2[_,j] := D[i,j]^2 |> sqrt
    @test D2 ≈ mapslices(norm, D, dims=1)

    # functions with underscores
    @tullio n2′[] := A[i]^2 |> (_)^0.5
    @test n2′[] ≈ norm(A,2)

    @tullio (max) E[i] := float(B[i,j]) |> atan(_, A[i]) # i is not reduced over
    @test E ≈ vec(atan.(maximum(B, dims=2), A))

    j = 2
    @tullio G[i'] := float(B[i',j]) |> atan(_, B[i',$j])
    @test G ≈ vec(atan.(sum(B, dims=2), B[:,j]))

    @test_throws LoadError @eval @tullio F[i] := B[i,j] |> (_ / A[j]) # wrong index
    C = randn(10^6)
    @test_throws String @tullio F[i] := B[i,j] |> (_ / C[i]) # wrong length

end

@testset "named dimensions" begin

    using NamedDims

    # reading
    N = NamedDimsArray(rand(Int8,3,10), (:r, :c))

    @tullio A[i,j] := N[i, j] + 100 * (1:10)[j]
    @test A == N .+ 100 .* (1:10)'

    @tullio B[i] := N[r=i, c=1]
    @test B == N[:,1]

    @tullio C[j,i] := N[c=j, r=i] + 100 * (1:10)[j]
    @test A == C'
    @test dimnames(C) == (:_, :_) # similar(parent(A)) avoids a bug

    # writing
    @tullio M[row=i, col=j, i=1] := (1:3)[i] // (1:7)[j]
    @test dimnames(M) == (:row, :col, :i)

end

@testset "options" begin

    # keyword threads accepts false or a positive integer
    @tullio A[i] := (1:10)[i]^2  threads=false
    @tullio A[i] := (1:10)[i]^2  threads=2^2
    # when using KernelAbstractions, something leaks from the 1st leading 2nd to error
    block = 64
    @tullio A[i] := (1:10)[i]^2  threads=block # Symbol
    @test_throws LoadError @macroexpand1 @tullio A[i] := (1:10)[i]^2  threads=:maybe

    # keyword verbose accepts values [true, false, 2, 3]
    @tullio A[i] := (1:10)[i]^2  verbose=1
    @tullio A[i] := (1:10)[i]^2  verbose=false
    @test_throws LoadError @macroexpand1 @tullio A[i] := (1:10)[i]^2  verbose=4

    # keyword grad accepts values [false, Base, Dual]
    @tullio A[i] := (1:10)[i]^2  grad=false
    @tullio A[i] := (1:10)[i]^2  grad=Base
    @test_throws LoadError @macroexpand1 @tullio A[i] := (1:10)[i]^2  grad=true

    # recognised keywords are [:threads, :verbose, :avx, :cuda, :grad]
    @test_throws LoadError @macroexpand1 @tullio A[i] := (1:10)[i]^2  key=nothing

end

@testset "bugs" begin

    # https://github.com/mcabbott/Tullio.jl/issues/10
    arr = [1 2; 3 4]
    function f10(arr)
        @tullio res1 = arr[i, k] - arr[i - 1, k]
        @tullio res2 = arr[i, k] - arr[i, k + 1]
        return res1 + res2
    end
    @test f10(arr) == 2

    let
        B = rand(3,3)
        @tullio tot = B[i, k] - B[i - 1, k]
        @test_throws UndefVarError 𝒜𝒸𝓉! isa Function
    end

    # https://github.com/mcabbott/Tullio.jl/issues/35
    a = [1 2 3; 4 5 6];
    b = [10,20,30];
    @test sum(b) == @tullio s := b[a[1,i]]

    # https://github.com/mcabbott/Tullio.jl/issues/36
    # final type real, intermediate complex... not fixed yet!
    xs = randn(1000)
    @test_throws InexactError @tullio z[i] := exp(im * xs[i] - xs[j]) |> abs2

    # https://github.com/mcabbott/Tullio.jl/issues/43
    P = rand(2,2,3); Diff = rand(3,3); n=4
    @test axes(
        @tullio dP[x,y,z] := Diff[a+2, b+2] * Diff[c+2, d+2] *
            P[mod(x+a+c), mod(y+b+d), z] * P[mod(x+a),mod(y+b),z] (a in -1:1,
            b in -1:1, c in -1:1, d in -1:1, z in 1:3, x in 1:n, y in 1:n)
        ) == (1:4, 1:4, 1:3)
    @test axes(
        @tullio out[x,y] := Diff[mod(x+(y+2z)),x] * Diff[y,clamp(x+1+2x+y)] z in 1:3
        ) == (1:3, 1:3)
    # unable to infer range of index z
    @test_throws LoadError @eval @tullio out[x,y] := Diff[mod(x+(y+2z)),x] * Diff[x,y]

    # https://discourse.julialang.org/t/unexpected-tullio-behavior/49371/6
    x = [rand(3), rand(2)]
    @test_throws String @tullio y[i,j] := x[i][j]
    n = 1; a = ones(4) # should now throw "elements x[\$n] must be of uniform size"
    @test_throws String @tullio a[j] = x[$n][j+0]  j in 1:3

    xt = [(a=1, b=rand(3)), (a=1, b=rand(2))] # version with field access
    @test_throws String @tullio y[i,j] := xt[i].b[j]

    @tullio a[j] = $(x[n])[j+0]  j in 1:3 # interpolation of an expression
    @test a == vcat(x[n],1) # requires postwalk -> prewalk

    # https://github.com/mcabbott/Tullio.jl/issues/46
    @test 1:4 == let s = [1,2,3,4], ndims=maximum(s), axes=size, eachindex=lastindex
       @tullio x[i] := s[i]
    end
    @test_broken [3,7] == let s = [1 2;3 4], zero=one # left because of #50
       @tullio x[i] := s[i,j]
    end

end


text:
x0 = 1
y0 = 2
n = 20
r0 = 5 .+ randn(n)/10
p = r0 .*exp.(im.*rand(0:0.1:2π,n))
x = real(p) .+ x0
y = imag(p) .+ y0
(a,b,r) =  CircleFit.taubin(x,y)
@test isapprox(a,x0,atol=0.1)
@test isapprox(b,y0,atol=0.1)
@test isapprox(r,5,atol=0.1)

ext: jl


text:
# module escape_dimer_one_hit
# export escape_exit_num_dimer
using DifferentialEquations
using Roots
using Printf
# using MATLAB
# using ProgressMeter

include("leap_frog_definitions.jl")
max_hit_q1_back = 1
barrier = 5

condition_max_hits_back_one(u, t, integrator) = u[3] == max_hit_q1_back
affect_stop!(integrator) = terminate!(integrator)
function condition_hits_PSS_q1(u, t, integrator) # Event when event_f(u,t) == 0
    u[1]
end

function affect_update_iterator_q1!(integrator)
    integrator.u[3] = integrator.u[3] + 1
end

callback_max_hits_back_one = DiscreteCallback(condition_max_hits_back_one, affect_stop!)
callback_hits_PSS_q1 = ContinuousCallback(condition_hits_PSS_q1,nothing,affect_update_iterator_q1!)
cb_neg = CallbackSet(callback_hits_PSS_q1, callback_max_hits_back_one)

function back_one(Q2,P2, t_end, H)
    P1 = P1_poly(Q2, P2, H)
    if ~isempty(P1)
        q0, p0 = [zeros(3) for i in 1:2]
        q0[1] = 0
        q0[2] = Q2
        q0[3] = 0
        p0[1] = P1[1]
        p0[2] = P2
        p0[3] = 0
       #constructor for ODE
        prob = HamiltonianProblem{true}(Hamiltonian_Dimer_Backwards, q0, p0, (0., t_end));
        #solve ode , save_everystep=false is important to prevent sol to include all points, not just event points
        # sol=solve(prob, RK4(),maxiters=1e20,callback=cb,save_start=false,save_end=true,save_everystep=false)

        # sol=solve(prob, RK4(),maxiters=1e20,reltol=1e-8,abstol=1e-10,callback=cb,save_start=true,save_end=true,save_everystep=false)
        sol=solve(prob, Vern9(),maxiters=1e20, reltol=1e-13,abstol=1e-15,callback=cb_neg,save_start=false,save_end=true,save_everystep=false)

        uf=zeros(6)
        uf[1]=sol[1,end]
        uf[2]=sol[2,end]
        uf[3]=sol[3,end]
        uf[4]=sol[4,end]
        uf[5]=sol[5,end]
        uf[6]=sol[6,end]
        q1=uf[1]
        q2=uf[2]
        p1=uf[4]
        p2=uf[5]
        q=[q1, q2]
        p=[p1, p2]
        diff=abs(Hamiltonian_Dimer(q,p,0)-H)

        if uf[3]==1 && diff<1e-10
            return q1,q2,p1,p2

        end
        # if dH<1e-5
        #     if is_it(uf,H,tol_dist)
        #         return sol.u[end][5]
        #     else
        #         return -1
        #     end
        # end
    end
end


text:
using LinearRegression
using Test

function check_type_stability(regressor, X)
    @testset "type-inference" begin
        @inferred coef(regressor)
        @inferred regressor(X)
        @inferred regressor(X[1, :])
    end
end

function test_linreg_multivariate(X, y, beta, X_element, y_element)
    intercept = size(X, 2) < size(beta, 1)

    @testset "weights $(repr(weights))" for weights in (nothing, ones(size(X, 1)))
        @testset "Solver $method" for method in (SolveQR(), SolveCholesky())
            regressor = linregress(X, y, weights; intercept=intercept, method=method)
            @test coef(regressor) ≈ beta
            @test regressor(X) ≈ y
            @test regressor(X_element) ≈ y_element
            check_type_stability(regressor, X)
        end
    end
end

@testset "LinearRegression.jl" begin
    @testset "univariate" begin
        X = [4 6]'
        y = [-2, 4]

        @testset "weights $(repr(weights))" for weights in (nothing, ones(2))
            @testset "Solver $method" for method in (SolveQR(), SolveCholesky())
                regressor = linregress(X, y, weights; intercept=true, method=method)
                @test regressor([5]) ≈ 1
                @test regressor([4 5 6]') ≈ [-2, 1, 4]
                @test coef(regressor) ≈ [3, -14]

                regressor = linregress(X, y, weights; intercept=false, method=method)
                @test regressor([0]) == 0
                check_type_stability(regressor, [4 5 6]')
            end
        end
    end

    @testset "multivariate" begin
        @testset "intercept $intercept" for intercept in (true, false)
            X = rand(100, 13)
            X0 = intercept ? [X ones(size(X, 1))] : X

            @testset "single-output" begin
                beta = rand(size(X0, 2))
                y = X0 * beta

                test_linreg_multivariate(X, y, beta, X[17, :], y[17])
            end

            @testset "multi-output" begin
                beta = rand(size(X0, 2), 5)
                y = X0 * beta

                test_linreg_multivariate(X, y, beta, X[17, :], y[17, :])
            end
        end
    end
end


text:
using APipe
using Test

function λ𝟎(; keyword = "default")
    keyword
end

function λ𝟏(arg; keyword = "default")
    arg, keyword
end

function λ(arg, default = "default"; keyword = "default")
    arg, default, keyword
end

@testset "APipe" begin

@testset "@> macro" begin
    @testset "applied to function" begin
        @test ("chained") == @> "chained" |> (:keyword, λ𝟎)
        @test ("chained", "default") == @> "chained" |> λ𝟏
        @test ("chained", "default") == @> "chained" |> (1, λ𝟏)
    end

    @testset "applied to function call" begin
        @test ("chained", "default", "default") == @> "chained" |> (1, λ)
        @test ("chained", "default", "default") == @> "chained" |> λ()
        @test ("chained", "default", "default") == @> "chained" |> (1, λ())
        @test ("chained", "passed", "default") == @> "chained" |> λ("passed")
        @test ("passed", "chained", "default") == @> "chained" |> (2, λ("passed"))
        @test ("passed", "default", "chained") == @> "chained" |> (:keyword, λ("passed"))
    end

    @testset "applied with dot fuse" begin
        strarray = [rand('A':'Z', 5) |> join for i ∈ 1:10]
        @test (strarray .== @> strarray .|> λ() .|> getindex(1)) |> all
    end
end

@testset "@>> macro" begin
    @testset "usual |> chaining" begin
        @test ("chained") == @>> "chained" x -> λ𝟎(; keyword = x)
        @test ("chained", "default") == @>> "chained" λ𝟏
        @test ("chained", "default") == @>> "chained" x -> λ𝟏(x)
    end

    @testset "@> chanining" begin
        @test ("chained", "default", "default") == @>> "chained" (1, λ)
        @test ("chained", "default", "default") == @>> "chained" λ()
        @test ("passed", "chained", "default") == @>> "chained" (2, λ("passed"))
        @test ("passed", "default", "chained") == @>> "chained" (:keyword, λ("passed"))
    end
end

@testset "@.>> macro" begin
    strarray = [rand('A':'Z', 5) |> join for i ∈ 1:10]
    @test (strarray .== @.>> strarray λ() getindex(1)) |> all
end

end


text:
function prepareRow(A::ITensor, op::ITensor, left_A::ITensor, right_A::ITensor, left_I::ITensor, right_I::ITensor, col::Int, Nx::Int)
    AA = A * op * dag(A')
    if col > 1
        ci  = commonindex(A, left_A)
        msi = multiply_side_ident(AA, ci, left_I)
        AA *= msi
    end
    if col < Nx
        ci  = commonindex(A, right_A)
        msi = multiply_side_ident(AA, ci, right_I)
        AA *= msi
    end
    return AA
end

function makeAncillaryIs(A::PEPS, L::Environments, R::Environments, col::Int)
    Ny, Nx   = size(A)
    is_cu   = is_gpu(A)
    dummy    = is_cu    ? cuITensor(1.0)  : ITensor(1.0)
    left_As  = [col > 1  ? A[row, col - 1] : dummy for row in 1:Ny]
    right_As = [col < Nx ? A[row, col + 1] : dummy for row in 1:Ny]
    col_site_inds = [findindex(x, "Site") for x in A[:, col]]
    ops = map(x -> spinI(x; is_gpu=is_cu), col_site_inds)
    AAs = [prepareRow(A[row, col], ops[row], left_As[row], right_As[row], L.I[row], R.I[row], col, Nx) for row in 1:Ny]
    Iabove = cumprod(reverse(AAs))
    return Iabove
end

function makeAncillaryIsBelow(A::PEPS, L::Environments, R::Environments, col::Int)
    Ny, Nx   = size(A)
    is_cu   = is_gpu(A)
    dummy    = is_cu    ? cuITensor(1.0)  : ITensor(1.0)
    left_As  = [col > 1  ? A[row, col - 1] : dummy for row in 1:Ny]
    right_As = [col < Nx ? A[row, col + 1] : dummy for row in 1:Ny]
    col_site_inds = [findindex(x, "Site") for x in A[:, col]]
    ops = map(x -> spinI(x; is_gpu=is_cu), col_site_inds)
    AAs = [prepareRow(A[row, col], ops[row], left_As[row], right_As[row], L.I[row], R.I[row], col, Nx) for row in 1:Ny]
    return cumprod(AAs)
end

function updateAncillaryIs(A::PEPS, Ibelow::Vector{ITensor}, L::Environments, R::Environments, row::Int, col::Int )
    Ny, Nx  = size(A)
    is_cu  = is_gpu(A)
    dummy   = is_cu   ? cuITensor(1.0)  : ITensor(1.0)
    left_A  = col > 1  ? A[row, col - 1] : dummy
    right_A = col < Nx ? A[row, col + 1] : dummy
    op      = spinI(findindex(A[row, col], "Site"); is_gpu=is_cu)
    AA      = prepareRow(A[row, col], op, left_A, right_A, L.I[row], R.I[row], col, Nx)
    AA     *= row > 1  ? Ibelow[row - 1] : dummy
    Ibelow[row] = AA
    return Ibelow
end

function makeAncillaryFs(A::PEPS, L::Environments, R::Environments, H, col::Int)
    Ny, Nx   = size(A)
    is_cu   = is_gpu(A)
    dummy    = is_cu    ? cuITensor(1.0)  : ITensor(1.0)
    left_As  = [col > 1  ? A[row, col - 1] : dummy for row in 1:Ny]
    right_As = [col < Nx ? A[row, col + 1] : dummy for row in 1:Ny]
    col_site_inds = [findindex(x, "Site") for x in A[:, col]]
    Fabove   = fill(Vector{ITensor}(), length(H))
    for opcode in 1:length(H)
        op_row      = H[opcode].sites[1][1]
        ops = map(x -> spinI(x; is_gpu=is_cu), col_site_inds)
        ops[op_row] = replaceindex!(copy(H[opcode].ops[1]), H[opcode].site_ind, col_site_inds[op_row])
        ops[op_row] = replaceindex!(copy(ops[op_row]), H[opcode].site_ind', col_site_inds[op_row]')
        ancFs = [prepareRow(A[row, col], ops[row], left_As[row], right_As[row], L.I[row], R.I[row], col, Nx) for row in 1:Ny]
        Fabove[opcode] = cumprod(reverse(ancFs))
    end
    return Fabove
end

function updateAncillaryFs(A::PEPS, Fbelow, Ibelow::Vector{ITensor}, L::Environments, R::Environments, H, row::Int, col::Int)
    Ny, Nx   = size(A)
    is_cu   = is_gpu(A)
    dummy    = is_cu    ? cuITensor(1.0)  : ITensor(1.0)
    left_As  = [col > 1  ? A[row, col - 1] : dummy for row in 1:Ny]
    right_As = [col < Nx ? A[row, col + 1] : dummy for row in 1:Ny]
    col_site_inds = [findindex(x, "Site") for x in A[:, col]]
    for opcode in 1:length(H)
        op_row      = H[opcode].sites[1][1]
        ops         = ITensor[spinI(spin_ind; is_gpu=is_cu) for spin_ind in col_site_inds]
        if op_row == row
            ops[op_row] = replaceindex!(copy(H[opcode].ops[1]), H[opcode].site_ind, col_site_inds[op_row])
            ops[op_row] = replaceindex!(ops[op_row], H[opcode].site_ind', col_site_inds[op_row]')
        end
        ancF = prepareRow(A[row, col], ops[row], left_As[row], right_As[row], L.I[row], R.I[row], col, Nx)
        if length(Fbelow[opcode]) > 0
            push!(Fbelow[opcode], Fbelow[opcode][end]*ancF)
        else
            push!(Fbelow[opcode], ancF)
        end
    end
    return Fbelow
end

function makeAncillaryVs(A::PEPS, L::Environments, R::Environments, H, col::Int)
    Ny, Nx   = size(A)
    is_cu   = is_gpu(A)
    dummy    = is_cu    ? cuITensor(1.0)  : ITensor(1.0)
    left_As  = [col > 1  ? A[row, col - 1] : dummy for row in 1:Ny]
    right_As = [col < Nx ? A[row, col + 1] : dummy for row in 1:Ny]
    col_site_inds = [findindex(x, "Site") for x in A[:, col]]
    Vabove   = fill(Vector{ITensor}(), length(H))
    for opcode in 1:length(H)
        op_row_a      = H[opcode].sites[1][1]
        op_row_b      = H[opcode].sites[2][1]
        ops           = map(x->spinI(x; is_gpu=is_cu), col_site_inds)
        ops[op_row_a] = replaceindex!(copy(H[opcode].ops[1]), H[opcode].site_ind, col_site_inds[op_row_a])
        ops[op_row_a] = replaceindex!(ops[op_row_a], H[opcode].site_ind', col_site_inds[op_row_a]')
        ops[op_row_b] = replaceindex!(copy(H[opcode].ops[2]), H[opcode].site_ind, col_site_inds[op_row_b])
        ops[op_row_b] = replaceindex!(ops[op_row_b], H[opcode].site_ind', col_site_inds[op_row_b]')
        ancVs         = [prepareRow(A[row, col], ops[row], left_As[row], right_As[row], L.I[row], R.I[row], col, Nx) for row in 1:Ny]
        Vabove[opcode] = cumprod(reverse(ancVs))#[1:op_row_a+1]
    end
    return Vabove
end

function makeAncillaryVsBelow(A::PEPS, L::Environments, R::Environments, H, col::Int)
    Ny, Nx   = size(A)
    is_cu   = is_gpu(A)
    dummy    = is_cu    ? cuITensor(1.0)  : ITensor(1.0)
    left_As  = [col > 1  ? A[row, col - 1] : dummy for row in 1:Ny]
    right_As = [col < Nx ? A[row, col + 1] : dummy for row in 1:Ny]
    col_site_inds = [findindex(x, "Site") for x in A[:, col]]
    Vbelow   = fill(Vector{ITensor}(), length(H))
    for opcode in 1:length(H)
        op_row_a      = H[opcode].sites[1][1]
        op_row_b      = H[opcode].sites[2][1]
        ops           = map(x->spinI(x; is_gpu=is_cu), col_site_inds)
        ops[op_row_a] = replaceindex!(copy(H[opcode].ops[1]), H[opcode].site_ind, col_site_inds[op_row_a])
        ops[op_row_a] = replaceindex!(ops[op_row_a], H[opcode].site_ind', col_site_inds[op_row_a]')
        ops[op_row_b] = replaceindex!(copy(H[opcode].ops[2]), H[opcode].site_ind, col_site_inds[op_row_b])
        ops[op_row_b] = replaceindex!(ops[op_row_b], H[opcode].site_ind', col_site_inds[op_row_b]')
        ancVs = [prepareRow(A[row, col], ops[row], left_As[row], right_As[row], L.I[row], R.I[row], col, Nx) for row in 1:Ny]
        Vbelow[opcode] = cumprod(ancVs)#[1:op_row_a+1]
    end
    return Vbelow
end

function updateAncillaryVs(A::PEPS, Vbelow, Ibelow::Vector{ITensor}, L::Environments, R::Environments, H, row::Int, col::Int)
    Ny, Nx   = size(A)
    is_cu   = is_gpu(A)
    dummy    = is_cu    ? cuITensor(1.0)  : ITensor(1.0)
    left_As  = [col > 1  ? A[row, col - 1] : dummy for row in 1:Ny]
    right_As = [col < Nx ? A[row, col + 1] : dummy for row in 1:Ny]
    col_site_inds = [findindex(x, "Site") for x in A[:, col]]
    for opcode in 1:length(H)
        op_row_a      = H[opcode].sites[1][1]
        op_row_b      = H[opcode].sites[2][1]
        ops           = map(x->spinI(x; is_gpu=is_cu), col_site_inds)
        ops[op_row_a] = replaceindex!(copy(H[opcode].ops[1]), H[opcode].site_ind, col_site_inds[op_row_a])
        ops[op_row_a] = replaceindex!(ops[op_row_a], H[opcode].site_ind', col_site_inds[op_row_a]')
        ops[op_row_b] = replaceindex!(copy(H[opcode].ops[2]), H[opcode].site_ind, col_site_inds[op_row_b])
        ops[op_row_b] = replaceindex!(ops[op_row_b], H[opcode].site_ind', col_site_inds[op_row_b]')
        if op_row_b < row
            AA  = prepareRow(A[row, col], ops[row], left_As[row], right_As[row], L.I[row], R.I[row], col, Nx)
            push!(Vbelow[opcode], Vbelow[opcode][end] * AA)
        elseif op_row_b == row
            Ib  = op_row_a > 1 ? Ibelow[op_row_a - 1] : dummy
            AA  = prepareRow(A[op_row_a, col], ops[op_row_a], left_As[op_row_a], right_As[op_row_a], L.I[op_row_a], R.I[op_row_a], col, Nx)
            AA *= prepareRow(A[op_row_b, col], ops[op_row_b], left_As[op_row_b], right_As[op_row_b], L.I[op_row_b], R.I[op_row_b], col, Nx)
            push!(Vbelow[opcode], AA*Ib)
        end
    end
    return Vbelow
end

function makeAncillarySide(A::PEPS, EnvIP::Environments, EnvIdent::Environments, H, col::Int, side::Symbol)
    Ny, Nx   = size(A)
    is_cu   = is_gpu(A)
    col_site_inds = [findindex(x, "Site") for x in A[:, col]]
    Sabove   = fill(Vector{ITensor}(), length(H))
    next_col = side == :left ? col + 1 : col - 1
    for opcode in 1:length(H)
        op_row      = side == :left ? H[opcode].sites[2][1] : H[opcode].sites[1][1]
        ops         = map(x->spinI(x; is_gpu=is_cu), col_site_inds)
        ops         = is_cu ? map(cuITensor, ops) : ops
        this_op     = side == :left ? H[opcode].ops[2] : H[opcode].ops[1]
        ops[op_row] = replaceindex!(copy(this_op), H[opcode].site_ind, col_site_inds[op_row])
        ops[op_row] = replaceindex!(ops[op_row], H[opcode].site_ind', col_site_inds[op_row]')
        AAs         = [A[row, col] * ops[row] * dag(A[row, col])' * EnvIP.InProgress[row, opcode] for row in 1:Ny]
        if (col > 1 && side == :right) || (col < Nx && side == :left)
            cis = [commonindex(A[row, col], A[row, next_col]) for row in 1:Ny]
            msi = [multiply_side_ident(AAs[row], cis[row], EnvIdent.I[row]) for row in 1:Ny]
            AAs = AAs .* msi
        end
        Sabove[opcode] = cumprod(reverse(AAs))
    end
    return Sabove
end

function makeAncillarySideBelow(A::PEPS, EnvIP::Environments, EnvIdent::Environments, H, col::Int, side::Symbol)
    Ny, Nx   = size(A)
    is_cu   = is_gpu(A)
    col_site_inds = [findindex(x, "Site") for x in A[:, col]]
    Sbelow   = fill(Vector{ITensor}(), length(H))
    next_col = side == :left ? col + 1 : col - 1
    for opcode in 1:length(H)
        op_row      = side == :left ? H[opcode].sites[2][1] : H[opcode].sites[1][1]
        ops         = map(x->spinI(x; is_gpu=is_cu), col_site_inds)
        ops         = is_cu ? map(cuITensor, ops) : ops
        this_op     = side == :left ? H[opcode].ops[2] : H[opcode].ops[1]
        ops[op_row] = replaceindex!(copy(this_op), H[opcode].site_ind, col_site_inds[op_row])
        ops[op_row] = replaceindex!(ops[op_row], H[opcode].site_ind', col_site_inds[op_row]')
        AAs         = [A[row, col] * ops[row] * dag(A[row, col])' * EnvIP.InProgress[row, opcode] for row in 1:Ny]
        if (col > 1 && side == :right) || (col < Nx && side == :left)
            cis = [commonindex(A[row, col], A[row, next_col]) for row in 1:Ny]
            msi = [multiply_side_ident(AAs[row], cis[row], EnvIdent.I[row]) for row in 1:Ny]
            AAs = AAs .* msi
        end
        Sbelow[opcode] = cumprod(AAs)
    end
    return Sbelow
end

function updateAncillarySide(A::PEPS, Sbelow, Ibelow::Vector{ITensor}, EnvIP::Environments, EnvIdent::Environments, H, row::Int, col::Int, side::Symbol)
    Ny, Nx   = size(A)
    is_cu    = is_gpu(A)
    col_site_inds = [findindex(x, "Site") for x in A[:, col]]
    next_col = side == :left ? col + 1 : col - 1
    prev_col = side == :left ? col + 1 : col - 1
    for opcode in 1:length(H)
        op_row        = side == :left ? H[opcode].sites[2][1] : H[opcode].sites[1][1]
        ops           = ITensor[spinI(spin_ind; is_gpu=is_cu) for spin_ind in col_site_inds]
        ops           = is_cu ? map(cuITensor, ops) : ops
        this_op       = side == :left ? H[opcode].ops[2] : H[opcode].ops[1]
        ops[op_row]   = replaceindex!(copy(this_op), H[opcode].site_ind, col_site_inds[op_row])
        ops[op_row]   = replaceindex!(ops[op_row], H[opcode].site_ind', col_site_inds[op_row]')
        AA            = A[row, col] * ops[row] * dag(A[row, col]')
        if (col > 1 && side == :right) || (col < Nx && side == :left)
            ci  = commonindex(A[row, col], A[row, prev_col])
            msi = multiply_side_ident(AA, ci, EnvIdent.I[row])
            AA  = AA * msi
        end
        AA *= EnvIP.InProgress[row, opcode]
        if row > 1
            AAinds = inds(AA)
            Sbinds = inds(Sbelow[opcode][row-1])
            for aI in AAinds, sI in Sbinds
                if hastags(aI, tags(sI)) && !hasindex(AAinds, sI)
                    Sbelow[opcode][row-1] = replaceindex!(Sbelow[opcode][row-1], sI, aI)
                end
            end
        end
        thisS = row >= 2 ? Sbelow[opcode][row-1] * AA : AA
        push!(Sbelow[opcode], thisS)
    end
    return Sbelow
end


text:
# Test using the Guerry's Moral statistics of France dataset
@testset "Guerry" begin

    guerry = sdataset("Guerry")
    W = polyneigh(guerry)
    Wbin = wtransform(W, :binary)

    # In this dataset there is no difference between Queen or Rook contiguity
    @testset "W contiguity" begin
        @test nobs(W) == 85
        @test mean(W) ≈ 4.941176 atol = 1e-5
        @test median(W) ≈ 5.0
        @test minimum(W) == 2
        @test maximum(W) == 8
        @test nislands(W) == 0
        @test islands(W) == []

        # Connectivity histogram
        @test_nowarn RecipesBase.apply_recipe(Dict{Symbol, Any}(), W)
    end

    @testset "Spatial lag" begin
        Wlitercy = W * guerry.Litercy
        Wlitercy2 = slag(W, guerry.Litercy)

        @test Wlitercy[1:5] ≈ [44.75; 54.5; 22.666667; 42.75; 39.0] atol = 1e-5
        @test Wlitercy[81:85] ≈ [31.75; 26.166667; 21.833333; 67.666667; 47.0] atol = 1e-5

        @test Wlitercy == Wlitercy
        @test typeof(Wlitercy) == Vector{Float64}
        @test typeof(Wlitercy2) == Vector{Float64}
    end

    @testset "Moran's I" begin
        mguerry = moran(guerry.Litercy, W, rng = StableRNG(1234567), permutations = 9999)

        @test score(mguerry) ≈ 0.717605 atol = 1e-5
        @test length(scoreperms(mguerry)) == 9999
        @test mean(mguerry) ≈ -0.012594 atol = 1e-5
        @test std(mguerry) ≈ 0.070790 atol = 1e-5
        @test zscore(mguerry) ≈ 10.315064 atol = 1e-5
        @test pvalue(mguerry) ≈ 0.0001 atol = 1e-5
        @test SpatialDependence.expected(mguerry) == - 1 / (length(guerry.Litercy) - 1)
        @test SpatialDependence.testname(mguerry) == "Moran's I"

        @test_nowarn show(IOBuffer(), mguerry)

        # Test Moran's I distribution
        @test_nowarn RecipesBase.apply_recipe(Dict{Symbol, Any}(), mguerry)
    end

    @testset "Geary's c" begin
        cguerry = geary(guerry.Litercy, W, rng = StableRNG(1234567), permutations = 9999)

        @test score(cguerry) ≈ 0.250202 atol = 1e-5
        @test length(scoreperms(cguerry)) == 9999
        @test mean(cguerry) ≈ 1.000411 atol = 1e-5
        @test std(cguerry) ≈ 0.072280 atol = 1e-5
        @test zscore(cguerry) ≈ -10.379150 atol = 1e-5
        @test pvalue(cguerry) ≈ 0.0001 atol = 1e-5
        @test SpatialDependence.expected(cguerry) == 1
        @test SpatialDependence.testname(cguerry) == "Geary's c"

        @test_nowarn show(IOBuffer(), cguerry)

        # Test Geary's c distribution
        @test_nowarn RecipesBase.apply_recipe(Dict{Symbol, Any}(), cguerry)
    end

    @testset "Polygon mean center" begin
        cx, cy = meancenter(guerry.geometry)
        @test cx[1] ≈ 836732.867469879565760 atol = 1e-5
        @test cy[1] ≈ 2125844.293172690551728 atol = 1e-5
        @test cx[50] ≈ 814647.356215213309042 atol = 1e-5
        @test cy[50] ≈ 2350156.217068645637482 atol = 1e-5
    end

    @testset "Polygon centroid" begin
        cx, cy = centroid(guerry.geometry)
        @test cx[1] ≈ 832852.278780025895685 atol = 1e-5
        @test cy[1] ≈ 2126600.575969800818712 atol = 1e-5
        @test cx[50] ≈ 815193.599669872317463 atol = 1e-5
        @test cy[50] ≈ 2349544.092820507008582 atol = 1e-5
    end

    @testset "Moran 5 KNN" begin
        cx, cy = centroid(guerry.geometry)
        Wk5 = knearneigh(cx, cy, k = 5)
        morank5 = moran(guerry.Litercy, Wk5)
        @test score(morank5) ≈ 0.678857 atol = 1e-5
    end

    @testset "Local Moran" begin
        lmguerry = localmoran(guerry.Litercy, W, rng = StableRNG(1234567), permutations = 9999)

        @test score(lmguerry)[1:5] ≈ [-0.039511; 0.599223; 1.416860; 0.081434; -0.013868] atol = 1e-5
        @test score(lmguerry)[81:85] ≈ [0.270916; 0.603623; 1.488529; 2.145242; 0.203191] atol = 1e-5
        @test size(scoreperms(lmguerry), 1) == 85
        @test size(scoreperms(lmguerry), 2) == 9999
        @test mean(lmguerry)[1:5] ≈ [9.1e-5; -0.006589; -0.025567; -0.002702; -0.029142] atol = 1e-5
        @test mean(lmguerry)[81:85] ≈ [-0.004464; -0.006655; -0.02424; -0.021923; -0.003153] atol = 1e-5
        @test std(lmguerry)[1:5] ≈ [0.060504; 0.270397; 0.588832; 0.193349; 0.96366] atol = 1e-5
        @test std(lmguerry)[81:85] ≈ [0.315696; 0.320665; 0.58715; 0.514386; 0.197625] atol = 1e-5
        @test zscore(lmguerry)[1:5] ≈ [-0.654535; 2.240454; 2.44964; 0.435151; 0.01585] atol = 1e-5
        @test zscore(lmguerry)[81:85] ≈ [0.872291; 1.903161; 2.576463; 4.213106; 1.044118] atol = 1e-5
        @test pvalue(lmguerry)[1:5] ≈ [0.2536; 0.0148; 0.0044; 0.3323; 0.4805] atol = 1e-4
        @test pvalue(lmguerry)[81:85] ≈ [0.2015; 0.0239; 0.0021; 0.0001; 0.1532] atol = 1e-4

        @test assignments(lmguerry)[1:5] == [:LH, :HH, :LL, :HH, :HL]
        @test assignments(lmguerry)[81:85] == [:LL, :LL, :LL, :HH, :HH]
        @test count(issignificant(lmguerry, 0.05, adjust = :none)) == 38
        @test count(issignificant(lmguerry, 0.01, adjust = :none)) == 21
        @test count(issignificant(lmguerry, 0.05, adjust = :bonferroni)) == 6
        @test count(issignificant(lmguerry, 0.05, adjust = :fdr)) == 23
        @test assignments(lmguerry, 0.05, adjust = :none)[1:5] == [:ns, :HH, :LL, :ns, :ns]
        @test assignments(lmguerry, 0.05, adjust = :none)[81:85] == [:ns, :LL, :LL, :HH, :ns]

        @test SpatialDependence.testname(lmguerry) == "Local Moran"
        @test_nowarn show(IOBuffer(), lmguerry)

        # Test LISA Cluster Map
        @test_nowarn RecipesBase.apply_recipe(Dict{Symbol, Any}(), guerry, lmguerry)
        @test_nowarn RecipesBase.apply_recipe(Dict{Symbol, Any}(), guerry.geometry, lmguerry)

        # Test dividing by n instead of (n - 1)
        lmguerry = localmoran(guerry.Litercy, W, rng = StableRNG(1234567), permutations = 0, corrected = false)
        @test score(lmguerry)[1:5] ≈ [-0.039981; 0.606357; 1.433727; 0.082403; -0.014033] atol = 1e-5
        @test score(lmguerry)[81:85] ≈ [0.274141; 0.610809; 1.506250; 2.170780; 0.205610] atol = 1e-5
    end

    @testset "Local Geary" begin
        lcguerry = localgeary(guerry.Litercy, W, rng = StableRNG(1234567), permutations = 9999)

        @test score(lcguerry)[1:5] ≈ [1.191784; 0.275807; 0.439757; 0.949973; 3.134225] atol = 1e-5
        @test score(lcguerry)[81:85] ≈ [0.272243; 0.325157; 0.446885; 0.20288; 0.734317] atol = 1e-5
        @test size(scoreperms(lcguerry), 1) == 85
        @test size(scoreperms(lcguerry), 2) == 9999
        @test mean(lcguerry)[1:5] ≈ [1.013359; 1.468334; 3.269937; 1.157147; 3.95656] atol = 1e-5
        @test mean(lcguerry)[81:85] ≈ [1.412199; 1.664964; 3.272061; 2.744478; 1.209385] atol = 1e-5
        @test std(lcguerry)[1:5] ≈ [0.564736; 0.55419; 1.388559; 0.522587; 1.82078] atol = 1e-5
        @test std(lcguerry)[81:85] ≈ [0.950857; 0.883169; 1.389196; 0.968423; 0.490151] atol = 1e-5
        @test zscore(lcguerry)[1:5] ≈ [0.315944; -2.151836; -2.038214; -0.396441; -0.451638] atol = 1e-5
        @test zscore(lcguerry)[81:85] ≈ [-1.198872; -1.517045; -2.033677; -2.624471; -0.969228] atol = 1e-5
        @test pvalue(lcguerry)[1:5] ≈ [0.3506; 0.0036; 0.003; 0.373; 0.3455] atol = 1e-4
        @test pvalue(lcguerry)[81:85] ≈ [0.0722; 0.0325; 0.0031; 0.0002; 0.1712] atol = 1e-4

        @test assignments(lcguerry)[1:5] == [:N, :P, :P, :P, :P]
        @test assignments(lcguerry)[81:85] == [:P, :P, :P, :P, :P]
        @test count(issignificant(lcguerry, 0.05, adjust = :none)) == 49
        @test count(issignificant(lcguerry, 0.01, adjust = :none)) == 29
        @test count(issignificant(lcguerry, 0.05, adjust = :bonferroni)) == 9
        @test count(issignificant(lcguerry, 0.05, adjust = :fdr)) == 40
        @test assignments(lcguerry, 0.05, adjust = :none)[1:5] == [:ns, :P, :P, :ns, :ns]
        @test assignments(lcguerry, 0.05, adjust = :none)[81:85] == [:ns, :P, :P, :P, :ns]

        @test SpatialDependence.testname(lcguerry) == "Local Geary"
        @test_nowarn show(IOBuffer(), lcguerry)

        # Test LISA Cluster Map
        @test_nowarn RecipesBase.apply_recipe(Dict{Symbol, Any}(), guerry, lcguerry)
        @test_nowarn RecipesBase.apply_recipe(Dict{Symbol, Any}(), guerry.geometry, lcguerry)

        # Test :moran categories
        lcguerry = localgeary(guerry.Litercy, W, rng = StableRNG(1234567), permutations = 9999, categories = :moran)
        @test assignments(lcguerry)[1:5] == [:NE, :HH, :LL, :HH, :OP]
        @test assignments(lcguerry)[81:85] == [:LL, :LL, :LL, :HH, :HH]
        @test_nowarn RecipesBase.apply_recipe(Dict{Symbol, Any}(), guerry, lcguerry)

        # Test :moran categories with the Donatn's variables to test the :NE label when plotting maps
        lcguerry = localgeary(guerry.Donatns, W, rng = StableRNG(1234567), permutations = 9999, categories = :moran)
        @test SpatialDependence.labelsnames(lcguerry) == ["Not Significant"; "High-High"; "Low-Low"; "Other Positive"; "Negative"]
        @test_nowarn RecipesBase.apply_recipe(Dict{Symbol, Any}(), guerry, lcguerry)

        # Test dividing by n instead of (n - 1)
        lcguerry = localgeary(guerry.Litercy, W, rng = StableRNG(1234567), permutations = 0, corrected = false)
        @test score(lcguerry)[1:5] ≈ [1.205972; 0.279091; 0.444992; 0.961282; 3.171538] atol = 1e-5
        @test score(lcguerry)[81:85] ≈ [0.275484; 0.329028; 0.452205; 0.205295; 0.743059] atol = 1e-5
    end

    @testset "Getis-Ord" begin
        goguerry = getisord(guerry.Litercy, W, rng = StableRNG(1234567), permutations = 9999, star = false)

        @test score(goguerry)[1:5] ≈ [0.013602; 0.016636; 0.00684; 0.01303; 0.011971] atol = 1e-5
        @test score(goguerry)[81:85] ≈ [0.009624; 0.007924; 0.006588; 0.020725; 0.014329] atol = 1e-5
        @test size(scoreperms(goguerry), 1) == 85
        @test size(scoreperms(goguerry), 2) == 9999
        @test mean(goguerry)[1:5] ≈ [0.011893; 0.011896; 0.011901; 0.011893; 0.011923] atol = 1e-5
        @test mean(goguerry)[81:85] ≈ [0.011901; 0.011897; 0.011896; 0.011899; 0.011896] atol = 1e-5
        @test std(goguerry)[1:5] ≈ [0.002611; 0.002116; 0.002066; 0.002612; 0.003011] atol = 1e-5
        @test std(goguerry)[81:85] ≈ [0.002611; 0.002087; 0.00206; 0.002095; 0.00233] atol = 1e-5
        @test zscore(goguerry)[1:5] ≈ [0.654535; 2.240454; -2.44964; 0.435151; 0.01585] atol = 1e-5
        @test zscore(goguerry)[81:85] ≈ [-0.872291; -1.903161; -2.576463; 4.213106; 1.044118] atol = 1e-5
        @test pvalue(goguerry)[1:5] ≈ [0.2628; 0.0147; 0.0042; 0.3323; 0.4833] atol = 1e-4
        @test pvalue(goguerry)[81:85] ≈ [0.1931; 0.0224; 0.0018; 0.0001; 0.157] atol = 1e-4

        @test assignments(goguerry)[1:5] == [:H, :H, :L, :H, :H]
        @test assignments(goguerry)[81:85] == [:L, :L, :L, :H, :H]
        @test count(issignificant(goguerry, 0.05, adjust = :none)) == 38
        @test count(issignificant(goguerry, 0.01, adjust = :none)) == 21
        @test count(issignificant(goguerry, 0.05, adjust = :bonferroni)) == 6
        @test count(issignificant(goguerry, 0.05, adjust = :fdr)) == 23
        @test assignments(goguerry, 0.05, adjust = :none)[1:5] == [:ns, :H, :L, :ns, :ns]
        @test assignments(goguerry, 0.05, adjust = :none)[81:85] == [:ns, :L, :L, :H, :ns]

        @test SpatialDependence.testname(goguerry) == "Getis-Ord Gi"
        @test_nowarn show(IOBuffer(), goguerry)

        # Test LISA Cluster Map
        @test_nowarn RecipesBase.apply_recipe(Dict{Symbol, Any}(), guerry, goguerry)

        # Getis-Ord Gi*
        goguerry = getisord(guerry.Litercy, W, rng = StableRNG(1234567), permutations = 9999, star = true)

        @test score(goguerry)[1:5] ≈ [0.012985; 0.016231; 0.006398; 0.013045; 0.013977] atol = 1e-5
        @test score(goguerry)[81:85] ≈ [0.009318; 0.007815; 0.006183; 0.020095; 0.014127] atol = 1e-5

        @test SpatialDependence.testname(goguerry) == "Getis-Ord Gi*"

        # Getis-Ord Gi* with binary weights
        goguerry = getisord(guerry.Litercy, Wbin, rng = StableRNG(1234567), permutations = 9999, star = true)

        @test score(goguerry)[1:5] ≈ [0.064923; 0.113616; 0.044785; 0.065224; 0.055906] atol = 1e-5
        @test score(goguerry)[81:85] ≈ [0.046589; 0.054704; 0.043282; 0.140667; 0.084761] atol = 1e-5
    end

end


text:
export mnist,
    sampled_mnist,
    twenty_datasets,
    twenty_dataset_names

using MLDatasets
using CSV
using Pkg.Artifacts
using DataFrames

####################
# Constants
#####################

const twenty_dataset_names = [
        "accidents", "ad", "baudio", "bbc", "bnetflix", "book", "c20ng", "cr52", "cwebkb",
        "dna", "jester", "kdd", "kosarek", "msnbc", "msweb", "nltcs", "plants", "pumsb_star", "tmovie", "tretail",
        "binarized_mnist"
];

#####################
# Data loaders
#####################

function mnist(labeled = false)
    # transposing makes slicing by variable much much faster
    # need to take a copy to physically move the data around
    train_x = collect(Float32, transpose(reshape(MNIST.traintensor(), 28*28, :)))
    test_x  = collect(Float32, transpose(reshape(MNIST.testtensor(), 28*28, :)))

    train = DataFrame(train_x)
    valid = nothing # why is there no validation set in `MLDataSets`??
    test = DataFrame(test_x)
    if (labeled)
        train_y::Vector{UInt8} = MNIST.trainlabels()
        test_y::Vector{UInt8}  = MNIST.testlabels()
        train.y = train_y
        test.y = test_y
    end
    return train, valid, test
end

sampled_mnist() = twenty_datasets("binarized_mnist")

function twenty_datasets(name)
    @assert in(name, twenty_dataset_names)
    data_dir = artifact"density_estimation_datasets"
    function load(type)
        dataframe = CSV.read(data_dir*"/Density-Estimation-Datasets-1.0.1/datasets/$name/$name.$type.data", DataFrame;
            header=false, truestrings=["1"], falsestrings=["0"], type=Bool, strict=true)
                 # make sure the data is backed by a `BitArray`
        DataFrame((BitArray(Base.convert(Matrix{Bool}, dataframe))))
    end
    train = load("train")
    valid = load("valid")
    test = load("test")
    return train, valid, test
end


text:
using Fontconfig
using Test

# select basic font
pattern = Fontconfig.Pattern(family="Helvetica", size=10, hinting=true)

# check that methods don't crash
match(pattern)
format(pattern)
list(pattern)


text:
using StaticArrays

function emptytest(lib::PolyhedraLibrary)
    vr = vrep(SVector{2, Float64}[])
    p = polyhedron(vr, lib)
    generator_fulltest(p, vr)
end


text:
"""
    create_saturated_ols_approximation(dd::DataFrame, y::Symbol, x_variables::Array{Symbol,1}, degree::Int; intercept::Bool = true,  bases::Dict{Symbol,Float64} = Dict{Symbol,Float64}(x_variables .=> repeat([0.0],length(x_variables))))
This creates MultivariateFunction from an OLS regression predicting some variable. You input a dataframe and specify what column in that dataframe
is to be predicted by inputting a symbol y. you also put in an array of what x_variables should be used in prediction. A saturated ols model is then calculated up to the specified degree which is returned as a MultivariateFunction.
"""
function create_saturated_ols_approximation(dd::DataFrame, y::Symbol, x_variables::Array{Symbol,1}, degree::Int; intercept::Bool = true,  bases::Dict{Symbol,Float64} = Dict{Symbol,Float64}(x_variables .=> repeat([0.0],length(x_variables))))
    model = Array{PE_Function,1}()
    if intercept
        append!(model, [PE_Function(1.0,Dict{Symbol,PE_Unit}())] )
    end
    if degree > 0
        number_of_variables = length(x_variables)
        linear_set = Array{PE_Function,1}(undef, number_of_variables)
        for i in 1:length(x_variables)
            linear_set[i] = PE_Function(1.0,Dict{Symbol,PE_Unit}(x_variables[i] => PE_Unit(0.0,bases[x_variables[i]],1) ))
        end
        higher_order_terms = Array{Array{PE_Function,1},1}(undef,degree)
        higher_order_terms[1] = linear_set
        for i in 2:degree
            degree_terms = Array{PE_Function,1}()
            for j in 1:number_of_variables
                append!(degree_terms, linear_set[j] .* hcat(higher_order_terms[i-1]))
            end
            higher_order_terms[i] = degree_terms
        end
        append!(model, vcat(higher_order_terms...))
    end
    sum_of_functions = Sum_Of_Functions(model) # We put it through here to remove duplicates.
    return create_ols_approximation(dd, y, sum_of_functions)
end

"""
    create_ols_approximation(dd::DataFrame, y::Symbol, model::MultivariateFunction; allowrankdeficient = true)
    create_ols_approximation(dd::DataFrame, y::Symbol, model::Sum_Of_Functions; allowrankdeficient = true)
    create_ols_approximation(dd::DataFrame, y::Symbol, model::Sum_Of_Piecewise_Functions; allowrankdeficient = true)
This creates MultivariationFunction from an OLS regression predicting some variable. You input a dataframe and specify what column in that dataframe
is to be predicted by inputting a symbol y. You also input the regression model. This is input as a Array{MultivariateFunction,1}.
Each function that is input will be multiplied by the ols coefficient and will return a new function with these coefficients
incorporated.
"""
function create_ols_approximation(dd::DataFrame, y::Symbol, model::Array; allowrankdeficient = true)
    X = hcat(evaluate.(model, Ref(dd))...)
    y = dd[y]
    reg = fit(LinearModel, X,y, allowrankdeficient)
    coefficients = reg.pp.beta0
    updated_model = Sum_Of_Piecewise_Functions(model .* coefficients)
    return updated_model, reg
end
function create_ols_approximation(dd::DataFrame, y::Symbol, model::Sum_Of_Functions)
    return create_ols_approximation(dd, y, model.functions_)
end
function create_ols_approximation(dd::DataFrame, y::Symbol, model::Sum_Of_Piecewise_Functions)
    return create_ols_approximation(dd, y, vcat( model.functions_, model.global_funcs_ ))
end
"""
    create_ols_approximation(y::Array{Float64,1}, x::Array{Float64,1}, degree::Int; intercept::Bool = true, dim_name::Symbol = default_symbol, base_x::Float64 = 0.0)
    create_ols_approximation(y::Array{Float64,1}, x::Array{Date,1}, degree::Int; intercept::Bool = true, dim_name::Symbol = default_symbol, base_date::Date = global_base_date)

This predicts a linear relationship between the y and x arrays and creates a MultivariateFunction containing the approximation function. The degree specifies how many higher
order terms of x should be used (for instance degree 2 implies x and x^2 are both used to predict y).
"""
function create_ols_approximation(y::Array{Float64,1}, x::Array{Float64,1}, degree::Int; intercept::Bool = true, dim_name::Symbol = default_symbol, base_x::Float64 = 0.0)
    dd = DataFrame()
    dd[dim_name] = x
    dd[:y]       = y
    base_dict = Dict{Symbol,Float64}(dim_name => base_x)
    return create_saturated_ols_approximation(dd, :y, [dim_name], degree; intercept = intercept, bases = base_dict)
end

function create_ols_approximation(y::Array{Float64,1}, x::Array{Date,1}, degree::Int; intercept::Bool = true, dim_name::Symbol = default_symbol, base_date::Date = global_base_date)
    return  create_ols_approximation(y, years_from_global_base.(x), degree; intercept = intercept, dim_name = dim_name, base_x = years_from_global_base(base_date))
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualMachineScaleSetManagedDiskParameters <: SwaggerModel
    storageAccountType::Any # spec type: Union{ Nothing, StorageAccountType } # spec name: storageAccountType
    diskEncryptionSet::Any # spec type: Union{ Nothing, DiskEncryptionSetParameters } # spec name: diskEncryptionSet

    function VirtualMachineScaleSetManagedDiskParameters(;storageAccountType=nothing, diskEncryptionSet=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetManagedDiskParameters, Symbol("storageAccountType"), storageAccountType)
        setfield!(o, Symbol("storageAccountType"), storageAccountType)
        validate_property(VirtualMachineScaleSetManagedDiskParameters, Symbol("diskEncryptionSet"), diskEncryptionSet)
        setfield!(o, Symbol("diskEncryptionSet"), diskEncryptionSet)
        o
    end
end # type VirtualMachineScaleSetManagedDiskParameters

const _property_map_VirtualMachineScaleSetManagedDiskParameters = Dict{Symbol,Symbol}(Symbol("storageAccountType")=>Symbol("storageAccountType"), Symbol("diskEncryptionSet")=>Symbol("diskEncryptionSet"))
const _property_types_VirtualMachineScaleSetManagedDiskParameters = Dict{Symbol,String}(Symbol("storageAccountType")=>"StorageAccountType", Symbol("diskEncryptionSet")=>"DiskEncryptionSetParameters")
Base.propertynames(::Type{ VirtualMachineScaleSetManagedDiskParameters }) = collect(keys(_property_map_VirtualMachineScaleSetManagedDiskParameters))
Swagger.property_type(::Type{ VirtualMachineScaleSetManagedDiskParameters }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualMachineScaleSetManagedDiskParameters[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetManagedDiskParameters }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetManagedDiskParameters[property_name]

function check_required(o::VirtualMachineScaleSetManagedDiskParameters)
    true
end

function validate_property(::Type{ VirtualMachineScaleSetManagedDiskParameters }, name::Symbol, val)
end


text:
# implicit scheme for advection-diffusion

using Revise
using AdFem
using ADCME
using PyPlot
using Statistics
using LinearAlgebra
using ADCMEKit

m = 100
n = 100
h = 1/n

function f1func(x,y)
    4.0*x^3*y^3 - 4.0*x^3*y*(1 - y^2) - 4.0*x*y^3*(1 - x^2) + 4.0*x*y*(1 - x^2)*(1 - y^2) - 2.0*x*(1 - x) - 4*y*(1 - y) - 1
end
function f2func(x,y)
    -20*x^2*y^2*(1 - x^2) - 10.0*x^2*y^2*(1 - y^2) + 4*x^2*(1 - x^2)*(1 - y^2) + 1.0*x*y - x*(1 - y) + 2.0*y^2*(1 - x^2)*(1 - y^2) - 1.0*y*(1 - x) + (1 - x)*(1 - y) - 1
end

function gfunc(x,y)
    -2*x^2*y^3*(1 - x^2) + 2*x^2*y*(1 - x^2)*(1 - y^2) - x*y*(1 - y) + y*(1 - x)*(1 - y)
end

hmat = zeros(4*m*n, 3, 3)
for i = 1:4*m*n
    hmat[i,:,:] = 2diagm(0=>ones(3))
end
hmat = constant(hmat)

K = compute_fem_stiffness_matrix(hmat, m, n, h)
B = constant(compute_interaction_matrix(m, n, h))

Z = [K -B'
B spzero(size(B,1))]

bd = bcnode("all", m, n, h)
Z, _ = fem_impose_Dirichlet_boundary_condition(Z, bd, m, n, h)

F1 = eval_f_on_gauss_pts(f1func, m, n, h)
F2 = eval_f_on_gauss_pts(f2func, m, n, h)
F = compute_fem_source_term(F1, F2, m, n, h)
xy = fvm_nodes(m, n, h)
G = gfunc.(xy[:,1], xy[:,2])
G = compute_fvm_source_term(G, m, n, h)
G = G .- mean(G)
rhs = [-F;G]
rhs[[bd;bd.+(m+1)*(n+1)]] .= 0.0

sol = Z[1:end-1, 1:end-1]\rhs[1:end-1]

sess = Session(); init(sess)
S = run(sess, sol)

# visualize_displacement(S[1:2(m+1)*(n+1)],m, n,h)
xy = fem_nodes(m, n, h)
U = @. (1-xy[:,1])*xy[:,1]*(1-xy[:,2])*xy[:,2]
figure(figsize=(12,5))
subplot(121)
visualize_scalar_on_fem_points(U, m, n, h)
title("Reference")
subplot(122)
visualize_scalar_on_fem_points(S[1:(m+1)*(n+1)], m, n, h)
title("Computed")
savefig("stokes1.png")

U = @. (1-xy[:,1]^2)*xy[:,1]^2*(1-xy[:,2]^2)*xy[:,2]^2
figure(figsize=(12,5))
subplot(121)
visualize_scalar_on_fem_points(U, m, n, h)
title("Reference")
subplot(122)
visualize_scalar_on_fem_points(S[(m+1)*(n+1)+1:2(m+1)*(n+1)], m, n, h)
title("Computed")
savefig("stokes2.png")

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
## Input/Output functions for networks
# import
import Base.print

# Printing function
function print(g::Network)
  n = length(g.tails) - 1
  m = 0
  str = string()
  for (id,start) in enumerate(g.tails[1:n])
    stop = g.tails[id+1] - 1
    str *= "node $id :"
    for a in g.links[start:stop]
      if a.cap > 0
        aux = (a.head,a.cap)
        str *= " $aux"
        m += 1
      end
    end
    str *= "\n"
  end
  str2 = "The network has $n nodes and $m links."
  str2 *= "The source node is $(g.source) and the sink node is $(g.target).\n"
  str = str2 * str
  print(str)
end

function printFull(g::Network)
  n = length(g.tails) - 1
  m = 0
  str = string()
  for (id,start) in enumerate(g.tails[1:n])
    stop = g.tails[id+1] - 1
    str *= "node $id :"
    for a in g.links[start:stop]
      aux = (a.head,a.cap)
      str *= " $aux"
      m += 1
    end
    str *= "\n"
  end
  str2 = "The network has $n nodes and $m links."
  str2 *= "The source node is $(g.source) and the sink node is $(g.target).\n"
  str = str2 * str
  print(str)
end

function printflow(g::Network, flows::Vector{Float64})
  for (id,e) in enumerate(g.links)
    esym = g.links[e.sym]
    println("e",id,": ",esym.head,"->",e.head,"=",flows[id],"/",e.cap)
  end
end

# Network input in natural graph description
function Network(links::Vector{Link}, directed::Bool, s::Int, t::Int)
  if directed
    links, n = addReverseLinks(links)
  else
    links, n = orientLinks(links)
  end
  sort!(links)
  l = length(links)
  tails = zeros(Int,n+1)
  tails[n+1] = l + 1
  arcs = Vector{Arc}()
  sizehint!(arcs, l)

  function findSymmetric(a::Link)
    start = tails[a[2]]
    stop = tails[a[2]+1] - 1
    for (id,e) in enumerate(arcs[start:stop])
      if e.head == a[1]
        return start + id - 1
      end
    end
    return 0
  end

  tail = 0
  for (id,a) in enumerate(links)
    if tail < a[1]
      for i in tail+1:a[1]
        tails[i] = id
      end
      tail = a[1]
    end
    if a[1] < a[2]
      push!(arcs,Arc(0,a[2],a[3]))
    else
      rev = findSymmetric(a)
      asym = arcs[rev]
      asym.sym = id
      push!(arcs,Arc(rev,a[2],a[3]))
    end
  end

  Network(s,t,tails,arcs)
end

function Network(links::Vector{SimpleLink}, directed::Bool, s::Int, t::Int)
  Network(map(weightLink, links), directed, s, t)
end

# DIMACS format input
function orderedDIMACS(fname::AbstractString, foutput = "")
  # Init Strings
  comments = Vector{AbstractString}()
  links = Vector{Link}()
  linklines = Vector{Int}()
  problem_line = 0
  source_line = 0
  sink_line = 0
  problem = ""
  source = ""
  sink = ""

  open(fname,"r") do file
    for (id,line) in enumerate(eachline(file))
      words = split(line)
      if words[1] == "p"
        problem = line
        problem_line = id
      elseif words[1] == "c"
        push!(comments,line)
      elseif words[1] == "n"
        if words[3] == "s"
          source = line
          source_line = id
        else words[3] == "t"
          sink = line
          sink_line = id
        end
      else # words[1] == 'a'
        tail = parse(Int,words[2])
        head = parse(Int,words[3])
        cap = parse(Float64,words[4])
        push!(links,(tail,head,cap))
        push!(linklines,id)
      end
    end
  end
  sort!(links)

  strout = ""
  if foutput == ""
    foutput = fname
  end
  open(foutput,"w") do file
    linecpt = 0
    while !isempty(links) || !isempty(comments)
      linecpt += 1
      if linecpt == problem_line
        strout *= problem
      elseif linecpt == source_line
        strout *= source
      elseif linecpt == sink_line
        strout *= sink
      elseif !isempty(links) && (linecpt == linklines[1])
        l = shift!(links)
        strout *= "a $(l[1]) $(l[2]) $(l[3])\n"
        shift!(linklines)
      else
        strout *= shift!(comments)
      end
    end
    write(file,strout)
  end
end

function dimacsNetwork(fname::AbstractString)
  links = Vector{Link}()
  source = 0
  sink = 0
  directed = true

  open(fname,"r") do file
    for line in eachline(file)
      words = split(line)
      if words[1] == "n"
        if words[3] == "s"
          source = parse(Int,words[2])
        else words[3] == "t"
          sink = parse(Int,words[2])
        end
      elseif words[1] == "a"
        tail = parse(Int,words[2])
        head = parse(Int,words[3])
        cap = parse(Float64,words[4])
        push!(links,(tail,head,cap))
      end
    end
  end

  Network(links, directed, source, sink)
end

function csvNetwork(fname::AbstractString, source::Int, sink::Int)
  links = Vector{Link}()
  directed = true

  data = readcsv(fname)
  for j in 1:size(data,1)
    push!(links, (round(Int,data[j,1]), round(Int,data[j,2]), data[j,3]))
  end
  Network(links,directed,source,sink)
end

function Network(fname::AbstractString,format::Symbol, st = (1,2) )
  if format == :DIMACS
    dimacsNetwork(fname)
  elseif format == :CSV
    source, sink = st
    csvNetwork(fname,source, sink)
  end
end


text:
@enum generator_inner_vars begin
    τe_var = 1
    τm_var = 2
    Vf_var = 3
    V_pss_var = 4
    VR_gen_var = 5
    VI_gen_var = 6
    ψd_var = 7
    ψq_var = 8
end

Base.to_index(ix::generator_inner_vars) = Int(ix)
include("generator_components/machine.jl")
include("generator_components/shaft.jl")
include("generator_components/pss.jl")
include("generator_components/avr.jl")
include("generator_components/turbinegov.jl")

"""
"""
mutable struct DynGenerator{M <: Machine,
                            S <: Shaft,
                            A <: AVR,
                            TG<: TurbineGov,
                            P <: PSS} <: DynInjection
    number::Int64
    name::Symbol
    bus::PSY.Bus
    ω_ref::Float64
    V_ref::Float64
    P_ref::Float64
    machine::M
    shaft::S
    avr::A
    tg::TG
    pss::P
    n_states::Int64
    states::Vector{Symbol}
    inner_vars::Vector{Float64}
    local_state_ix::Dict{GeneratorComponent,Vector{Int64}}
    input_port_mapping::Dict{GeneratorComponent,Vector{Int64}}

        function DynGenerator(number::Int64,
                              name::Symbol,
                              bus::PSY.Bus,
                              ω_ref::Float64,
                              V_ref::Float64,
                              P_ref::Float64,
                              machine::M,
                              shaft::S,
                              avr::A,
                              tg::TG,
                              pss::P) where {M <: Machine,
                                             S <: Shaft,
                                             A <: AVR,
                                             TG <: TurbineGov,
                                             P <: PSS}

            n_states = (machine.n_states +
                       shaft.n_states +
                       avr.n_states +
                       tg.n_states +
                       pss.n_states)

            states = vcat(machine.states,
                          shaft.states,
                          avr.states,
                          tg.states,
                          pss.states)

            components = [machine, shaft, avr, tg, pss]

            new{M,S,A,TG,P}(number,
                        name,
                        bus,
                        ω_ref,
                        V_ref,
                        P_ref,
                        machine,
                        shaft,
                        avr,
                        tg,
                        pss,
                        n_states,
                        states,
                        zeros(8),
                        _make_state_mapping(components, states),
                        _make_port_mapping(components, states))

        end
  end


text:
using DiffEqFlux, Test

u0 = Float32[2.0; 0.0]
dudt2 = FastChain((x, p) -> x.^3,
                  FastDense(2, 50, tanh),
                  FastDense(50, 2))
p = initial_params(dudt2)
@test !DiffEqSensitivity.hasbranching(dudt2,u0,p)

dudt2 = FastChain((x, p) -> x.^3,
                  StaticDense(2, 4, tanh),
                  StaticDense(4, 2))
p = initial_params(dudt2)
@test !DiffEqSensitivity.hasbranching(dudt2,u0,p)


text:
# project x1 on to x2
project(x1, x2) = (x1' * x2) / (x2' * x2) * x2

# compute intercept and slope from two points
function get_coefs(x1, x2)
     β1 = (x2[2] - x1[2]) / (x2[1] - x1[1])
     return x2[2] - β1 * x2[1], β1
end

# make a line
make_line(β0, β1, x) = β0 .+ β1 .* x


text:
isinstalled(pkg::AbstractString) =
  pkg != "METADATA" && pkg != "REQUIRE" && pkg[1] != '.' && Pkg.cd(isdir, pkg)

asl_pkgs = ["BulletCollision", "AstrobeeRobot", "PandaRobot"]
for pkg in asl_pkgs
  if isinstalled(pkg)
    println("$pkg already installed")
    continue
  else
    println("Cloning out $pkg")
  end
  Pkg.clone("git://github.com/StanfordASL/$pkg.jl.git")
  Pkg.build(pkg)
end


text:
# Calculate the Sum of Distance between the Cities

function sumLen(city)

    n = size(city)[1]
    len = 0

    for ii = 1: (n - 1)

        len += euclidean(city[ii, :], city[ii + 1, :])

    end

    len += euclidean(city[1, :], city[n, :])

    return len

end

max_stars_repo_name: mizu-bai/Try-Julia

max_issues_repo_name: mizu-bai/Data-Science-Demos

max_forks_repo_name: mizu-bai/Data-Science-Demos


text:
Pkg.init()
Pkg.update()
Pkg.add("HTTP")
Pkg.resolve()

max_stars_repo_name: bcbi/PhytoKB-API

max_issues_repo_name: bcbi/PhytoKB-API

max_forks_repo_name: bcbi/PhytoKB-API


text:
n = 10000;
a = rand(Float32,(n,n))
b = rand(Float32,(n,n))
c = zeros(Float32,(n,n));
t1 = time();

for i = 1:n
  for j = 1:n
    c[i,j] = a[i,j] + b[i,j]
  end
end

t2 = time();

for i = 1:10
c = a+ b;
end

t3 = time();

dt1 = t2 - t1
dt2 = t3 - t2
n2 = n*n

mflops1 = n2 / dt1 / 1e6;
mflops2 = n2 / dt2 / 1e5;

println("times $dt1,$dt2, mflops $mflops1,$mflops2")

max_stars_repo_name: oeclint/BIDMach

max_issues_repo_name: oeclint/BIDMach

max_forks_repo_name: oeclint/BIDMach


text: activate_github("QuantEcon/QuantEconLectureAllPackages", tag = "v0.9.5", add_default_environment = true)

max_stars_repo_name: aderdzyan/QuantEcon-lecture-source-jl

max_issues_repo_name: aderdzyan/QuantEcon-lecture-source-jl

max_forks_repo_name: aderdzyan/QuantEcon-lecture-source-jl


text:
# segment types
"""
Finite set {0, .., k-1}
"""
struct Segment{T} <: AbstractSegment
    itr :: UnitRange{T}
end

"""
Finite set {c, .., c + k-1}
"""
struct ShiftedSegment{T} <: AbstractSegment
    itr :: UnitRange{T}
end

# Segment initializers

"""
Construct a Segment [0..k-1]
"""
Segment(k :: Int) = Segment(0 : (k-1))

"""
Construct a Segment [start..finish]
"""
Segment(start, finish) = ShiftedSegment(start : finish)

# INTERFACE FUNCTIONS
iterator(s :: Segment) = s.itr

max_stars_repo_name: kirtsar/DiscreteFunctions

max_issues_repo_name: kirtsar/DiscreteFunctions


text:
include("Particles.jl")
import Base: length, size

"Abstract struct for all the type of packings to be used."
abstract type AbstractPacking{d, T} end

#########################################################################################################
#########################################################################################################
### Functions for generic (i.e. any type of packings) packings
#########################################################################################################
#########################################################################################################

get_positions(packing::AbstractPacking) = getfield.(packing.Particles, :X)
length(packing::AbstractPacking) = length(packing.Particles)
size(packing::AbstractPacking{d, T}) where {d, T<:Real} = (d, length(packing))

"Obtain the set of stable particles in a packing."
function get_non_rattlers(packing::AbstractPacking)
    N = length(packing)
    return (1:N)[is_stable_particle.(packing.Particles)]
end

"Obtain the set of rattlers (i.e. unstable particles) in a packing."
function get_rattlers(packing::AbstractPacking)
    N = length(packing)
    return (1:N)[is_a_rattler.(packing.Particles)]
end

function get_coordination_number(packing::AbstractPacking; only_stable::Bool=false)
    if only_stable # so only the contacts of non-rattlers are considered
        non_rattlers = get_non_rattlers(packing)
        return get_coordination_number.(packing.Particles[non_rattlers])
    else
        return get_coordination_number.(packing.Particles)
    end
end

"Test whether a given packing is isostatic or not. The output is a boolean"
function is_isostatic(packing::AbstractPacking{d, T}; only_stable::Bool=true)::Tuple{Bool, Int64, Int64} where {d, T<:Real}
    Nnr = length(get_non_rattlers(packing)); # amount of non-rattlers (i.e. stable particles)
    N_dof = d*(Nnr-1)+1 # number of degrees of freedom
    Nc = 0.5*sum(get_coordination_number(packing; only_stable=only_stable)) # number of contacts; if 'only_stable=true' only non-rattlers are considered
    return Nc==N_dof, Nc, Nnr
end

"Output an array of `StaticVector`s corresponding to the sum of forces acting on each particle of the packing."
function total_force(packing::AbstractPacking{d, T})::Vector{SVector{d, T}} where {d, T<:Real}
    total_force.(packing.Particles)
end

distances_between_centers(packing::AbstractPacking) = distances_between_centers(get_positions(packing))

function difference_in_packings(pack1::AbstractPacking, pack2::AbstractPacking)
    Xs1 = get_positions(pack1); Xs2 = get_positions(pack2)
    isostatic, Nc, Nnr = is_isostatic(pack1)
    if isostatic
        non_rattlers = get_non_rattlers(pack1)
        norms = norm.(Xs1-Xs2)[non_rattlers]
    else
        norms = norm.(Xs1-Xs2)
    end
    return mean(norms), std(norms)
end

#########################################################################################################
#########################################################################################################
### Monodisperse packing type and related constructors
#########################################################################################################
#########################################################################################################

#When  dealing with monodisperse packings it's more efficient to store the particles' radius as a field of
#the packing itself rather than of each particle. Of course, in case of polydisperse configurations, that would change.
"""
    MonoPacking{d,T}

Structure used to store all the relevant info of a *monodisperse* hard-spheres packing in d-dimensions.

Its fields are:
1. 'Particles': An vector of `MonoParticle` that contains all the particles that make the packing.
2. 'R': The radius of *all* particles.
3. 'mechanical_equilibrium': A boolean that specifies whether force balances is satisfied
all across the packing (i.e. for each particle).
4. 'jammed': A boolean that specifies whether a packing is jammed or not.

Note that a packing might not be necessarily jammed nor in mechanical_equilibrium. That is
why specifying such fields is important.

See also: [`MonoParticle`](@ref), [`PolyPacking`](@ref).
"""
mutable struct MonoPacking{d,T} <: AbstractPacking{d, T}
    Particles::Vector{MonoParticle{d, T}} # Vector of particles, each with info of position, contacts, forces, and indices of contacts
    R::T  # The radius of all the particles in a packing
    mechanical_equilibrium::Bool  # boolean to indicate if the packing is in mechanical equilibrium (it should be, in general)
    jammed::Bool # Boolean to indicate if the packing corresponds to a jammed state or not
end
MonoPacking([MonoParticle(rand(3), rand(), rand(3,3), rand(3), collect(1:3)), MonoParticle(rand(3), rand(), rand(3,3), rand(3), collect(1:3))], rand(), false, false)

#### Format how a `MonoPacking` is shown in IO
function Base.show(io::IO, packing::MonoPacking{d,T}) where {d, T<:Real}
    N = length(packing)
    iso, Nc, Nnr = is_isostatic(packing)
    fr = (N-Nnr)/Nnr
    println(io, d, "d Monodisperse packing\t of N= ", N, " particles \t of radius R= ", packing.R)
    println(io, Nnr, " stable particles; \t fraction of rattlers = ", round(fr, digits=3))
    println(io, "jammed: ", packing.jammed, "\t isostatic: ", iso ,"\t and in mechanical equilibrium: ", packing.mechanical_equilibrium)
end

#####################
# Now we define some useful constructors
#####################

# Empty packing
"Create an empty d-dimensional packing of N particles"
MonoPacking(d::Int64, N::Int64=1) = MonoPacking(Vector{MonoParticle{d, Float64}}(undef, N), 0.0, false, false)
MonoPacking(5)

@doc raw"""
    MonoPacking(Xs::Vector{SVector{d, PeriodicNumber{T}}}, contact_vecs::Vector{Vector{SVector{d,T}}}, fs::Vector{Vector{T}}, neighbours::Vector{Vector{Int64}}, R::T, jammed::Bool=false; <keyword arguments>)
    MonoPacking(Xs::Vector{SVector{d, PeriodicNumber{T}}}, contact_vecs::Vector{Matrix{T}}, fs::Vector{Vector{T}}, neighbours::Vector{Vector{Int64}}, R::T, jammed::Bool=false; <keyword arguments>)

Create a d-dimensional packing, where the attributes of each particles are inferred from the
 elements of the input arrays.

The number of particles (N) is inferred from the length of the input arrays.
Then a Vector{MonoParticle{d,T}} of size N but undefined elements is constructed. Each of its
 elements is then defined by calling 'MonoParticle(Xs[i], contact_vecs[i], fs[i], neighbours[i])'.

The constructors also asses whether force balance for each particle is satisfied, within a
given precision 'tol_mechanical_equilibrium' (that defaults to `default_tol_force_equilibrium=1e-12`).
When this condition is not met, it throws a warning, but the packing is created.
"""
MonoPacking()

function MonoPacking(Xs::Vector{SVector{d, PeriodicNumber{T}}}, contact_vecs::Vector{Vector{SVector{d,T}}}, fs::Vector{Vector{T}}, neighbours::Vector{Vector{Int64}}, R::T, jammed::Bool=false; tol_mechanical_equilibrium::Float64=default_tol_force_equilibrium, verbose::Bool=true) where {d, T<:Real}
    N = length(Xs) # number of particles
    # th construct the packing all the input arrays should have the same length
    if N==length(fs) && N==length(contact_vecs) && N==length(neighbours)
        particles = Vector{MonoParticle{d,T}}(undef, N) # empty array of particles
        mech_eq = true # default value of mechanical equilibrium condition

        # create particles with the elements of the input arrays
        for i in 1:N
            particles[i] = MonoParticle(Xs[i], contact_vecs[i], fs[i], neighbours[i])
        end

    # test whether the mechanical equilibrium condition is satisfied for all particles, at least within a 'tol_mechanical_equilibrium' precision
        any(norm.(total_force.(particles)).>tol_mechanical_equilibrium) && (mech_eq=false)

        # if this is not the case throw a warning, but create the packing anyway
        if !mech_eq && verbose
            fmismatch = maximum(norm.(total_force.(particles)))
            @warn "Force balance condition is NOT satisfied! Max force mismatch = $fmismatch ;\tCreating the packing anyway."
        end

        return MonoPacking(particles, R, mech_eq, jammed) # return the packing created

    else # if there is a mismatch in the dimensions of the input array throw an error (and some info)
        error("Number of dimensions mismatch! There are: ", N, " particles \t ", length(contact_vecs) ," contact vectors \t", length(fs), " forces and \t ", length(neighbours), " list of contacts indices")
    end
end

function MonoPacking(Xs::Vector{SVector{d, PeriodicNumber{T}}}, contact_vecs::Vector{Matrix{T}}, fs::Vector{Vector{T}}, neighbours::Vector{Vector{Int64}}, R::T, jammed::Bool=false; tol_mechanical_equilibrium::Float64=default_tol_force_equilibrium, verbose::Bool=true) where {d, T<:Real}
    MonoPacking(Xs, svectors.(contact_vecs, Val(d)), fs, neighbours, R, jammed; tol_mechanical_equilibrium=tol_mechanical_equilibrium, verbose=verbose)
end

function MonoPacking(particles::Vector{MonoParticle{d, T}}, R::T, jammed::Bool=false;  tol_mechanical_equilibrium::Float64=default_tol_force_equilibrium, verbose::Bool=true) where {d, T<:Real}
    force_balance = force_equilibrium(particles; tol_mechanical_equilibrium=tol_mechanical_equilibrium)
    if !force_balance && verbose
        fmismatch = maximum(norm.(total_force.(particles)))
        @warn "Force balance condition is NOT satisfied! Max force mismatch = $fmismatch ; \tCreating the packing anyway."
    end

    MonoPacking(particles, R, force_balance, jammed)
end

"""
    distances_between_particles(packing::MonoPacking)

Compute the distance between the *borders* of all pairs of particles in 'packing'.

For a given pair, the distance between the centers of the two particles is computed, by calling
`distances_between_centers`, and then the diameter is subtracted. In this way, the result is
really the distance between the border of each pair of particles.

See also [`distances_between_centers`](@ref)
"""
function distances_between_particles(packing::MonoPacking)
    distances_between_centers(packing) .- (2*packing.R)
end

"""
    check_for_overlaps(packing::MonoPacking, tolerance::Float64)

Apply `check_for_overlaps` to all the particles in 'packing'.
"""
function check_for_overlaps(packing::MonoPacking{d, T}, tolerance::T) where {d, T<:AbstractFloat}
    R = packing.R
    Xs = get_positions(packing)
    check_for_overlaps(Xs, R, tolerance)
end

#################################################
# Other useful functions
################################################

"""
    packing_fraction(d::Int64, R::Real, N::Int64, L::Real=1.0)

Compute the packing fraction of N d-dimensional hyperspheres of the same radius, 'R', inside
 a box of size 'L'.

See also [`volume`](@ref), [`volume_d_ball`](@ref)
 """
function packing_fraction(d::Int64, R::Real, N::Int64, L::Real=1.0)::Float64
    return N*volume_d_ball(d, R)/L^d
end
packing_fraction(3, 0.5, 2, 1.0)

"""
    packing_fraction(packing::MonoPacking{d, T})

Compute the packing fraction of a monodisperse packing.

See also [`volume`](@ref), [`volume_d_ball`](@ref)
"""
function packing_fraction(packing::MonoPacking{d, T})::Float64 where {d, T}
    R = packing.R
    N = length(packing.Particles)
    L = packing.Particles[1].X[1].L
    packing_fraction(d, R, N, L)
end

#########################################################################################################
#########################################################################################################
### Polydisperse packing type and related constructors
#########################################################################################################
#########################################################################################################
"""
    PolyPacking{d,T}

Structure used to store all the relevant info of a *monodisperse* hard-spheres packing in d-dimensions.

Its fields are:
1. 'Particles': An vector of `Particle`s that contains all the particles that make the packing.
2. 'mechanical_equilibrium': A boolean that specifies whether force balances is satisfied
all across the packing (i.e. for each particle).
3. 'jammed': A boolean that specifies whether a packing is jammed or not.

Note that a packing might not be necessarily jammed nor in mechanical_equilibrium. That is
why specifying such fields is important.

See also: [`Particle`](@ref), [`MonoPacking`(@ref)].
"""
mutable struct PolyPacking{d,T} <: AbstractPacking{d, T}
    Particles::Vector{Particle{d, T}} # Vector of particles, each with info of position, contacts, forces, and indices of contacts
    mechanical_equilibrium::Bool  # boolean to indicate if the packing is in mechanical equilibrium (it should be, in general)
    jammed::Bool # Boolean to indicate if the packing corresponds to a jammed state or not
end
PolyPacking([Particle(rand(3), 0.1, rand(), rand(3,3), rand(3), collect(1:3)), Particle(rand(3), 0.1, rand(), rand(3,3), rand(3), collect(1:3))],  false, false)

function get_radii(packing::PolyPacking{d,T}) where {d, T<:AbstractFloat}
    getfield.(packing.Particles, :R)
end
get_radii(PolyPacking([Particle(rand(3), 0.1, rand(), rand(3,3), rand(3), collect(1:3)), Particle(rand(3), 0.5, rand(), rand(3,3), rand(3), collect(1:3))],  false, false))

#### Format how a `PolyPacking` is shown in IO
function Base.show(io::IO, packing::PolyPacking{d,T}) where {d, T<:Real}
    N = length(packing)
    iso, Nc, Nnr = is_isostatic(packing)
    fr = (N-Nnr)/Nnr
    println(io, d, "d Polydisperse packing\t of N= ", N, " particles \t of radii = ", get_radii(packing))
    println(io, Nnr, " stable particles; \t fraction of rattlers = ", round(fr, digits=3))
    println(io, "jammed: ", packing.jammed, "\t isostatic: ", iso ,"\t and in mechanical equilibrium: ", packing.mechanical_equilibrium)
end

#####################
# Now we define some useful constructors
#####################

# Empty packing
"Create an empty d-dimensional packing of N particles"
PolyPacking(d::Int64, N::Int64=1) = PolyPacking(Vector{Particle{d, Float64}}(undef, N), false, false)
PolyPacking(5)

@doc raw"""
    PolyPacking(Xs::Vector{SVector{d, PeriodicNumber{T}}}, contact_vecs::Vector{Vector{SVector{d,T}}}, fs::Vector{Vector{T}}, neighbours::Vector{Vector{Int64}}, jammed::Bool=false; <keyword arguments>)
    PolyPacking(Xs::Vector{SVector{d, PeriodicNumber{T}}}, contact_vecs::Vector{Matrix{T}}, fs::Vector{Vector{T}}, neighbours::Vector{Vector{Int64}}, jammed::Bool=false; <keyword arguments>)

Create a d-dimensional polydisperse packing, where the attributes of each particles are inferred from the
 elements of the input arrays.

The number of particles (N) is inferred from the length of the input arrays.
Then a Vector{Particle{d,T}} of size N but undefined elements is constructed. Each of its
 elements is then defined by calling 'Particle(Xs[i], contact_vecs[i], fs[i], neighbours[i])'.

The constructors also asses whether force balance for each particle is satisfied, within a
given precision 'tol_mechanical_equilibrium' (that defaults to `default_tol_force_equilibrium=1e-12`).
When this condition is not met, it throws a warning, but the packing is created.
"""
PolyPacking()

function PolyPacking(Xs::Vector{SVector{d, PeriodicNumber{T}}}, contact_vecs::Vector{Vector{SVector{d,T}}}, fs::Vector{Vector{T}}, neighbours::Vector{Vector{Int64}}, jammed::Bool=false; tol_mechanical_equilibrium::Float64=default_tol_force_equilibrium, verbose::Bool=true) where {d, T<:Real}
    N = length(Xs) # number of particles
    # th construct the packing all the input arrays should have the same length
    if N==length(fs) && N==length(contact_vecs) && N==length(neighbours)
        particles = Vector{Particle{d,T}}(undef, N) # empty array of particles
        mech_eq = true # default value of mechanical equilibrium condition

        # create particles with the elements of the input arrays
        for i in 1:N
            particles[i] = PolyParticle(Xs[i], contact_vecs[i], fs[i], neighbours[i])
        end

    # test whether the mechanical equilibrium condition is satisfied for all particles, at least within a 'tol_mechanical_equilibrium' precision
        any(norm.(total_force.(particles)).>tol_mechanical_equilibrium) && (mech_eq=false)

        # if this is not the case throw a warning, but create the packing anyway
        if !mech_eq && verbose
            fmismatch = maximum(norm.(total_force.(particles)))
            @warn "Force balance condition is NOT satisfied! Max force mismatch = $fmismatch ;\tCreating the packing anyway."
        end

        return PolyPacking(particles, mech_eq, jammed) # return the packing created

    else # if there is a mismatch in the dimensions of the input array throw an error (and some info)
        error("Number of dimensions mismatch! There are: ", N, " particles \t ", length(contact_vecs) ," contact vectors \t", length(fs), " forces and \t ", length(neighbours), " list of contacts indices")
    end
end

function PolyPacking(Xs::Vector{SVector{d, PeriodicNumber{T}}}, contact_vecs::Vector{Matrix{T}}, fs::Vector{Vector{T}}, neighbours::Vector{Vector{Int64}}, jammed::Bool=false; tol_mechanical_equilibrium::Float64=default_tol_force_equilibrium, verbose::Bool=true) where {d, T<:Real}
    PolyPacking(Xs, svectors.(contact_vecs, Val(d)), fs, neighbours, jammed; tol_mechanical_equilibrium=tol_mechanical_equilibrium, verbose=verbose)
end

function PolyPacking(particles::Vector{Particle{d, T}}, jammed::Bool=false;  tol_mechanical_equilibrium::Float64=default_tol_force_equilibrium, verbose::Bool=true) where {d, T<:Real}
    force_balance = force_equilibrium(particles; tol_mechanical_equilibrium=tol_mechanical_equilibrium)
    if !force_balance && verbose
        fmismatch = maximum(norm.(total_force.(particles)))
        @warn "Force balance condition is NOT satisfied! Max force mismatch = $fmismatch ; \tCreating the packing anyway."
    end

    PolyPacking(particles, force_balance, jammed)
end

"""
    distances_between_particles(packing::PolyPacking)

Compute the distance between the *borders* of all pairs of particles in 'packing'.

For a given pair, the distance between the centers of the two particles is computed, by calling
`distances_between_centers`, and then the sum of their radii is subtracted. In this way, the result is
really the distance between the border of each pair of particles.

See also [`distances_between_centers`](@ref)
"""
function distances_between_particles(packing::PolyPacking)
    N = length(packing)
    Rs = get_radii(packing)
    dists = distances_between_centers(packing).data
    for i in 1:N, j in i+1:N
        dists[i, j] -= Rs[i]+Rs[j]
    end
    Symmetric(dists)
end

"""
    check_for_overlaps(packing::PolyPacking{d,T}, tolerance::T)

Apply `check_for_overlaps` to all the particles in 'packing'.
"""
function check_for_overlaps(packing::PolyPacking{d, T}, tolerance::T) where {d, T<:AbstractFloat}
    Rs = get_radii(packing)
    Xs = get_positions(packing)
    check_for_overlaps(Xs, Rs, tolerance)
end
# function check_for_overlaps(packing::PolyPacking{d, T}, tolerance::T) where {d, T<:AbstractFloat}
#     parts = packing.Particles
#     N = length(packing)
#     overlap = false; message = "No overlap is present"; ovlp_indices = (0,0)
#     for i in 1:N, j in i+1:N
#         overlap, gap = check_for_overlaps(parts[i], parts[j], tolerance)
#         if overlap
#             Xi = value.(parts[i].X); Ri = parts[i].R
#             Xj = value.(parts[j].X); Rj = parts[j].R

#             message = "Overlap between particles $i and $j;\nCenters at: $Xi; \t $Xj. \n\t\tradii:  ($R1, $R2); \t Overlap = $gap ."
#             ovlp_indices = (i,j)
#             break
#         end
#     end
#     return overlap, message, ovlp_indices
# end

#################################################
# Other useful functions
################################################
"""
    packing_fraction(d::Int64, R1::Float64, N1::Int64, R2::Float64, N2::Int64, L::Float64=1.0)

Compute the packing fraction of a configuration of d-dimensional hyperspheres, of which 'N1' have radius 'R1', and 'N2' have radius 'R2'. The configuration is assumed to be inside  a box of size 'L'.

See also [`volume`](@ref), [`volume_d_ball`](@ref)
 """
function packing_fraction(d::Int64, R1::T, N1::Int64, R2::T, N2::Int64, L::T=1.0) where T<:AbstractFloat
    return (N1*volume_d_ball(d, R1)+ N2*volume_d_ball(d, R2))/L^d
end
packing_fraction(3, 0.1, 5, 0.3, 8, 1.0)

"""
    packing_fraction(d::Int64, Rs::Vector{Float64}, L::Real=1.0)

Compute the packing fraction of N d-dimensional hyperspheres with radii 'Rs', inside
 a box of size 'L'.

See also [`volume`](@ref), [`volume_d_ball`](@ref)
 """
function packing_fraction(d::Int64, Rs::Vector{Float64}, L::Real=1.0)::Float64
    return sum(volume_d_ball.(d, Rs))/L^d
end
packing_fraction(3, rand(5), 1.0)

"""
    packing_fraction(packing::PolyPacking{d, T})

Compute the packing fraction of a monodisperse packing.

See also [`volume`](@ref), [`volume_d_ball`](@ref)
"""
function packing_fraction(packing::PolyPacking{d, T})::Float64 where {d, T<:AbstractFloat}
    Rs = get_radii(packing)
    L = packing.Particles[1].X[1].L
    packing_fraction(d, Rs, L)
end

max_stars_repo_name: rdhr/CALiPPSO

max_issues_repo_name: rdhr/CALiPPSO

max_forks_repo_name: rdhr/CALiPPSO


text:
"""
Project density on to a staggerred grid using simple interpolation
"""
function get_rhovxI(rhoI::AbstractArray{Float64})
	rhovxI = zero(rhoI);
	get_rhovxI!(rhovxI, rhoI)
	return rhovxI
end
function get_rhovxI!(rhovxI, rhoI::AbstractArray{Float64})
	for ix = 1:size(rhoI, 2)-1
		for iz = 1:size(rhoI,1)
			rhovxI[iz, ix] = 0.5e0 *(rhoI[iz,ix+1] + rhoI[iz,ix])
		end
	end
	# last point, circular boundary conditions
	nx=size(rhoI, 2)
	for iz = 1:size(rhoI,1)
		rhovxI[iz, nx] = 0.5e0 *(rhoI[iz,1] + rhoI[iz,nx])
	end
	return nothing
end # get_rhovxI

"""
Project density on to a staggerred grid using simple interpolation
"""
function get_rhovzI(rhoI::AbstractArray{Float64})
	rhovzI = zero(rhoI);
	get_rhovzI!(rhovzI, rhoI)
	return rhovzI
end
function get_rhovzI!(rhovzI, rhoI::AbstractArray{Float64})
	nz=size(rhoI, 1)
	for ix = 1: size(rhoI, 2)
		for iz = 1: size(rhoI, 1)-1
			rhovzI[iz, ix] =  0.5e0 * (rhoI[iz+1,ix] + rhoI[iz,ix])
		end
		# last point, circular boundary conditions
		rhovzI[nz, ix] =0.5e0 * (rhoI[1,ix] + rhoI[nz,ix])
	end
	return nothing
end # get_rhovzI

function grad_modrr_sprayrr!(grad_modrr_stack,grad_modrhovxI_stack,grad_modrhovzI_stack)
	@simd for i in eachindex(grad_modrr_stack)
		@inbounds grad_modrr_stack[i] = (grad_modrhovxI_stack[i] + grad_modrhovzI_stack[i]) * (0.5)
	end
end

function grad_modrr_sprayrhovxI!(grad_modrr_stack,grad_modrhovxI_stack)
	for ix=1:size(grad_modrr_stack,2)-1
		for iz=1:size(grad_modrr_stack,1)
			@inbounds grad_modrr_stack[iz,ix+1] +=  0.5e0 * grad_modrhovxI_stack[iz,ix]
		end
	end
	nx=size(grad_modrr_stack,2)
	# circular boundary
	for iz=1:size(grad_modrr_stack,1)
		@inbounds grad_modrr_stack[iz,1] +=  0.5e0 * grad_modrhovxI_stack[iz,nx]
	end
end
function grad_modrr_sprayrhovzI!(grad_modrr_stack,grad_modrhovzI_stack)
	nz=size(grad_modrr_stack,1)
	for ix=1:size(grad_modrr_stack,2)
		for iz=1:size(grad_modrr_stack,1)-1
			@inbounds grad_modrr_stack[iz+1,ix] +=  0.5e0 * grad_modrhovzI_stack[iz,ix]
		end
		# circular boundary
		grad_modrr_stack[1,ix] +=  0.5e0 * grad_modrhovzI_stack[nz,ix]
	end
end


text:
using ShipMMG
using Documenter

DocMeta.setdocmeta!(ShipMMG, :DocTestSetup, :(using ShipMMG); recursive = true)

makedocs(;
    modules = [ShipMMG],
    authors = "Taiga MITSUYUKI",
    repo = "https://github.com/ShipMMG/ShipMMG.jl/blob/{commit}{path}#{line}",
    sitename = "ShipMMG.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://ShipMMG.github.io/ShipMMG.jl",
        assets = String[],
    ),
    pages = ["Home" => "index.md"],
)

deploydocs(; repo = "github.com/ShipMMG/ShipMMG.jl")

max_stars_repo_name: mitsuyukiLab/JuliaShipMMG

max_issues_repo_name: mitsuyukiLab/JuliaShipMMG


text: include("penalty.jl")