text:
#Array test

## basics

@test length([1, 2, 3]) == 3
@test nfilled([1, 2, 3]) == 3

a = ones(4)
b = a+a
@test b[1]==2. && b[2]==2. && b[3]==2. && b[4]==2.

@test length((1,)) == 1
@test length((1,2)) == 2

@test isequal(1.+[1,2,3], [2,3,4])
@test isequal([1,2,3].+1, [2,3,4])
@test isequal(1.-[1,2,3], [0,-1,-2])
@test isequal([1,2,3].-1, [0,1,2])

@test isequal(5*[1,2,3], [5,10,15])
@test isequal([1,2,3]*5, [5,10,15])
@test isequal(1./[1,2,5], [1.0,0.5,0.2])
@test isequal([1,2,3]/5, [0.2,0.4,0.6])

a = ones(2,2)
a[1,1] = 1
a[1,2] = 2
a[2,1] = 3
a[2,2] = 4
b = a'
@test a[1,1] == 1. && a[1,2] == 2. && a[2,1] == 3. && a[2,2] == 4.
@test b[1,1] == 1. && b[2,1] == 2. && b[1,2] == 3. && b[2,2] == 4.

a = Array(Float64, 2, 2, 2, 2, 2)
a[1,1,1,1,1] = 10
a[1,2,1,1,2] = 20
a[1,1,2,2,1] = 30

@test a[1,1,1,1,1] == 10
@test a[1,2,1,1,2] == 20
@test a[1,1,2,2,1] == 30

b = reshape(a, (32,))
@test b[1]  == 10
@test b[19] == 20
@test b[13] == 30

b = rand(32)
a = reshape(b, (2, 2, 2, 2, 2))
@test ndims(a) == 5
@test a[2,1,2,2,1] == b[14]
@test a[2,2,2,2,2] == b[end]

sz = (5,8,7)
A = reshape(1:prod(sz),sz...)
@test A[2:6] == [2:6]
tmp = A[1:3,2,2:4]
@test tmp == cat(3,46:48,86:88,126:128)
@test A[:,7:-3:1,5] == [191 176 161; 192 177 162; 193 178 163; 194 179 164; 195 180 165]
tmp = A[:,3:9]
@test tmp == reshape(11:45,5,7)
rng = (2,2:3,2:2:5)
tmp = zeros(Int,map(maximum,rng)...)
tmp[rng...] = A[rng...]
@test  tmp == cat(3,zeros(Int,2,3),[0 0 0; 0 47 52],zeros(Int,2,3),[0 0 0; 0 127 132])

x = rand(2,2)
b = x[1,:]
@test isequal(size(b), (1, 2))
b = x[:,1]
@test isequal(size(b), (2,))

x = rand(5,5)
b = x[2:3,2]
@test b[1] == x[2,2] && b[2] == x[3,2]

B = zeros(4,5)
B[:,3] = 1:4
@test B == [0 0 1 0 0; 0 0 2 0 0; 0 0 3 0 0; 0 0 4 0 0]
B[2,:] = 11:15
@test B == [0 0 1 0 0; 11 12 13 14 15; 0 0 3 0 0; 0 0 4 0 0]
B[[3,1],[2,4]] = [21 22; 23 24]
@test B == [0 23 1 24 0; 11 12 13 14 15; 0 21 3 22 0; 0 0 4 0 0]
B[4,[2,3]] = 7
@test B == [0 23 1 24 0; 11 12 13 14 15; 0 21 3 22 0; 0 7 7 0 0]

@test isequal(reshape(1:27, 3, 3, 3)[1,:], [1  4  7  10  13  16  19  22  25])

a = [3, 5, -7, 6]
b = [4, 6, 2, -7, 1]
ind = findin(a, b)
@test ind == [3,4]

rt = Base.return_types(setindex!, (Array{Int32, 3}, Uint8, Vector{Int}, Float64, Range1{Int}))
@test length(rt) == 1 && rt[1] == Array{Int32, 3}

# sub
A = reshape(1:120, 3, 5, 8)
sA = sub(A, 2, 1:5, :)
@test parent(sA) == A
@test parentindexes(sA) == (2:2, 1:5, 1:8)
@test Base.parentdims(sA) == [1:3]
@test size(sA) == (1, 5, 8)
@test_throws sA[2, 1:8]
@test sA[1, 2, 1:8][:] == [5:15:120]
sA[2:5:end] = -1
@test all(sA[2:5:end] .== -1)
@test all(A[5:15:120] .== -1)
@test strides(sA) == (1,3,15)
@test stride(sA,3) == 15
@test stride(sA,4) == 120
sA = sub(A, 1:3, 1:5, 5)
@test Base.parentdims(sA) == [1:2]
sA[1:3,1:5] = -2
@test all(A[:,:,5] .== -2)
sA[:] = -3
@test all(A[:,:,5] .== -3)
@test strides(sA) == (1,3)
sA = sub(A, 1:3, 3, 2:5)
@test Base.parentdims(sA) == [1:3]
@test size(sA) == (3,1,4)
@test sA == A[1:3,3,2:5]
@test sA[:] == A[1:3,3,2:5][:]
sA = sub(A, 1:2:3, 1:3:5, 1:2:8)
@test Base.parentdims(sA) == [1:3]
@test strides(sA) == (2,9,30)
@test sA[:] == A[1:2:3, 1:3:5, 1:2:8][:]

# sub logical indexing #4763
A = sub([1:10], 5:8)
@test A[A.<7] == [5, 6]
B = reshape(1:16, 4, 4)
sB = sub(B, 2:3, 2:3)
@test sB[sB.>8] == [10, 11]

# slice
A = reshape(1:120, 3, 5, 8)
sA = slice(A, 2, :, 1:8)
@test parent(sA) == A
@test parentindexes(sA) == (2, 1:5, 1:8)
@test Base.parentdims(sA) == [2:3]
@test size(sA) == (5, 8)
@test strides(sA) == (3,15)
@test sA[2, 1:8][:] == [5:15:120]
@test sA[:,1] == [2:3:14]
@test sA[2:5:end] == [5:15:110]
sA[2:5:end] = -1
@test all(sA[2:5:end] .== -1)
@test all(A[5:15:120] .== -1)
sA = slice(A, 1:3, 1:5, 5)
@test Base.parentdims(sA) == [1:2]
@test size(sA) == (3,5)
@test strides(sA) == (1,3)
sA = slice(A, 1:2:3, 3, 1:2:8)
@test Base.parentdims(sA) == [1,3]
@test size(sA) == (2,4)
@test strides(sA) == (2,30)
@test sA[:] == A[sA.indexes...][:]

a = [5:8]
@test parent(a) == a
@test parentindexes(a) == (1:4,)

# issue #4335
@test_throws slice(A, 1:2)
@test_throws slice(A, 1:2, 3:4)
@test_throws slice(A, 1:2, 3:4, 5:6, 7:8)

# Out-of-bounds construction. See #4044
A = rand(7,7)
rng = 1:4
sA = sub(A, 2, rng-1)
@test_throws sA[1,1]
@test sA[1,2] == A[2,1]
sA = sub(A, 2, rng)
B = sub(sA, 1, rng-1)
C = sub(B, 1, rng+1)
@test C == sA
sA = slice(A, 2, rng-1)
@test_throws sA[1]
@test sA[2] == A[2,1]
sA = slice(A, 2, rng)
B = slice(sA, rng-1)
C = sub(B, rng+1)
@test C == sA

# issue #6218 - logical indexing
A = rand(2, 2, 3)
msk = ones(Bool, 2, 2)
msk[2,1] = false
sA = sub(A, :, :, 1)
sA[msk] = 1.0
@test sA[msk] == ones(countnz(msk))

# get
let
    A = reshape(1:24, 3, 8)
    x = get(A, 32, -12)
    @test x == -12
    x = get(A, 14, -12)
    @test x == 14
    x = get(A, (2,4), -12)
    @test x == 11
    x = get(A, (4,4), -12)
    @test x == -12
    X = get(A, -5:5, nan(Float32))
    @test eltype(X) == Float32
    @test isnan(X) == [trues(6),falses(5)]
    @test X[7:11] == [1:5]
    X = get(A, (2:4, 9:-2:-13), 0)
    Xv = zeros(Int, 3, 12)
    Xv[1:2, 2:5] = A[2:3, 7:-2:1]
    @test X == Xv
    X2 = get(A, Vector{Int}[[2:4], [9:-2:-13]], 0)
    @test X == X2
end

## arrays as dequeues
l = {1}
push!(l,2,3,8)
@test l[1]==1 && l[2]==2 && l[3]==3 && l[4]==8
v = pop!(l)
@test v == 8
v = pop!(l)
@test v == 3
@test length(l)==2
unshift!(l,4,7,5)
@test l[1]==4 && l[2]==7 && l[3]==5 && l[4]==1 && l[5]==2
v = shift!(l)
@test v == 4
@test length(l)==4

# concatenation
@test isequal([ones(2,2)  2*ones(2,1)], [1. 1 2; 1 1 2])
@test isequal([ones(2,2), 2*ones(1,2)], [1. 1; 1 1; 2 2])

# typed array literals
X = Float64[1 2 3]
Y = [1. 2. 3.]
@test size(X) == size(Y)
for i = 1:3 @test X[i] === Y[i] end
X = Float64[1;2;3]
Y = [1.,2.,3.]
@test size(X) == size(Y)
for i = 1:3 @test X[i] === Y[i] end
X = Float64[1 2 3; 4 5 6]
Y = [1. 2. 3.; 4. 5. 6.]
@test size(X) == size(Y)
for i = 1:length(X) @test X[i] === Y[i] end

# "end"
X = [ i+2j for i=1:5, j=1:5 ]
@test X[end,end] == 15
@test X[end]     == 15  # linear index
@test X[2,  end] == 12
@test X[end,  2] == 9
@test X[end-1,2] == 8
Y = [2, 1, 4, 3]
@test X[Y[end],1] == 5
@test X[end,Y[end]] == 11

## find, findfirst ##
a = [0,1,2,3,0,1,2,3]
@test find(a) == [2,3,4,6,7,8]
@test find(a.==2) == [3,7]
@test find(isodd,a) == [2,4,6,8]
@test findfirst(a) == 2
@test findfirst(a.==0) == 1
@test findfirst(a.==5) == 0
@test findfirst([1,2,4,1,2,3,4], 3) == 6
@test findfirst(isodd, [2,4,6,3,9,2,0]) == 4
@test findfirst(isodd, [2,4,6,2,0]) == 0

## findn ##

b = findn(ones(2,2,2,2))
@test (length(b[1]) == 16)
@test (length(b[2]) == 16)
@test (length(b[3]) == 16)
@test (length(b[4]) == 16)

#hand made case
a = ([2,1,2],[1,2,2],[2,2,2])
z = zeros(2,2,2)
for i = 1:3
    z[a[1][i],a[2][i],a[3][i]] = 10
end
@test isequal(a,findn(z))

#argmin argmax
@test indmax([10,12,9,11]) == 2
@test indmin([10,12,9,11]) == 3
@test findmin([NaN,3.2,1.8]) == (1.8,3)
@test findmax([NaN,3.2,1.8]) == (3.2,2)
@test findmin([NaN,3.2,1.8,NaN]) == (1.8,3)
@test findmax([NaN,3.2,1.8,NaN]) == (3.2,2)
@test findmin([3.2,1.8,NaN,2.0]) == (1.8,2)
@test findmax([3.2,1.8,NaN,2.0]) == (3.2,1)

## permutedims ##

#keeps the num of dim
p = randperm(5)
q = randperm(5)
a = rand(p...)
b = permutedims(a,q)
@test isequal(size(b), tuple(p[q]...))

#hand made case
y = zeros(1,2,3)
for i = 1:6
    y[i]=i
end

z = zeros(3,1,2)
for i = 1:3
    z[i] = i*2-1
    z[i+3] = i*2
end

#permutes correctly
@test isequal(z,permutedims(y,[3,1,2]))
@test isequal(z,permutedims(y,(3,1,2)))

# of a subarray
a = rand(5,5)
s = sub(a,2:3,2:3)
p = permutedims(s, [2,1])
@test p[1,1]==a[2,2] && p[1,2]==a[3,2]
@test p[2,1]==a[2,3] && p[2,2]==a[3,3]

## ipermutedims ##

tensors = {rand(1,2,3,4),rand(2,2,2,2),rand(5,6,5,6),rand(1,1,1,1)}
for i = tensors
    perm = randperm(4)
    @test isequal(i,ipermutedims(permutedims(i,perm),perm))
    @test isequal(i,permutedims(ipermutedims(i,perm),perm))
end

## unique across dim ##

# All rows and columns unique
A = ones(10, 10)
A[diagind(A)] = shuffle!([1:10])
@test unique(A, 1) == A
@test unique(A, 2) == A

# 10 repeats of each row
B = A[shuffle!(repmat(1:10, 10)), :]
C = unique(B, 1)
@test sortrows(C) == sortrows(A)
@test unique(B, 2) == B
@test unique(B.', 2).' == C

# Along third dimension
D = cat(3, B, B)
@test unique(D, 1) == cat(3, C, C)
@test unique(D, 3) == cat(3, B)

# With hash collisions
immutable HashCollision
    x::Float64
end
Base.hash(::HashCollision) = uint(0)
@test map(x->x.x, unique(map(HashCollision, B), 1)) == C

## reduce ##

z = zeros(2,2,2,2)
for i=1:16
    z[i] = i
end

@test sum(z) == sum(z,(1,2,3,4))[1] == 136

# check variants of summation for type-stability and other issues (#6069)
sum2(itr) = invoke(sum, (Any,), itr)
plus(x,y) = x + y
plus() = 0
sum3(A) = reduce(plus, A)
sum4(itr) = invoke(reduce, (Function, Any), plus, itr)
sum5(A) = reduce(plus, 0, A)
sum6(itr) = invoke(reduce, (Function, Int, Any), plus, 0, itr)
sum7(A) = mapreduce(x->x, plus, A)
sum8(itr) = invoke(mapreduce, (Function, Function, Any), x->x, plus, itr)
sum9(A) = mapreduce(x->x, plus, 0, A)
sum10(itr) = invoke(mapreduce, (Function, Function, Int, Any), x->x,plus,0,itr)
for f in (sum2, sum3, sum4, sum5, sum6, sum7, sum8, sum9, sum10)
    @test sum(z) == f(z)
    @test sum(Int[]) == f(Int[]) == 0
    @test sum(Int[7]) == f(Int[7]) == 7
    if f == sum3 || f == sum4 || f == sum7 || f == sum8
        @test typeof(f(Int8[])) == typeof(f(Int8[1 7]))
    else
        @test typeof(f(Int8[])) == typeof(f(Int8[1])) == typeof(f(Int8[1 7]))
    end
end
@test typeof(sum(Int8[])) == typeof(sum(Int8[1])) == typeof(sum(Int8[1 7]))

prod2(itr) = invoke(prod, (Any,), itr)
@test prod(Int[]) == prod2(Int[]) == 1
@test prod(Int[7]) == prod2(Int[7]) == 7
@test typeof(prod(Int8[])) == typeof(prod(Int8[1])) == typeof(prod(Int8[1 7])) == typeof(prod2(Int8[])) == typeof(prod2(Int8[1])) == typeof(prod2(Int8[1 7]))

v = cell(2,2,1,1)
v[1,1,1,1] = 28.0
v[1,2,1,1] = 36.0
v[2,1,1,1] = 32.0
v[2,2,1,1] = 40.0

@test isequal(v,sum(z,(3,4)))

z = rand(10^6)
let es = sum_kbn(z), es2 = sum_kbn(z[1:10^5])
    @test (es - sum(z)) < es * 1e-13
    cs = cumsum(z)
    @test (es - cs[end]) < es * 1e-13
    @test (es2 - cs[10^5]) < es2 * 1e-13
end
@test_approx_eq sum(sin(z)) sum(sin, z)

@test maximum([4, 3, 5, 2]) == 5
@test minimum([4, 3, 5, 2]) == 2
@test extrema([4, 3, 5, 2]) == (2, 5)

@test isnan(maximum([NaN]))
@test isnan(minimum([NaN]))
@test isequal(extrema([NaN]), (NaN, NaN))

@test maximum([4., 3., NaN, 5., 2.]) == 5.
@test minimum([4., 3., NaN, 5., 2.]) == 2.
@test extrema([4., 3., NaN, 5., 2.]) == (2., 5.)

@test extrema(1:5) == (1,5)

@test any([true false; false false], 2) == [true false]'
@test any([true false; false false], 1) == [true false]

@test all([true true; false true], 2) == [true false]'
@test all([true false; false true], 1) == [false false]

## large matrices transpose ##

for i = 1 : 3
    a = rand(200, 300)

    @test isequal(a', permutedims(a, [2, 1]))
end

## cumsum, cummin, cummax

@test isequal(cummin([1, 2, 5, -1, 3, -2]), [1, 1, 1, -1, -1, -2])
@test isequal(cummax([1, 2, 5, -1, 3, -2]), [1, 2, 5, 5, 5, 5])

@test isequal(cummax([1 0; 0 1], 1), [1 0; 1 1])
@test isequal(cummax([1 0; 0 1], 2), [1 1; 0 1])
@test isequal(cummin([1 0; 0 1], 1), [1 0; 0 0])
@test isequal(cummin([1 0; 0 1], 2), [1 0; 0 0])

@test sum_kbn([1,1e100,1,-1e100]) == 2

begin
    local A, A1, A2, A3, v, v2, cv, cv2, c, R, T
    A = ones(Int,2,3,4)
    A1 = reshape(repmat([1,2],1,12),2,3,4)
    A2 = reshape(repmat([1 2 3],2,4),2,3,4)
    A3 = reshape(repmat([1 2 3 4],6,1),2,3,4)
    @test isequal(cumsum(A),A1)
    @test isequal(cumsum(A,1),A1)
    @test isequal(cumsum(A,2),A2)
    @test isequal(cumsum(A,3),A3)

    R = repeat([1, 2], inner = [1], outer = [1])
    @test R == [1, 2]
    R = repeat([1, 2], inner = [2], outer = [1])
    @test R == [1, 1, 2, 2]
    R = repeat([1, 2], inner = [1], outer = [2])
    @test R == [1, 2, 1, 2]
    R = repeat([1, 2], inner = [2], outer = [2])
    @test R == [1, 1, 2, 2, 1, 1, 2, 2]
    R = repeat([1, 2], inner = [1, 1], outer = [1, 1])
    @test R == [1, 2]''
    R = repeat([1, 2], inner = [2, 1], outer = [1, 1])
    @test R == [1, 1, 2, 2]''
    R = repeat([1, 2], inner = [1, 2], outer = [1, 1])
    @test R == [1 1; 2 2]
    R = repeat([1, 2], inner = [1, 1], outer = [2, 1])
    @test R == [1, 2, 1, 2]''
    R = repeat([1, 2], inner = [1, 1], outer = [1, 2])
    @test R == [1 1; 2 2]

    R = repeat([1 2;
                3 4], inner = [1, 1], outer = [1, 1])
    @test R == [1 2;
                  3 4]
    R = repeat([1 2;
                3 4], inner = [1, 1], outer = [2, 1])
    @test R == [1 2;
                  3 4;
                  1 2;
                  3 4]
    R = repeat([1 2;
                3 4], inner = [1, 1], outer = [1, 2])
    @test R == [1 2 1 2;
                  3 4 3 4]
    R = repeat([1 2;
                3 4], inner = [1, 1], outer = [2, 2])
    @test R == [1 2 1 2;
                  3 4 3 4;
                  1 2 1 2;
                  3 4 3 4]
    R = repeat([1 2;
                3 4], inner = [2, 1], outer = [1, 1])
    @test R == [1 2;
                  1 2;
                  3 4;
                  3 4]
    R = repeat([1 2;
                3 4], inner = [2, 1], outer = [2, 1])
    @test R == [1 2;
                  1 2;
                  3 4;
                  3 4;
                  1 2;
                  1 2;
                  3 4;
                  3 4]
    R = repeat([1 2;
                3 4], inner = [2, 1], outer = [1, 2])
    @test R == [1 2 1 2;
                  1 2 1 2;
                  3 4 3 4;
                  3 4 3 4;]
    R = repeat([1 2;
                3 4], inner = [2, 1], outer = [2, 2])
    @test R == [1 2 1 2;
                  1 2 1 2;
                  3 4 3 4;
                  3 4 3 4;
                  1 2 1 2;
                  1 2 1 2;
                  3 4 3 4;
                  3 4 3 4]
    R = repeat([1 2;
                3 4], inner = [1, 2], outer = [1, 1])
    @test R == [1 1 2 2;
                  3 3 4 4]
    R = repeat([1 2;
                3 4], inner = [1, 2], outer = [2, 1])
    @test R == [1 1 2 2;
                  3 3 4 4;
                  1 1 2 2;
                  3 3 4 4]
    R = repeat([1 2;
                3 4], inner = [1, 2], outer = [1, 2])
    @test R == [1 1 2 2 1 1 2 2;
                  3 3 4 4 3 3 4 4]
    R = repeat([1 2;
                3 4], inner = [1, 2], outer = [2, 2])
    @test R == [1 1 2 2 1 1 2 2;
                  3 3 4 4 3 3 4 4;
                  1 1 2 2 1 1 2 2;
                  3 3 4 4 3 3 4 4]
    R = repeat([1 2;
                3 4], inner = [2, 2], outer = [1, 1])
    @test R == [1 1 2 2;
                  1 1 2 2;
                  3 3 4 4;
                  3 3 4 4]
    R = repeat([1 2;
                3 4], inner = [2, 2], outer = [2, 1])
    @test R == [1 1 2 2;
                  1 1 2 2;
                  3 3 4 4;
                  3 3 4 4;
                  1 1 2 2;
                  1 1 2 2;
                  3 3 4 4;
                  3 3 4 4]
    R = repeat([1 2;
                3 4], inner = [2, 2], outer = [1, 2])
    @test R == [1 1 2 2 1 1 2 2;
                  1 1 2 2 1 1 2 2;
                  3 3 4 4 3 3 4 4;
                  3 3 4 4 3 3 4 4]
    R = repeat([1 2;
                3 4], inner = [2, 2], outer = [2, 2])
    @test R == [1 1 2 2 1 1 2 2;
                  1 1 2 2 1 1 2 2;
                  3 3 4 4 3 3 4 4;
                  3 3 4 4 3 3 4 4;
                  1 1 2 2 1 1 2 2;
                  1 1 2 2 1 1 2 2;
                  3 3 4 4 3 3 4 4;
                  3 3 4 4 3 3 4 4]

    A = reshape([1:8], 2, 2, 2)
    R = repeat(A, inner = [1, 1, 2], outer = [1, 1, 1])
    T = reshape([1:4, 1:4, 5:8, 5:8], 2, 2, 4)
    @test R == T
    A = Array(Int, 2, 2, 2)
    A[:, :, 1] = [1 2;
                  3 4]
    A[:, :, 2] = [5 6;
                  7 8]
    R = repeat(A, inner = [2, 2, 2], outer = [2, 2, 2])
    @test R[1, 1, 1] == 1
    @test R[2, 2, 2] == 1
    @test R[3, 3, 3] == 8
    @test R[4, 4, 4] == 8
    @test R[5, 5, 5] == 1
    @test R[6, 6, 6] == 1
    @test R[7, 7, 7] == 8
    @test R[8, 8, 8] == 8

    A = rand(4,4)
    for s in {A[1:2:4, 1:2:4], sub(A, 1:2:4, 1:2:4)}
        c = cumsum(s, 1)
        @test c[1,1] == A[1,1]
        @test c[2,1] == A[1,1]+A[3,1]
        @test c[1,2] == A[1,3]
        @test c[2,2] == A[1,3]+A[3,3]

        c = cumsum(s, 2)
        @test c[1,1] == A[1,1]
        @test c[2,1] == A[3,1]
        @test c[1,2] == A[1,1]+A[1,3]
        @test c[2,2] == A[3,1]+A[3,3]
    end

    v   = [1,1e100,1,-1e100]*1000
    v2  = [1,-1e100,1,1e100]*1000

    cv  = [1,1e100,1e100,2]*1000
    cv2 = [1,-1e100,-1e100,2]*1000

    @test isequal(cumsum_kbn(v), cv)
    @test isequal(cumsum_kbn(v2), cv2)

    A = [v reverse(v) v2 reverse(v2)]

    c = cumsum_kbn(A, 1)

    @test isequal(c[:,1], cv)
    @test isequal(c[:,3], cv2)
    @test isequal(c[4,:], [2.0 2.0 2.0 2.0]*1000)

    c = cumsum_kbn(A, 2)

    @test isequal(c[1,:], cv2')
    @test isequal(c[3,:], cv')
    @test isequal(c[:,4], [2.0,2.0,2.0,2.0]*1000)

end

@test (1:5)[[true,false,true,false,true]] == [1,3,5]
@test [1:5][[true,false,true,false,true]] == [1,3,5]
@test_throws (1:5)[[true,false,true,false]]
@test_throws (1:5)[[true,false,true,false,true,false]]
@test_throws [1:5][[true,false,true,false]]
@test_throws [1:5][[true,false,true,false,true,false]]
a = [1:5]
a[[true,false,true,false,true]] = 6
@test a == [6,2,6,4,6]
a[[true,false,true,false,true]] = [7,8,9]
@test a == [7,2,8,4,9]
@test_throws (a[[true,false,true,false,true]] = [7,8,9,10])
A = reshape(1:15, 3, 5)
@test A[[true, false, true], [false, false, true, true, false]] == [7 10; 9 12]
@test_throws A[[true, false], [false, false, true, true, false]]
@test_throws A[[true, false, true], [false, true, true, false]]
@test_throws A[[true, false, true, true], [false, false, true, true, false]]
@test_throws A[[true, false, true], [false, false, true, true, false, true]]
A = ones(Int, 3, 5)
@test_throws A[2,[true, false, true, true, false]] = 2:5
A[2,[true, false, true, true, false]] = 2:4
@test A == [1 1 1 1 1; 2 1 3 4 1; 1 1 1 1 1]
@test_throws A[[true,false,true], 5] = [19]
@test_throws A[[true,false,true], 5] = 19:21
A[[true,false,true], 5] = 7
@test A == [1 1 1 1 7; 2 1 3 4 1; 1 1 1 1 7]

B = cat(3, 1, 2, 3)
@test B[:,:,[true, false, true]] == reshape([1,3], 1, 1, 2)  # issue #5454

# issue #2342
@test isequal(cumsum([1 2 3]), [1 2 3])

# set-like operations
@test isequal(union([1,2,3], [4,3,4]), [1,2,3,4])
@test isequal(union(['e','c','a'], ['b','a','d']), ['e','c','a','b','d'])
@test isequal(union([1,2,3], [4,3], [5]), [1,2,3,4,5])
@test isequal(union([1,2,3]), [1,2,3])
@test isequal(union([1,2,3], Int64[]), Int64[1,2,3])
@test isequal(union([1,2,3], Float64[]), Float64[1.0,2,3])
@test isequal(union(Int64[], [1,2,3]), Int64[1,2,3])
@test isequal(union(Int64[]), Int64[])
@test isequal(intersect([1,2,3], [4,3,4]), [3])
@test isequal(intersect(['e','c','a'], ['b','a','d']), ['a'])
@test isequal(intersect([1,2,3], [3,1], [2,1,3]), [1,3])
@test isequal(intersect([1,2,3]), [1,2,3])
@test isequal(intersect([1,2,3], Int64[]), Int64[])
@test isequal(intersect([1,2,3], Float64[]), Float64[])
@test isequal(intersect(Int64[], [1,2,3]), Int64[])
@test isequal(intersect(Int64[]), Int64[])
@test isequal(setdiff([1,2,3,4], [2,5,4]), [1,3])
@test isequal(setdiff([1,2,3,4], [7,8,9]), [1,2,3,4])
@test isequal(setdiff([1,2,3,4], Int64[]), Int64[1,2,3,4])
@test isequal(setdiff([1,2,3,4], [1,2,3,4,5]), Int64[])
@test isequal(symdiff([1,2,3], [4,3,4]), [1,2,4])
@test isequal(symdiff(['e','c','a'], ['b','a','d']), ['e','c','b','d'])
@test isequal(symdiff([1,2,3], [4,3], [5]), [1,2,4,5])
@test isequal(symdiff([1,2,3,4,5], [1,2,3], [3,4]), [3,5])
@test isequal(symdiff([1,2,3]), [1,2,3])
@test isequal(symdiff([1,2,3], Int64[]), Int64[1,2,3])
@test isequal(symdiff([1,2,3], Float64[]), Float64[1.0,2,3])
@test isequal(symdiff(Int64[], [1,2,3]), Int64[1,2,3])
@test isequal(symdiff(Int64[]), Int64[])

# mapslices
begin
    local a,h,i
    a = rand(5,5)
    h = mapslices(v -> hist(v,0:0.1:1)[2], a, 1)
    H = mapslices(v -> hist(v,0:0.1:1)[2], a, 2)
    s = mapslices(sort, a, [1])
    S = mapslices(sort, a, [2])
    for i = 1:5
        @test h[:,i] == hist(a[:,i],0:0.1:1)[2]
        @test vec(H[i,:]) == hist(vec(a[i,:]),0:0.1:1)[2]
        @test s[:,i] == sort(a[:,i])
        @test vec(S[i,:]) == sort(vec(a[i,:]))
    end

    # issue #3613
    b = mapslices(sum, ones(2,3,4), [1,2])
    @test size(b) === (1,1,4)
    @test all(b.==6)

    # issue #5141
    ## Update Removed the version that removes the dimensions when dims==1:ndims(A)
    c1 = mapslices(x-> maximum(-x), a, [])
    @test c1 == -a

    # other types than Number
    @test mapslices(prod,["1" "2"; "3" "4"],1) == ["13" "24"]

    # issue #5177

    c = ones(2,3,4)
    m1 = mapslices(x-> ones(2,3), c, [1,2])
    m2 = mapslices(x-> ones(2,4), c, [1,3])
    m3 = mapslices(x-> ones(3,4), c, [2,3])
    @test size(m1) == size(m2) == size(m3) == size(c)

    n1 = mapslices(x-> ones(6), c, [1,2])
    n2 = mapslices(x-> ones(6), c, [1,3])
    n3 = mapslices(x-> ones(6), c, [2,3])
    n1a = mapslices(x-> ones(1,6), c, [1,2])
    n2a = mapslices(x-> ones(1,6), c, [1,3])
    n3a = mapslices(x-> ones(1,6), c, [2,3])
    @test size(n1a) == (1,6,4) && size(n2a) == (1,3,6)  && size(n3a) == (2,1,6)
    @test size(n1) == (6,1,4) && size(n2) == (6,3,1)  && size(n3) == (2,6,1)
end

# single multidimensional index
let
    a = rand(6,6)
    I = [1 4 5; 4 2 6; 5 6 3]
    a2 = a[I]
    @test size(a2) == size(I)
    for i = 1:length(a2)
        @test a2[i] == a[I[i]]
    end
    a = [1,3,5]
    b = [1 3]
    a[b] = 8
    @test a == [8,3,8]
end

# assigning an array into itself
a = [1,3,5]
b = [3,1,2]
a[b] = a
@test a == [3,5,1]

# lexicographic comparison
@test lexcmp([1.0], [1]) == 0
@test lexcmp([1], [1.0]) == 0
@test lexcmp([1, 1], [1, 1]) == 0
@test lexcmp([1, 1], [2, 1]) == -1
@test lexcmp([2, 1], [1, 1]) == 1
@test lexcmp([1, 1], [1, 2]) == -1
@test lexcmp([1, 2], [1, 1]) == 1
@test lexcmp([1], [1, 1]) == -1
@test lexcmp([1, 1], [1]) == 1

# sort on arrays
begin
    local a = rand(3,3)

    asr = sortrows(a)
    @test lexless(asr[1,:],asr[2,:])
    @test lexless(asr[2,:],asr[3,:])

    asc = sortcols(a)
    @test lexless(asc[:,1],asc[:,2])
    @test lexless(asc[:,2],asc[:,3])

    asr = sortrows(a, rev=true)
    @test lexless(asr[2,:],asr[1,:])
    @test lexless(asr[3,:],asr[2,:])

    asc = sortcols(a, rev=true)
    @test lexless(asc[:,2],asc[:,1])
    @test lexless(asc[:,3],asc[:,2])

    as = sort(a, 1)
    @test issorted(as[:,1])
    @test issorted(as[:,2])
    @test issorted(as[:,3])

    as = sort(a, 2)
    @test issorted(as[1,:])
    @test issorted(as[2,:])
    @test issorted(as[3,:])
end

# fill
@test fill!(Array(Float64,1),-0.0)[1] === -0.0
A = ones(3,3)
S = sub(A, 2, 1:3)
fill!(S, 2)
S = sub(A, 1:2, 3)
fill!(S, 3)
@test A == [1 1 3; 2 2 3; 1 1 1]
rt = Base.return_types(fill!, (Array{Int32, 3}, Uint8))
@test length(rt) == 1 && rt[1] == Array{Int32, 3}

# splice!
for idx in {1, 2, 5, 9, 10, 1:0, 2:1, 1:1, 2:2, 1:2, 2:4, 9:8, 10:9, 9:9, 10:10,
            8:9, 9:10, 6:9, 7:10}
    for repl in {[], [11], [11,22], [11,22,33,44,55]}
        a = [1:10]; acopy = copy(a)
        @test splice!(a, idx, repl) == acopy[idx]
        @test a == [acopy[1:(first(idx)-1)], repl, acopy[(last(idx)+1):end]]
    end
end

# deleteat!
for idx in {1, 2, 5, 9, 10, 1:0, 2:1, 1:1, 2:2, 1:2, 2:4, 9:8, 10:9, 9:9, 10:10,
            8:9, 9:10, 6:9, 7:10}
    a = [1:10]; acopy = copy(a)
    @test deleteat!(a, idx) == [acopy[1:(first(idx)-1)], acopy[(last(idx)+1):end]]
end
a = [1:10]
@test deleteat!(a, [1,3,5,7:10...]) == [2,4,6]

# comprehensions
X = [ i+2j for i=1:5, j=1:5 ]
@test X[2,3] == 8
@test X[4,5] == 14
@test isequal(ones(2,3) * ones(2,3)', [3. 3.; 3. 3.])
@test isequal([ [1,2] for i=1:2, : ], [1 2; 1 2])
# where element type is a Union. try to confuse type inference.
foo32_64(x) = (x<2) ? int32(x) : int64(x)
boo32_64() = [ foo32_64(i) for i=1:2 ]
let a36 = boo32_64()
    @test a36[1]==1 && a36[2]==2
end
@test isequal([1,2,3], [a for (a,b) in enumerate(2:4)])
@test isequal([2,3,4], [b for (a,b) in enumerate(2:4)])

@test_throws (10.^[-1])[1] == 0.1
@test (10.^[-1.])[1] == 0.1

# reverse
@test reverse([2,3,1]) == [1,3,2]
@test reverse([1:10],1,4) == [4,3,2,1,5,6,7,8,9,10]
@test reverse([1:10],3,6) == [1,2,6,5,4,3,7,8,9,10]
@test reverse([1:10],6,10) == [1,2,3,4,5,10,9,8,7,6]
@test reverse(1:10,1,4) == [4,3,2,1,5,6,7,8,9,10]
@test reverse(1:10,3,6) == [1,2,6,5,4,3,7,8,9,10]
@test reverse(1:10,6,10) == [1,2,3,4,5,10,9,8,7,6]
@test reverse!([1:10],1,4) == [4,3,2,1,5,6,7,8,9,10]
@test reverse!([1:10],3,6) == [1,2,6,5,4,3,7,8,9,10]
@test reverse!([1:10],6,10) == [1,2,3,4,5,10,9,8,7,6]

# flipdim
@test isequal(flipdim([2,3,1], 1), [1,3,2])
@test isequal(flipdim([2,3,1], 2), [2,3,1])
@test isequal(flipdim([2 3 1], 1), [2 3 1])
@test isequal(flipdim([2 3 1], 2), [1 3 2])
@test_throws flipdim([2,3,1], -1)
@test isequal(flipdim(1:10, 1), 10:-1:1)
@test isequal(flipdim(1:10, 2), 1:10)
@test_throws flipdim(1:10, -1)
@test isequal(flipdim(Array(Int,0,0),1), Array(Int,0,0))  # issue #5872

# issue 4228
A = [[i i; i i] for i=1:2]
@test cumsum(A) == Any[[1 1; 1 1], [3 3; 3 3]]
@test cumprod(A) == Any[[1 1; 1 1], [4 4; 4 4]]

# PR #4627
A = [1,2]
@test append!(A, A) == [1,2,1,2]
@test prepend!(A, A) == [1,2,1,2,1,2,1,2]

# cases where shared arrays can/can't be grown
A = [1 3;2 4]
B = reshape(A, 4)
@test push!(B,5) == [1,2,3,4,5]
@test pop!(B) == 5
C = reshape(B, 1, 4)
@test_throws push!(C, 5)

A = [NaN]; B = [NaN]
@test !(A==A)
@test isequal(A,A)
@test A===A
@test !(A==B)
@test isequal(A,B)
@test A!==B

# complete testsuite for reducedim

include("reducedim.jl")
# Inferred types
Nmax = 3 # TODO: go up to CARTESIAN_DIMS+2 (currently this exposes problems)
for N = 1:Nmax
    #indexing with (Range1, Range1, Range1)
    args = ntuple(N, d->Range1{Int})
    @test Base.return_types(getindex, tuple(Array{Float32, N}, args...)) == {Array{Float32, N}}
    @test Base.return_types(getindex, tuple(BitArray{N}, args...)) == {BitArray{N}}
    @test Base.return_types(setindex!, tuple(Array{Float32, N}, Array{Int, 1}, args...)) == {Array{Float32, N}}
    # Indexing with (Range1, Range1, Float64)
    args = ntuple(N, d->d<N ? Range1{Int} : Float64)
    N > 1 && @test Base.return_types(getindex, tuple(Array{Float32, N}, args...)) == {Array{Float32, N-1}}
    N > 1 && @test Base.return_types(getindex, tuple(BitArray{N}, args...)) == {BitArray{N-1}}
    N > 1 && @test Base.return_types(setindex!, tuple(Array{Float32, N}, Array{Int, 1}, args...)) == {Array{Float32, N}}
end

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
## This file autogenerated by BinaryProvider.write_deps_file().
## Do not edit.
##
## Include this file within your main top-level source, and call
## `check_deps()` from within your module's `__init__()` method

if isdefined((@static VERSION < v"0.7.0-DEV.484" ? current_module() : @__MODULE__), :Compat)
    import Compat.Libdl
elseif VERSION >= v"0.7.0-DEV.3382"
    import Libdl
end
const libextraf = joinpath(dirname(@__FILE__), "usr/lib/libextraf.so")
function check_deps()
    global libextraf
    if !isfile(libextraf)
        error("$(libextraf) does not exist, Please re-run Pkg.build(\"MakeFshared\"), and restart Julia.")
    end

    if Libdl.dlopen_e(libextraf) in (C_NULL, nothing)
        error("$(libextraf) cannot be opened, Please re-run Pkg.build(\"MakeFshared\"), and restart Julia.")
    end

end

ext: jl


text:
using DynamicalSystems, MakieLayout
export interactive_poincaresos

"""
    interactive_poincaresos(cds, plane, idxs, complete; kwargs...)
Launch an interactive application for exploring a Poincaré surface of section (PSOS)
of the continuous dynamical system `cds` using `poincaresos` from DynamicalSystems.jl.

The `plane` can only be the `Tuple` type accepted by `DynamicalSystems.poincaresos`,
i.e. `(i, r)` for the `i`th variable crossing the value `r`. `idxs` gives the two
indices of the variables to be displayed, since the PSOS plot is always a 2D scatterplot.
I.e. `idxs = (1, 2)` will plot the 1st versus 2nd variable of the PSOS. It follows
that `plane[1] ∉ idxs` must be true.

`complete` is a three-argument **function** that completes the new initial state
during interactive use, see below.

The function returns: an observable containing the latest initial `state`
and the scene.

## Keyword Arguments
* `direction, rootkw` : Same use as in `DynamicalSystems.poincaresos`.
* `tfinal = (1000.0, 10.0^4)` : A 2-element tuple for the range of values
  for the total integration time (chosen interactively).
* `color` : A **function** of the system's initial condition, that returns a color to
  plot the new points with. The color must be `RGBf0/RGBAf0`.
   A random color is chosen by default.
* `labels = ("u₁" , "u₂")` : Axis labels (you can change them youself `ax.xlabel = ...`)
* `scatterkwargs = ()`: Named tuple of keywords passed to `scatter`.
* `diffeq...` : Any extra keyword arguments are passed into `init` of DiffEq.

## Interaction
The application is a standard scatterplot, which shows the PSOS of the system,
initially using the system's `u0`. Two sliders control the total evolution time
and the size of the marker points (which is always in pixels).

Upon clicking within the bounds of the scatter plot your click is transformed into
a new initial condition, which is further evolved and its PSOS is computed and then
plotted into the scatter plot.

Your click is transformed into a full `D`-dimensional initial condition through
the function `complete`. The first two arguments of the function are the positions
of the click on the PSOS. The third argument is the value of the variable the PSOS
is defined on. To be more exact, this is how the function is called:
```julia
x, y = mouseclick; z = plane[2]
newstate = complete(x, y, z)
```
The `complete` function can throw an error for ill-conditioned `x, y, z`.
This will be properly handled instead of breaking the application.
This `newstate` is also given to the function `color` that
gets a new color for the new points.
"""
function interactive_poincaresos(ds::ContinuousDynamicalSystem{IIP, S, D}, plane, idxs, complete;
         # PSOS kwargs:
         direction = -1,
         tfinal = (1000.0, 10.0^4),
         rootkw = (xrtol = 1e-6, atol = 1e-6),
         # AbstractPlotting kwargs:
         color = randomcolor,
         scatterkwargs = (),
         labels = ("u₁" , "u₂"),
         # DiffEq kwargs:
         diffeq...) where {IIP, S, D}

    @assert typeof(plane) <: Tuple
    @assert length(idxs) == 2
    @assert eltype(idxs) == Int
    @assert plane[1] ∉ idxs
    u0 = get_state(ds)

    # This is the low-level call of poincaresos:
    ChaosTools._check_plane(plane, D)
    integ = integrator(ds, u0; diffeq...)
    planecrossing = PlaneCrossing(plane, direction > 0)
    i = SVector{2, Int}(idxs)

    scene, layout = layoutscene(resolution = (1000, 800))

    T_slider, m_slider = _add_psos_controls!(scene, layout, tfinal)
    ax = layout[0, :] = LAxis(scene)

    # Initial Section
    data = poincaresos(integ, planecrossing, T_slider[], 0.0, i, rootkw)
    length(data) == 0 && error(ChaosTools.PSOS_ERROR)

    positions_node = Observable(data)
    colors = (c = color(u0); [c for i in 1:length(data)])
    colors_node = Observable(colors)
    scplot = scatter!(
        ax, positions_node, color = colors_node,
        markersize = lift(o -> o*px, m_slider), marker = MARKER, scatterkwargs...
    )

    ax.xlabel, ax.ylabel = labels
    laststate = Observable(u0)

    # Interactive clicking on the psos:
    spoint = select_point(ax.scene)
    on(spoint) do pos
        x, y = pos; z = plane[2] # third variable comes from plane
        newstate = try
           complete(x, y, z)
        catch err
           @error "Could not get state, got error: " exception=err
           return
        end

        reinit!(integ, newstate)
        data = poincaresos(integ, planecrossing, T_slider[], 0.0, i, rootkw)
        positions = positions_node[]; colors = colors_node[]
        append!(positions, data)
        c = color(newstate)
        append!(colors, fill(c, length(data)))
        # update all the observables with Array as value:
        positions_node[], colors_node[], laststate[] = positions, colors, newstate
    end
    display(scene)
    return laststate, scene
end

function _add_psos_controls!(scene, layout, tfinal)
    T_slider = LSlider(scene, range = range(tfinal[1], tfinal[2], length = 1000))
    T_text_prev = LText(scene, "T =", halign = :right)
    T_text_after = LText(scene, lift(a -> "$(round(a))", T_slider.value), halign = :left)
    m_slider = LSlider(scene, range = 10.0 .^ range(-1, 2, length = 100), startvalue = 10)
    m_text_prev = LText(scene, "ms =", halign=:right)
    sublayout = GridLayout()
    layout[1, :] = sublayout
    sublayout[:h] = [T_text_prev, T_slider, T_text_after, LText(scene, "   "),
                    m_text_prev, m_slider]

    colsize!(sublayout, 6, Auto(0.5))
    return T_slider.value, m_slider.value
end


text:
export zoo_vtree, zoo_vtree_file,
    VtreeFormat, DotFormat

struct VtreeFormat <: FileFormat end
struct DotFormat <: FileFormat end

##############################################
# Read Vtrees
##############################################

zoo_vtree_file(name) =
    artifact"circuit_model_zoo" * zoo_version * "/vtrees/$name"

"""
    zoo_vtree(name)

Loads Vtree file with given name from model zoo. See https://github.com/UCLA-StarAI/Circuit-Model-Zoo.
"""
function zoo_vtree(name, ::Type{V}=PlainVtree)::V where V<:Vtree
    read(zoo_vtree_file(name), V)
end

const vtree_grammar = raw"""
    start: header (_NL node)+ _NL?

    header : "vtree" _WS INT

    node : "L" _WS INT _WS INT -> leaf
         | "I" _WS INT _WS INT _WS INT -> inode

    %import common.INT
    %import common.WS_INLINE -> _WS
    %import common.NEWLINE -> _NL
    """ * dimacs_comments

const vtree_parser = Lark(vtree_grammar)

mutable struct VtreeParse{V <: Vtree} <: JuiceTransformer
    nodes::Dict{String,V}
    VtreeParse{V}() where V = new{V}(Dict{String,V}())
end

@inline_rule header(t::VtreeParse,x) = Base.parse(Int,x)

@rule leaf(t::VtreeParse{V},x) where V = begin
    t.nodes[x[1]] = V(Base.parse(Var,x[2]))
end

@rule inode(t::VtreeParse{V},x) where V =
    t.nodes[x[1]] = V(t.nodes[x[2]], t.nodes[x[3]])

@rule start(t::VtreeParse, x) = begin
    @assert num_nodes(x[end]) == x[1]
    x[end]
end

function Base.parse(::Type{V}, str,
                    ::VtreeFormat = VtreeFormat()) where V <: Vtree
    ast = Lerche.parse(vtree_parser, str)
    Lerche.transform(VtreeParse{V}(), ast)
end

function Base.parse(::Type{Dict{String,V}}, str,
                    ::VtreeFormat = VtreeFormat()) where V <: Vtree
    ast = Lerche.parse(vtree_parser, str)
    transformer = VtreeParse{V}()
    Lerche.transform(transformer, ast)
    transformer.nodes
end

Base.read(io::IO, ::Type{V},
          ::VtreeFormat = VtreeFormat()) where V <: Vtree =
    parse(V, read(io, String), VtreeFormat())

Base.read(io::IO, ::Type{Dict{String,V}},
          ::VtreeFormat = VtreeFormat()) where V <: Vtree =
    parse(Dict{String,V}, read(io, String), VtreeFormat())

##############################################
# Write Vtrees
##############################################

const VTREE_FORMAT = """c this file was saved by LogicCircuits.jl
c ids of vtree nodes start at 0
c ids of variables start at 1
c vtree nodes appear bottom-up, children before parents
c
c file syntax:
c vtree number-of-nodes-in-vtree
c L id-of-leaf-vtree-node id-of-variable
c I id-of-internal-vtree-node id-of-left-child id-of-right-child
c"""
#end vtree format header

"""
    write(file::AbstractString, vtree::PlainVtree)

Saves a vtree in the given file path based on file format.
Supported formats:
* ".vtree" for Vtree files
* ".dot" for dot files
"""
function Base.write(file::AbstractString, vtree::Vtree)
    if endswith(file,".vtree")
        write(file, vtree, VtreeFormat())
    elseif endswith(file, ".dot")
        write(file, vtree, DotFormat())
    else
        throw("Unsupported file extension in $file: choose either *.vtree or *.dot")
    end
end

function Base.write(io::IO, vtree::Vtree, format = VtreeFormat())

    labeling = label_nodes(vtree)
    map!(x -> x-1, values(labeling)) # vtree nodes are 0-based indexed

    if format == VtreeFormat()
        println(io, VTREE_FORMAT)
        println(io, "vtree $(num_nodes(vtree))")
        foreach(vtree) do n
            if isleaf(n)
                println(io, "L $(labeling[n]) $(n.var)")
            else
                @assert isinner(n)
                println(io, "I $(labeling[n]) $(labeling[n.left]) $(labeling[n.right])")
            end
        end

    elseif format == DotFormat()
        println(io,"strict graph vtree { node [shape=point]; splines=false;")
        # reverse order is better for dot
        foreach_down(vtree) do n
            if isleaf(n)
                println(io, "$(labeling[n]) [label=$(n.var), shape=\"plaintext\"]")
            else
                @assert isinner(n)
                println(io, "$(labeling[n]) -- $(labeling[n.right])")
                println(io, "$(labeling[n]) -- $(labeling[n.left])")
            end
        end
        println(io, "}")
    else
        error("Unknown vtree format: $format")
    end

    # return the labeling to be reused by consumers of this vtree file
    labeling
end

max_issues_repo_name: Juice-jl/LogicCircuits

max_forks_repo_name: Juice-jl/LogicCircuits


text:
using StructuresKit

#A Ix Iy J Cw xo yo
SectionProperties = [(272.0,363370.0,64100.0, 188.0, 122720891.0, -32.59, 0.0)]

#E  ν
MaterialProperties = [(200000,0.30)]

#kx ky kϕ hx hy
Springs = [(0.0, 0.0, 0.0, 0.0, 0.0)]

#member information
#L dL SectionProperties MaterialProperties Springs
MemberDefinitions = [(2438.0,2438.0/12,1,1,1)]

#P qx qy ax ay
Loads = [(10000*ones(13)),(2*ones(13)), (0.0*ones(13)),(0.0*ones(13)),(-1.0*ones(13))]

#end boundary conditions
#type=1 u''=v''=ϕ''=0 (simply supported), type=2 u'=v'=ϕ'=0  (fixed), type=3 u''=v''=ϕ''=u'''=v'''=ϕ'''=0 (free end, e.g., a cantilever)
EndBoundaryConditions = [1 1]

#supports
#location where u=v=ϕ=0
Supports = [0.0 2438.0]

u, v, ϕ, properties = BeamColumn.solve(MemberDefinitions, SectionProperties, MaterialProperties, Loads, Springs, EndBoundaryConditions, Supports)

using Plots

plot(properties.z, u)


text:
module GDBM

export DBM

include(joinpath("..", "deps", "deps.jl"))

__init__() = check_deps()

include("types.jl")
include("low_level.jl")
include("dict.jl")

end


text:
# Submit coverage to Coveralls.io
export Coveralls
"""
Coverage.Coveralls Module

This module provides functionality to push coverage information to the Coveralls
web service. It exports the `submit` and `submit_token` methods.
"""
module Coveralls
    using Coverage
    using HTTP
    using JSON
    using Compat
    using Compat.LibGit2
    using MbedTLS

    export submit, submit_token

    #=
    JSON structure for Coveralls
    Accessed 2015/07/24:
    https://coveralls.zendesk.com/hc/en-us/articles/201774865-API-Introduction
    {
      "service_job_id": "1234567890",
      "service_name": "travis-ci",
      "source_files": [
        {
          "name": "example.rb",
          "source": "def four\n  4\nend",
          "coverage": [null, 1, null]
        },
        {
          "name": "lib/two.rb",
          "source": "def seven\n  eight\n  nine\nend",
          "coverage": [null, 1, 0, null]
        }
      ]
    }
    =#

    # to_json
    # Convert a FileCoverage instance to its Coveralls JSON representation
    to_json(fc::FileCoverage) = Dict("name"          => fc.filename,
                                     "source_digest" => digest(MD_MD5, fc.source, "secret"),
                                     "coverage"      => fc.coverage)

    # Format the body argument to HTTP.post
    makebody(data::Dict) =
        Dict("json_file" => HTTP.Multipart("json_file", IOBuffer(JSON.json(data)),
                                           "application/json"))

    """
        submit(fcs::Vector{FileCoverage}; kwargs...)

    Take a vector of file coverage results (produced by `process_folder`),
    and submits them to Coveralls. Assumes that this code is being run
    on TravisCI or AppVeyor. If running locally, use `submit_token`.
    """
    function submit(fcs::Vector{FileCoverage}; kwargs...)
        verbose = true
        for (k,v) in kwargs
            if k == :verbose
                verbose = v
            end
        end

        if lowercase(get(ENV, "APPVEYOR", "false")) == "true"
            # Submission from AppVeyor requires a REPO_TOKEN environment variable
            data = Dict("service_job_id"    => ENV["APPVEYOR_JOB_ID"],
                        "service_name"      => "appveyor",
                        "source_files"      => map(to_json, fcs),
                        "repo_token"        => ENV["REPO_TOKEN"])

        elseif lowercase(get(ENV, "TRAVIS", "false")) == "true"
            data = Dict("service_job_id"    => ENV["TRAVIS_JOB_ID"],
                        "service_name"      => "travis-ci",
                        "source_files"      => map(to_json, fcs))
        elseif lowercase(get(ENV, "JENKINS", "false")) == "true"
            data = Dict("service_job_id"    => ENV["BUILD_ID"],
                        "service_name"      => "jenkins-ci",
                        "source_files"      => map(to_json, fcs),
                        "repo_token"        => ENV["REPO_TOKEN"])

            data["git"] = query_git_info()

            # get the name of the branch if not a pull request
            if get(ENV, "CI_PULL_REQUEST", "false") == "false"
                data["git"]["branch"] = split(ENV["GIT_BRANCH"], "/")[2]
            end
        else
            error("No compatible CI platform detected")
        end

        if verbose
            println("Submitting data to Coveralls...")
        end

        url = "https://coveralls.io/api/v1/jobs"
        body = HTTP.Form(makebody(data))
        headers = ["Content-Type" => "multipart/form-data; boundary=$(body.boundary)"]
        req = HTTP.post(url, headers, body)

        if verbose
            println("Result of submission:")
            println(String(req.body))
        end
    end

    # query_git_info
    # Pulls information about the repository that isn't available if we
    # are running somewhere other than TravisCI
    function query_git_info(dir=pwd())
        repo            = LibGit2.GitRepo(dir)
        head            = LibGit2.head(repo)
        head_cmt        = LibGit2.peel(head)
        head_oid        = LibGit2.GitHash(head_cmt)
        commit_sha      = string(head_oid)
        author_name     = string(LibGit2.author(head_cmt).name)
        author_email    = string(LibGit2.author(head_cmt).email)
        committer_name  = string(LibGit2.committer(head_cmt).name)
        committer_email = string(LibGit2.committer(head_cmt).email)
        message         = LibGit2.message(head_cmt)
        remote          = ""
        branch          = LibGit2.shortname(head)

        if branch != "HEAD" # if repo is not in detached state
            LibGit2.with(LibGit2.get(LibGit2.GitRemote, repo, branch)) do remote
                remote = LibGit2.url(remote)
            end
        end
        LibGit2.close(repo)

        return Dict(
            "branch"    => branch,
            "remotes"   => [
                Dict(
                    "name"  => "origin",
                    "url"   => remote
                )
            ],
            "head" => Dict(
                "id" => commit_sha,
                "author_name"       => author_name,
                "author_email"      => author_email,
                "committer_name"    => committer_name,
                "committer_email"   => committer_email,
                "message"           => message
            )
        )
    end

    """
        submit_token(fcs::Vector{FileCoverage}, git_info=query_git_info; kwargs...)

    Take a `Vector` of file coverage results (produced by `process_folder`),
    and submits them to Coveralls. For submissions not from TravisCI.

    git_info can be either a `Dict` or a function that returns a `Dict`.
    """
    function submit_token(fcs::Vector{FileCoverage}, git_info=query_git_info; kwargs...)
        data = Dict("repo_token" => ENV["REPO_TOKEN"],
                    "source_files" => map(to_json, fcs))

        verbose = true
        for (k,v) in kwargs
            if k == :verbose
                verbose = v
            end
        end

        # Attempt to parse git info via git_info, unless the user explicitly disables it by setting git_info to nothing
        try
            if isa(git_info, Function)
                data["git"] = git_info()
            elseif isa(git_info, Dict)
                data["git"] = git_info
            end
        catch
        end

        r = HTTP.post("https://coveralls.io/api/v1/jobs", body=makebody(data))

        if verbose
            println("Result of submission:")
            println(String(r.body))
        end
    end
end  # module Coveralls


text:
using Random
using DABenchmarks
using LaTeXStrings

Random.seed!(0)

using Gurobi
using LShapedSolvers

function prepare_dbenchmark(nscenarios::Integer, nsamples::Integer; timeout::Integer = 1000)
    dls = LShapedSolver(GurobiSolver(OutputFlag=0), log=false)
    dtr = LShapedSolver(GurobiSolver(OutputFlag=0), regularization = :tr, log=false)
    solvers = [dls,dtr]
    solvernames = ["Distributed L-shaped","Distributed L-shaped with trust-region"]
    # Create Day-ahead benchmark
    return Scaling(solvernames, solvers, nscenarios, nsamples; timeout = timeout)
end
@info "Running L-shaped benchmarks on 1 worker..."
dab = prepare_dbenchmark(1000, 10); warmup_benchmarks!(dab); run_benchmarks!(dab); save_results!(dab, "ls_scaling_1.json")
exit()

max_stars_repo_name: martinbiel/mpc-reproducibility

max_issues_repo_name: martinbiel/mpc-reproducibility

max_forks_repo_name: martinbiel/mpc-reproducibility


text:
export  w_legendre,
        build_w_jacobi,
        build_w_jacobi01,
        w_jacobi,
        w_jacobi01,
        w_laguerre,
        w_hermite,
        build_w_genhermite,
        build_w_genlaguerre,
        w_meixner_pollaczek,
        build_w_meixner_pollaczek,

        w_gaussian,
        w_uniform01,
        w_logistic,
        w_genhermite,
        build_w_beta,
        build_w_gamma

function w_legendre(t)
    -1. <= t <= 1. ? 1. : throw(error("$t not in support"))
end

function build_w_jacobi(a,b)
    @assert a > -1. && b > -1. "Invalid shape parameters"
    # return w(t) = (1-Float64(t))^a*(1+Float64(t))
    return t->w_jacobi(t,a,b)
end
function w_jacobi(t,a,b)
    -1. <=t<= 1. ? ((1-t)^a*(1+t)^b) : throw(error("$t not in support"))
end

function build_w_jacobi01(a,b)
    @assert a > -1. && b > -1. "Invalid shape parameters"
    # return w(t) = (1-Float64(t))^a*(1+Float64(t))
    return t->w_jacobi01(t,a,b)
end
function w_jacobi01(t,a,b)
    0. <=t<= 1. ? ((1-t)^a*t^b) : throw(error("$t not in support"))
end

function w_hermite(t)
    exp(-t^2)
end

function build_w_genhermite(mu)
    @assert mu>-0.5 "invalid parameter ($mu !> -0.5)"
    return t->w_genhermite(t,mu)
end
function w_genhermite(t,μ::Float64)
    abs(t)^(2*μ)*exp(-t^2)
end
function build_w_genlaguerre(a)
    @assert a>-1 "Invalid shape parameter"
    # return w(t) = Float64(t)^a*exp(-Float64(t))
    return t->w_genlaguerre(t,a)
end
function w_genlaguerre(t,a)
    t>=0. ? Float64(t)^a*exp(-Float64(t)) : throw(error("$t not in support"))
end

function w_laguerre(t)
    t>=0. ? exp(-t) : throw(error("$t not in support"))
end

function w_meixner_pollaczek(t,lambda,phi)
    1/(2*pi)*exp((2*phi-pi)*t)*abs(gamma(lambda+im*t))^2
end

function build_w_meixner_pollaczek(lambda,phi)
    @assert lambda>0 "lambda has to be positive"
    @assert 0<phi<pi "phi has to be between 0 and pi"
    t->w_meixner_pollaczek(t,lambda,phi)
end

##################################################
# probability density functions
function w_gaussian(t)
    1/(sqrt(2*pi))*exp(-0.5*t^2)
end

function build_w_beta(α,β)
    @assert α>-0. && β>-0. "Invalid shape parameters"
    return t->w_beta(t,α,β)
end
function w_beta(t,α,β)
    # 0. <=t<= 1. ? (t^(α-1)*(1-t)^(β-1)/beta(α,β)) : (error("$t not in support"))
    t^(α-1)*(1-t)^(β-1)/beta(α,β)
end

function build_w_gamma(α)
    @assert α>0 "Invalid shape parameter"
    return t->w_gamma(t,α)
end
function w_gamma(t,α)
    t>=0. ? (1/gamma(α)*Float64(t)^(α-1)*exp(-t)) : (error("$t not in support"))
end

function w_uniform01(t)
    0. <=t<= 1. ? 1. : (error("$t not in support"))
end

function w_logistic(t::Float64)
    0.25*sech(0.5t)^2
end


text:
export SchobelZhuParams, evaluateCharFunc, evaluateLogCharFunc, computeCumulants

struct SchobelZhuParams{T}
    v0::T
    κ::T
    θ::T
    ρ::T
    σ::T
end
DefaultCharFunc(params::SchobelZhuParams{Float64}) = DefaultCharFunc{SchobelZhuParams{Float64},Complex}(params)

cinf(params::SchobelZhuParams{T}, τ::T) where {T} = (params.v0^2/params.σ+τ*params.σ)/2 *sqrt(1-params.ρ^2)

#CC the Nemo arbField or Complex type. z: Nemo complex number or Complex number.
@inline function evaluateCharFunc(p::CharFunc{SchobelZhuParams{T},CR}, z::CT, τ::T)::CR where {T,CR,CT}
    return exp(evaluateLogCharFunc(p, z, τ))
end

evaluateLogCharFunc(p::CharFunc{SchobelZhuParams{T},CR}, z::CT, τ::T) where {T,CR,CT} = evaluateLogCharFuncZhu(p, z, τ)

function evaluateLogCharFuncLK(cf::CharFunc{SchobelZhuParams{T},CR}, z::CT, τ::T)::CR where {T,CR,CT}
    #Lord and Kahl formulation
    p = model(cf) #model(cf), abstract cf with concrete types: LK, Zhu.
    v0 = p.v0
    κ = p.κ
    θ = p.θ
    ρ = p.ρ
    σ = p.σ
    cc1 = oneim(cf)
    α = -(z * (z + cc1)) / 2
    β = 2 * κ - 2 * σ * ρ * cc1 * z
    γ = 2 * σ^2
    D = sqrt(β^2 - 4 * α * γ)
    G = (β - D) / (β + D)
    Dt = D * τ
    eDT = exp(-Dt)
    oneGeDT = (1 - G * eDT)
    l = log(oneGeDT / (1 - G))
    sqrteDT = exp(-D * τ / 2) #do not use sqrt(eDT) as it may be the other root
    As =
        (β - D) * (κ * θ)^2 / (2 * σ^2 * D^3) * (
            β * (Dt - 4) +
            D * (Dt - 2) +
            4 * sqrteDT * ((D^2 - 2 * β^2) / (β + D) * sqrteDT + 2 * β) / oneGeDT
        )
    A = As + (β - D) * τ / 4 - l / 2
    Bs = 2 * κ * θ * (β - D) * (1 - sqrteDT)^2 / (D * γ * oneGeDT)
    Bv = (1 - eDT) / oneGeDT * (β - D) / (2 * γ)
    return A + Bs * v0 + Bv * v0^2
end

function evaluateLogCharFuncZhu(cf::CharFunc{SchobelZhuParams{T},CR}, z::CT, τ::T)::CR where {T,CR,CT}
    #Zhu formulation: faster but issue with complex log. Here we use Cui log adj.
    p = model(cf)
    v0 = p.v0
    κ = p.κ
    θ = p.θ
    ρ = p.ρ
    σ = p.σ
    cc1 = oneim(cf)
    iu = cc1 * z
    s12 = -iu / 2 * (iu * (1 - ρ^2) - 1 + 2 * ρ * κ / σ)
    s22 = ρ * κ * θ / σ * iu
    s32 = ρ / (2 * σ) * iu
    γ1 = sqrt(2 * σ^2 * s12 + κ^2)
    γ3 = κ^2 * θ - s22 * σ^2
    γ2 = (κ - 2 * σ^2 * s32) / γ1
    ch1 = cosh(γ1 * τ)
    sh1 = sinh(γ1 * τ)
    egt = ch1 + sh1
    γ4 = ch1 + γ2 * sh1
    H3 = (κ - γ1 * (sh1 + γ2 * ch1) / γ4) / σ^2
    H4 = ((κ * θ * γ1 - γ2 * γ3) + γ3 * (sh1 + γ2 * ch1)) / (γ4 * γ1 * σ^2) - κ * θ / σ^2
    #H5 = -(γ1*τ-log(γ1)+ log(γ1* γ4/egt) )/ 2
    H5 = -(γ1 * τ + log(γ4 / egt)) / 2 # here we adjust the log to always use the main branch.
    #original SZ
    #H5 = -log(γ4)/2
    H5 += κ * τ / 2 + (sh1 / γ4 - γ1 * τ) / (2 * γ1^3 * σ^2) * ((κ^2 * θ^2 * γ1^2 - γ3^2))
    H5 += γ3 / (γ1^3 * γ4 * σ^2) * (κ * θ * γ1 - γ2 * γ3) * (ch1 - 1)
    return -s32 * v0^2 - iu * ρ * σ * τ / 2 + H3 / 2 * v0^2 + H4 * v0 + H5
end

function computeCumulants(p::SchobelZhuParams{T}, τ::T) where {T}
    κ = p.κ
    θ = p.θ
    v0 = p.v0
    σ = p.σ
    ρ = p.ρ

    shk = sinh(κ * τ)
    chk = cosh(κ * τ)
    c1H3 = (chk * (chk - shk) * σ + ρ * κ - σ) / (σ * κ)
    c1H4 = (chk - 1) * θ * (shk - chk) / κ
    c1H5 =
        (
            (2 * κ * θ^2 - σ^2) * chk^2 +
            ((-2 * κ * θ^2 + σ^2) * shk - 4 * κ * θ^2) * chk +
            4 * shk * κ * θ^2 +
            2 * τ * (ρ * σ - θ^2) * κ^2 +
            (-σ^2 * τ + 2 * θ^2) * κ +
            σ^2
        ) / (4 * κ^2)
    c1 = -ρ / (2 * σ) * v0^2 - ρ * σ * τ / 2 + c1H3 / 2 * v0^2 + c1H4 * v0 + c1H5
    c2H3 =
        (
            (-2 * κ * ρ * σ + 2 * κ^2) * chk^2 -
            2 * (κ * ρ * σ - 1 / 2 * σ^2 - κ^2) * shk * chk +
            (2 * ρ * σ * τ - 2) * κ^2 +
            (-σ^2 * τ + 2 * ρ * σ) * κ
        ) / (κ^3 * (2 * chk * shk + 2 * chk^2 - 1))
    c2H4 =
        2 *
        (
            (-2 * κ * ρ * σ + κ^2 + σ^2) * chk^2 +
            (
                (-2 * κ * ρ * σ + 1 / 2 * σ^2 + κ^2) * shk +
                (ρ * σ * τ - 1) * κ^2 +
                (-1 / 2 * σ^2 * τ + 2 * ρ * σ) * κ - σ^2
            ) * chk +
            ((ρ * σ * τ - 1) * κ^2 + (-1 / 2 * σ^2 * τ + 2 * ρ * σ) * κ - σ^2 / 2) * shk -
            κ * σ * (ρ * κ - σ / 2) * τ
        ) *
        θ / (κ^3 * (2 * chk * shk + 2 * chk^2 - 1))
    c2H5 =
        (
            (-8 * κ * σ^2 * θ^2 + 2 * σ^4) * chk^4 -
            2 * σ^2 * ((-4 * κ * θ^2 + σ^2) * shk - 8 * κ * θ^2) * chk^3 -
            8 *
            κ *
            (
                2 * σ^2 * θ^2 * shk +
                (-2 * ρ * σ * τ * θ^2 + θ^2) * κ^2 +
                σ * (ρ * τ * σ^2 + (τ * θ^2 - 1 / 2) * σ - 3 * ρ * θ^2) * κ +
                σ^2 * (-1 / 2 * σ^2 * τ + ρ * σ + 2 * θ^2)
            ) *
            chk^2 +
            (
                (
                    (-16 * ρ * σ * τ * θ^2 + 8 * θ^2) * κ^3 +
                    (8 * ρ * τ * σ^3 + (8 * τ * θ^2 - 4) * σ^2 - 24 * ρ * σ * θ^2) * κ^2 +
                    (-4 * σ^4 * τ + 8 * ρ * σ^3 + 20 * σ^2 * θ^2) * κ - σ^4
                ) * shk -
                16 *
                ((ρ * σ * τ - 1) * κ^2 + (-1 / 2 * σ^2 * τ + 3 * ρ * σ) * κ - σ^2) *
                κ *
                θ^2
            ) * chk +
            16 *
            κ *
            θ^2 *
            ((ρ * σ * τ - 1) * κ^2 + (-1 / 2 * σ^2 * τ + 3 * ρ * σ) * κ - (3 * σ^2) / 2) *
            shk +
            8 * κ^4 * τ * θ^2 +
            (-24 * ρ * σ * τ * θ^2 + 4 * σ^2 * τ - 8 * θ^2) * κ^3 +
            ((12 * τ * θ^2 - 4) * σ^2 + 24 * ρ * σ * θ^2) * κ^2 +
            (-σ^4 * τ + 8 * ρ * σ^3 - 8 * σ^2 * θ^2) * κ - 2 * σ^4
        ) / (8 * κ^4)
    c2 = c2H3 / 2 * v0^2 + c2H4 * v0 + c2H5

    c4H3 =
        32 *
        σ^2 *
        (
            (
                (3 * ρ^2 + 3 / 4) * κ^4 - 3 * σ * ρ * (ρ^2 + 3 / 2) * κ^3 +
                9 * (ρ^2 + 1 / 4) * σ^2 * κ^2 / 2 - (15 * κ * ρ * σ^3) / 8 + (3 * σ^4) / 16
            ) * chk^4 -
            3 *
            (
                (-ρ^2 - 1 / 4) * κ^4 + σ * ρ * (ρ^2 + 3 / 2) * κ^3 -
                3 * (ρ^2 + 1 / 4) * σ^2 * κ^2 / 2 + (5 * κ * ρ * σ^3) / 8 - (3 * σ^4) / 32
            ) *
            shk *
            chk^3 +
            (
                ρ^2 * τ^2 * (-3 / 2 + ρ * σ * τ) * κ^6 +
                (3 * τ * (-ρ^2 * σ^2 * τ^2 + ρ * τ * (ρ^2 + 2) * σ - ρ^2 - 1 / 2) * κ^5) /
                2 +
                (
                    (3 * ρ * σ^3 * τ^3) / 4 - 3 * τ^2 * (ρ^2 + 3 / 8) * σ^2 +
                    (3 * ρ * τ * (ρ^2 + 2) * σ) / 2 - 3 / 4 - (15 * ρ^2) / 4
                ) * κ^4 +
                (
                    15 *
                    σ *
                    (
                        -σ^3 * τ^3 / 30 +
                        (2 * ρ * τ^2 * σ^2) / 5 +
                        (-3 / 5 * ρ^2 * τ - 3 / 20 * τ) * σ +
                        ρ^3 +
                        (6 * ρ) / 5
                    ) *
                    κ^3
                ) / 4 -
                (9 * σ^2 * (1 / 24 * σ^2 * τ^2 - 5 / 24 * ρ * σ * τ + ρ^2 + 1 / 4) * κ^2) /
                2 + 15 * (-σ * τ / 10 + ρ) * σ^3 * κ / 8 - (3 * σ^4) / 16
            ) * chk^2 +
            (
                ρ^2 * τ^2 * (-3 / 2 + ρ * σ * τ) * κ^6 +
                (3 * τ * (-ρ^2 * σ^2 * τ^2 + ρ * τ * (ρ^2 + 2) * σ - ρ^2 - 1 / 2) * κ^5) /
                2 +
                (
                    (3 * ρ * σ^3 * τ^3) / 4 - 3 * τ^2 * (ρ^2 + 3 / 8) * σ^2 +
                    (3 * ρ * τ * (ρ^2 + 2) * σ) / 2 - 3 / 8 - (9 * ρ^2) / 4
                ) * κ^4 +
                (
                    9 *
                    σ *
                    (
                        -σ^3 * τ^3 / 18 +
                        (2 * ρ * τ^2 * σ^2) / 3 +
                        (-ρ^2 * τ - 1 / 4 * τ) * σ +
                        ρ^3 +
                        ρ
                    ) *
                    κ^3
                ) / 4 -
                9 * (1 / 12 * σ^2 * τ^2 - 1 / 6 * ρ * σ * τ + ρ^2 + 1 / 8) * σ^2 * κ^2 / 4 +
                (15 * σ^3 * (σ * τ / 5 + ρ) * κ) / 32 - (3 * σ^4) / 64
            ) *
            shk *
            chk -
            (
                (
                    ρ^2 * τ^2 * (-3 / 2 + ρ * σ * τ) * κ^5 +
                    (
                        3 *
                        τ *
                        (-ρ^2 * σ^2 * τ^2 + ρ * τ * (ρ^2 + 2) * σ - ρ^2 - 1 / 2) *
                        κ^4
                    ) / 2 +
                    (
                        (3 * ρ * σ^3 * τ^3) / 4 - (3 * τ^2 * (ρ^2 + 3 / 4) * σ^2) / 2 +
                        (3 * ρ * τ * (ρ^2 + 1) * σ) / 2 - (3 * ρ^2) / 2
                    ) * κ^3 +
                    (-1 / 8 * σ^4 * τ^3 + 3 / 2 * ρ^3 * σ + 9 / 16 * σ^2 * τ) * κ^2 -
                    (15 * σ^3 * (-σ * τ / 5 + ρ) * τ * κ) / 16 + (3 * σ^4 * τ) / 32
                ) * κ
            ) / 2
        ) / (κ^7 * (8 * chk^4 + 8 * chk^3 * shk - 8 * chk^2 - 4 * chk * shk + 1))
    c4H4 =
        -32 *
        σ^2 *
        θ *
        (
            -(3 * chk^8 * σ^4) / 4 +
            (3 * σ^4 * (shk + 1) * chk^7) / 4 +
            (
                9 *
                σ^2 *
                (-σ^2 * shk / 6 + κ * ((ρ * σ * τ - 1 / 2) * κ + σ * (-σ * τ / 2 + ρ))) *
                chk^6
            ) / 2 -
            9 *
            (
                ((ρ * σ * τ - 1 / 2) * κ^2 + σ * (-σ * τ / 2 + ρ) * κ - σ^2 / 12) * shk +
                5 * ((ρ * σ * τ - 3 / 5) * κ - σ^2 * τ / 2 + (6 * ρ * σ) / 5) * κ / 6
            ) *
            σ^2 *
            chk^5 / 2 +
            (
                15 *
                (
                    (ρ * σ * τ - 3 / 5) * κ^2 + (-1 / 2 * σ^2 * τ + 6 / 5 * ρ * σ) * κ -
                    σ^2 / 10
                ) *
                σ^2 *
                shk / 4 -
                6 *
                κ *
                (
                    (ρ^2 * σ^2 * τ^2 - ρ * σ * τ + 1 / 8) * κ^3 +
                    2 * σ * (-ρ * τ^2 * σ^2 / 2 + τ * (ρ^2 + 3 / 8) * σ - ρ / 2) * κ^2 +
                    (1 / 4 * τ^2 * σ^4 - 7 / 8 * ρ * τ * σ^3 + 5 / 4 * ρ^2 * σ^2) * κ -
                    σ^4 * τ / 8
                )
            ) * chk^4 +
            (
                (
                    (6 * ρ^2 * σ^2 * τ^2 - 6 * ρ * σ * τ + 3 / 4) * κ^4 +
                    12 * σ * (-ρ * τ^2 * σ^2 / 2 + τ * (ρ^2 + 3 / 8) * σ - ρ / 2) * κ^3 +
                    (15 * σ^2 * (1 / 5 * σ^2 * τ^2 - ρ * σ * τ + ρ^2 + 3 / 20) * κ^2) / 2 -
                    (9 * σ^3 * (-σ * τ / 6 + ρ) * κ) / 4 + (9 * σ^4) / 32
                ) * shk +
                (
                    27 *
                    κ *
                    (
                        (ρ^2 * σ^2 * τ^2 - 4 / 3 * ρ * σ * τ + 2 / 9) * κ^3 +
                        (
                            8 *
                            σ *
                            (
                                -(3 * ρ * τ^2 * σ^2) / 8 + τ * (ρ^2 + 3 / 8) * σ -
                                (2 * ρ) / 3
                            ) *
                            κ^2
                        ) / 3 +
                        (1 / 4 * τ^2 * σ^4 - 10 / 9 * ρ * τ * σ^3 + 20 / 9 * ρ^2 * σ^2) * κ -
                        (7 * σ^4 * τ) / 36
                    )
                ) / 8
            ) * chk^3 +
            (
                (
                    (-27 / 8 * ρ^2 * σ^2 * τ^2 + 9 / 2 * ρ * σ * τ - 3 / 4) * κ^4 -
                    9 *
                    σ *
                    (-(3 * ρ * τ^2 * σ^2) / 8 + τ * (ρ^2 + 3 / 8) * σ - (2 * ρ) / 3) *
                    κ^3 -
                    (
                        15 *
                        σ^2 *
                        (9 / 80 * σ^2 * τ^2 - 3 / 4 * ρ * σ * τ + ρ^2 + 3 / 20) *
                        κ^2
                    ) / 2 + 9 * (-σ * τ / 8 + ρ) * σ^3 * κ / 4 - (9 * σ^4) / 32
                ) * shk +
                κ * (
                    ρ^2 * τ^2 * (-3 / 2 + ρ * σ * τ) * κ^5 +
                    3 *
                    τ *
                    (-ρ^2 * σ^2 * τ^2 / 2 + ρ * τ * (ρ^2 + 1) * σ - ρ^2 - 1 / 4) *
                    κ^4 +
                    (
                        (3 * ρ * σ^3 * τ^3) / 4 - (9 * σ^2 * τ^2) / 8 +
                        9 * (ρ^2 + 1 / 3) * τ * ρ * σ / 2 - (9 * ρ^2) / 4
                    ) * κ^3 +
                    3 *
                    σ *
                    (
                        -σ^3 * τ^3 / 24 - (7 * ρ * τ^2 * σ^2) / 8 +
                        τ * (ρ^2 / 2 + 9 / 16) * σ +
                        ρ^3
                    ) *
                    κ^2 - 27 * (-(5 * σ * τ) / 18 + ρ) * σ^3 * τ * κ / 8 +
                    (3 * σ^4 * τ) / 8
                )
            ) * chk^2 +
            (
                (
                    (-ρ^3 * σ * τ^3 + 3 / 2 * ρ^2 * τ^2) * κ^6 -
                    3 *
                    τ *
                    (-ρ^2 * σ^2 * τ^2 / 2 + ρ * τ * (ρ^2 + 1) * σ - ρ^2 - 1 / 4) *
                    κ^5 +
                    (
                        -(3 * ρ * σ^3 * τ^3) / 4 + 3 * τ^2 * (ρ^2 + 3 / 8) * σ^2 -
                        (9 * ρ * τ * (ρ^2 + 1) * σ) / 2 +
                        3 / 8 +
                        (9 * ρ^2) / 4
                    ) * κ^4 -
                    3 *
                    (
                        -σ^3 * τ^3 / 24 +
                        ρ * τ^2 * σ^2 / 8 +
                        τ * (-(3 * ρ^2) / 2 - 3 / 16) * σ +
                        ρ^3 +
                        ρ
                    ) *
                    σ *
                    κ^3 +
                    15 *
                    (-1 / 20 * σ^2 * τ^2 - 1 / 4 * ρ * σ * τ + ρ^2 + 9 / 40) *
                    σ^2 *
                    κ^2 / 4 - (27 * σ^3 * (-σ * τ / 18 + ρ) * κ) / 16 + (15 * σ^4) / 64
                ) * shk -
                κ * (
                    (ρ^3 * σ * τ^3 - 3 * ρ^2 * τ^2) * κ^5 +
                    6 *
                    (-ρ^2 * σ^2 * τ^2 / 4 + ρ * τ * (ρ^2 + 1) * σ - 2 * ρ^2 - 1 / 2) *
                    τ *
                    κ^4 +
                    (
                        (3 * ρ * σ^3 * τ^3) / 4 +
                        τ^2 * (9 * ρ^2 - 9 / 4) * σ^2 +
                        18 * ρ^3 * σ * τ - 18 * ρ^2
                    ) * κ^3 +
                    (
                        -σ^4 * τ^3 / 8 - (57 * ρ * τ^2 * σ^3) / 4 +
                        τ * (18 * ρ^2 + 45 / 4) * σ^2 +
                        24 * ρ^3 * σ
                    ) * κ^2 +
                    (33 / 8 * τ^2 * σ^4 - 45 / 2 * ρ * τ * σ^3) * κ +
                    (21 * σ^4 * τ) / 8
                ) / 8
            ) * chk +
            (
                ρ^2 * τ^2 * (ρ * σ * τ - 3) * κ^6 / 8 +
                3 *
                (-ρ^2 * σ^2 * τ^2 / 4 + ρ * τ * (ρ^2 + 1) * σ - 2 * ρ^2 - 1 / 2) *
                τ *
                κ^5 / 4 +
                (
                    (3 * ρ * σ^3 * τ^3) / 32 - 9 * (ρ^2 + 1 / 2) * τ^2 * σ^2 / 16 +
                    (9 * ρ * τ * (ρ^2 + 1) * σ) / 4 - 3 / 8 - (9 * ρ^2) / 4
                ) * κ^4 +
                3 *
                σ *
                (
                    -σ^3 * τ^3 / 192 - ρ * τ^2 * σ^2 / 32 +
                    τ * (-(3 * ρ^2) / 4 - 3 / 32) * σ +
                    ρ^3 +
                    ρ
                ) *
                κ^3 -
                (
                    15 *
                    σ^2 *
                    (-1 / 40 * σ^2 * τ^2 - 1 / 8 * ρ * σ * τ + ρ^2 + 9 / 40) *
                    κ^2
                ) / 4 + 27 * (-σ * τ / 36 + ρ) * σ^3 * κ / 16 - (15 * σ^4) / 64
            ) * shk -
            κ *
            (
                ρ^2 * τ * (-3 / 2 + ρ * σ * τ) * κ^5 +
                (-(3 * ρ^2 * σ^2 * τ^2) / 2 + 3 * ρ * τ * (ρ^2 + 1) * σ - 3 * ρ^2 - 3 / 4) *
                κ^4 +
                (
                    9 *
                    σ *
                    (ρ * τ^2 * σ^2 / 6 + τ * (-ρ^2 - 1 / 4) * σ + ρ^3 + (4 * ρ) / 3) *
                    κ^3
                ) / 2 -
                15 * (1 / 60 * σ^2 * τ^2 - 1 / 4 * ρ * σ * τ + ρ^2 + 9 / 40) * σ^2 * κ^2 /
                2 + (27 * σ^3 * (-σ * τ / 18 + ρ) * κ) / 8 - (15 * σ^4) / 32
            ) *
            τ / 2
        ) / κ^7
    c4H5 =
        -8 *
        σ^2 *
        (
            -(3 * σ^4 * (-8 * κ * θ^2 + σ^2) * chk^8) / 16 +
            (3 * σ^4 * ((-8 * κ * θ^2 + σ^2) * shk - 16 * κ * θ^2) * chk^7) / 16 +
            (
                3 *
                κ *
                σ^2 *
                (
                    2 * σ^2 * θ^2 * shk +
                    (-6 * ρ * σ * τ * θ^2 + 3 * θ^2) * κ^2 +
                    σ * (ρ * τ * σ^2 + (3 * τ * θ^2 - 1 / 2) * σ - 7 * ρ * θ^2) * κ +
                    σ^2 * (-1 / 2 * σ^2 * τ + ρ * σ + 2 * θ^2)
                ) *
                chk^6
            ) / 2 -
            (
                3 *
                σ^2 *
                (
                    (
                        (-6 * ρ * σ * τ * θ^2 + 3 * θ^2) * κ^3 +
                        σ * (ρ * τ * σ^2 + (3 * τ * θ^2 - 1 / 2) * σ - 7 * ρ * θ^2) * κ^2 +
                        σ^2 * (-1 / 2 * σ^2 * τ + ρ * σ + 5 / 2 * θ^2) * κ - σ^4 / 16
                    ) * shk -
                    10 *
                    (
                        (ρ * σ * τ - 3 / 5) * κ^2 + (-1 / 2 * σ^2 * τ + 7 / 5 * ρ * σ) * κ -
                        σ^2 / 5
                    ) *
                    κ *
                    θ^2
                ) *
                chk^5
            ) / 2 -
            3 *
            κ *
            (
                5 *
                σ^2 *
                (
                    (ρ * σ * τ - 3 / 5) * κ^2 + (-1 / 2 * σ^2 * τ + 7 / 5 * ρ * σ) * κ -
                    (3 * σ^2) / 10
                ) *
                θ^2 *
                shk - 4 * (ρ^2 * σ^2 * τ^2 - ρ * σ * τ + 1 / 8) * θ^2 * κ^4 +
                σ *
                (
                    ρ^2 * σ^3 * τ^2 +
                    (4 * ρ * τ^2 * θ^2 - ρ * τ) * σ^2 +
                    (1 / 8 + (-10 * ρ^2 - 3) * τ * θ^2) * σ +
                    5 * ρ * θ^2
                ) *
                κ^3 +
                2 *
                (
                    -ρ * τ^2 * σ^3 / 2 +
                    τ * (-τ * θ^2 / 2 + ρ^2 + 3 / 8) * σ^2 +
                    (17 / 4 * ρ * τ * θ^2 - 1 / 2 * ρ) * σ - 15 * (ρ^2 + 2 / 5) * θ^2 / 4
                ) *
                σ^2 *
                κ^2 +
                (
                    σ^6 * τ^2 / 4 - ρ * σ^5 * τ +
                    (-(3 * τ * θ^2) / 2 + (5 * ρ^2) / 4) * σ^4 +
                    7 * ρ * σ^3 * θ^2
                ) * κ - θ^2 * σ^4 - σ^6 * τ / 16
            ) *
            chk^4 +
            (
                (
                    -12 * (ρ^2 * σ^2 * τ^2 - ρ * σ * τ + 1 / 8) * θ^2 * κ^5 +
                    3 *
                    σ *
                    (
                        ρ^2 * σ^3 * τ^2 +
                        (4 * ρ * τ^2 * θ^2 - ρ * τ) * σ^2 +
                        (1 / 8 + (-10 * ρ^2 - 3) * τ * θ^2) * σ +
                        5 * ρ * θ^2
                    ) *
                    κ^4 +
                    6 *
                    σ^2 *
                    (
                        -ρ * τ^2 * σ^3 / 2 +
                        τ * (-τ * θ^2 / 2 + ρ^2 + 3 / 8) * σ^2 +
                        (5 * ρ * τ * θ^2 - 1 / 2 * ρ) * σ - 15 * (ρ^2 + 1 / 2) * θ^2 / 4
                    ) *
                    κ^3 +
                    15 *
                    (
                        σ^3 * τ^2 / 5 - ρ * τ * σ^2 +
                        (-(9 * τ * θ^2) / 5 + ρ^2 + 1 / 10) * σ +
                        7 * ρ * θ^2
                    ) *
                    σ^3 *
                    κ^2 / 4 -
                    (3 * σ^4 * (-1 / 4 * σ^2 * τ + ρ * σ + 27 / 4 * θ^2) * κ) / 4 +
                    (9 * σ^6) / 128
                ) * shk -
                (
                    27 *
                    κ *
                    (
                        (ρ^2 * σ^2 * τ^2 - 4 / 3 * ρ * σ * τ + 2 / 9) * κ^4 +
                        (
                            -ρ * τ^2 * σ^3 + τ * ((10 * ρ^2) / 3 + 1) * σ^2 -
                            (20 * ρ * σ) / 9
                        ) * κ^3 +
                        (
                            10 *
                            σ^2 *
                            (3 / 40 * σ^2 * τ^2 - 19 / 30 * ρ * σ * τ + ρ^2 + 1 / 5) *
                            κ^2
                        ) / 3 - (14 * σ^3 * (-(5 * σ * τ) / 56 + ρ) * κ) / 9 +
                        (2 * σ^4) / 9
                    ) *
                    θ^2
                ) / 2
            ) * chk^3 +
            κ *
            (
                (
                    27 *
                    θ^2 *
                    (
                        (ρ^2 * σ^2 * τ^2 - 4 / 3 * ρ * σ * τ + 2 / 9) * κ^4 +
                        (
                            -ρ * τ^2 * σ^3 + τ * ((10 * ρ^2) / 3 + 1) * σ^2 -
                            (20 * ρ * σ) / 9
                        ) * κ^3 +
                        (
                            τ^2 * σ^4 / 4 - (8 * ρ * τ * σ^3) / 3 +
                            ((10 * ρ^2) / 3 + 1) * σ^2
                        ) * κ^2 - (7 * σ^3 * (-(5 * σ * τ) / 28 + ρ) * κ) / 3 +
                        (5 * σ^4) / 12
                    ) *
                    shk
                ) / 2 +
                (-2 * ρ^3 * σ * τ^3 * θ^2 + 3 * ρ^2 * τ^2 * θ^2) * κ^6 +
                (
                    ρ^3 * σ^3 * τ^2 + 3 * τ * ρ^2 * (τ * θ^2 - 1 / 2) * σ^2 -
                    9 * (ρ^2 + 2 / 3) * τ * ρ * θ^2 * σ + 9 * (ρ^2 + 1 / 6) * θ^2
                ) *
                τ *
                κ^5 +
                (
                    -(3 * ρ^2 * σ^4 * τ^3) / 2 +
                    3 * (-τ * θ^2 / 2 + ρ^2 + 1) * τ^2 * ρ * σ^3 +
                    9 * (ρ^2 + 1 / 4) * τ * (τ * θ^2 - 1 / 3) * σ^2 -
                    18 * ρ * τ * θ^2 * (ρ^2 + 1) * σ + (9 * ρ^2 + 3) * θ^2
                ) * κ^4 +
                (
                    9 *
                    σ *
                    (
                        ρ * σ^4 * τ^3 / 6 - ((-τ * θ^2 / 9 + ρ^2 + 1 / 2) * τ^2 * σ^3) / 2 +
                        ρ * τ * (-τ * θ^2 / 3 + ρ^2 + 2 / 3) * σ^2 +
                        (τ * ((20 * ρ^2) / 3 + 5 / 4) * θ^2 - ρ^2 / 2) * σ -
                        (10 * ρ * θ^2 * (ρ^2 + 2)) / 3
                    ) *
                    κ^3
                ) / 2 +
                3 *
                σ^2 *
                (
                    -σ^4 * τ^3 / 24 - ρ * τ^2 * σ^3 / 8 -
                    ((τ * θ^2 / 4 + ρ^2 - 3 / 8) * τ * σ^2) / 2 +
                    (-27 / 8 * ρ * τ * θ^2 + ρ^3) * σ +
                    (15 * ρ^2 + 3) * θ^2
                ) *
                κ^2 +
                (-21 * ρ * σ^3 * θ^2 - 9 / 8 * ρ * σ^5 * τ + 3 / 8 * σ^6 * τ^2) * κ +
                3 * θ^2 * σ^4 +
                (3 * σ^6 * τ) / 32
            ) *
            chk^2 +
            (
                (
                    (2 * ρ^3 * σ * τ^3 * θ^2 - 3 * ρ^2 * τ^2 * θ^2) * κ^7 -
                    (
                        ρ^3 * σ^3 * τ^2 + 3 * τ * ρ^2 * (τ * θ^2 - 1 / 2) * σ^2 -
                        9 * (ρ^2 + 2 / 3) * τ * ρ * θ^2 * σ + 9 * (ρ^2 + 1 / 6) * θ^2
                    ) *
                    τ *
                    κ^6 +
                    (
                        (3 * ρ^2 * σ^4 * τ^3) / 2 -
                        3 * (-τ * θ^2 / 2 + ρ^2 + 1) * τ^2 * ρ * σ^3 +
                        (-15 * τ^2 * (ρ^2 + 3 / 20) * θ^2 + τ * (3 * ρ^2 + 3 / 4)) * σ^2 +
                        18 * (ρ^2 + 4 / 3) * τ * ρ * θ^2 * σ - 9 * (ρ^2 + 5 / 12) * θ^2
                    ) * κ^5 -
                    (
                        9 *
                        σ *
                        (
                            ρ * σ^4 * τ^3 / 6 -
                            (5 * τ^2 * (-τ * θ^2 / 15 + ρ^2 + 3 / 10) * σ^3) / 6 +
                            ρ * τ * (-(5 * τ * θ^2) / 3 + ρ^2 + 1) * σ^2 +
                            (τ * (10 * ρ^2 + 9 / 4) * θ^2 - ρ^2 / 2 - 1 / 24) * σ -
                            10 * (ρ^2 + 5 / 2) * ρ * θ^2 / 3
                        ) *
                        κ^4
                    ) / 2 -
                    3 *
                    σ^2 *
                    (
                        -σ^4 * τ^3 / 24 + (3 * ρ * τ^2 * σ^3) / 8 -
                        3 * (-τ * θ^2 / 4 + ρ^2 + 1 / 8) * τ * σ^2 / 2 +
                        ρ * (-(35 * τ * θ^2) / 4 + ρ^2 + 1 / 2) * σ +
                        75 * (ρ^2 + 27 / 100) * θ^2 / 4
                    ) *
                    κ^3 +
                    (
                        15 *
                        σ^3 *
                        (
                            -ρ * τ * σ^2 / 2 +
                            (-(21 * τ * θ^2) / 10 + ρ^2 + 3 / 20) * σ +
                            (189 * ρ * θ^2) / 10
                        ) *
                        κ^2
                    ) / 8 - (9 * σ^4 * (-1 / 6 * σ^2 * τ + ρ * σ + 65 / 6 * θ^2) * κ) / 16 +
                    (15 * σ^6) / 256
                ) * shk +
                κ *
                (
                    (ρ^3 * σ * τ^3 - 3 * ρ^2 * τ^2) * κ^6 +
                    9 *
                    (-ρ^2 * σ^2 * τ^2 / 6 + ρ * τ * (ρ^2 + 2 / 3) * σ - 2 * ρ^2 - 1 / 3) *
                    τ *
                    κ^5 +
                    (
                        (3 * ρ * σ^3 * τ^3) / 4 +
                        3 * τ^2 * (ρ^2 - 3 / 4) * σ^2 +
                        (36 * ρ^3 + 18 * ρ) * τ * σ - 36 * ρ^2 - 6
                    ) * κ^4 +
                    60 *
                    σ *
                    (
                        -σ^3 * τ^3 / 480 - (7 * ρ * τ^2 * σ^2) / 40 +
                        τ * (-ρ^2 / 4 + 3 / 80) * σ +
                        ρ^3 +
                        ρ
                    ) *
                    κ^3 -
                    90 *
                    σ^2 *
                    (-3 / 80 * σ^2 * τ^2 + 1 / 8 * ρ * σ * τ + ρ^2 + 1 / 5) *
                    κ^2 + (27 / 8 * σ^4 * τ + 42 * ρ * σ^3) * κ - 6 * σ^4
                ) *
                θ^2 / 2
            ) * chk -
            κ *
            (
                (ρ^3 * σ * τ^3 - 3 * ρ^2 * τ^2) * κ^6 +
                9 *
                (-ρ^2 * σ^2 * τ^2 / 6 + ρ * τ * (ρ^2 + 2 / 3) * σ - 2 * ρ^2 - 1 / 3) *
                τ *
                κ^5 +
                (
                    (3 * ρ * σ^3 * τ^3) / 4 - 21 * (ρ^2 + 3 / 14) * τ^2 * σ^2 / 2 +
                    36 * ρ * τ * (ρ^2 + 1) * σ - 36 * ρ^2 - 9
                ) * κ^4 +
                60 *
                (
                    -σ^3 * τ^3 / 480 +
                    ρ * τ^2 * σ^2 / 20 +
                    τ * (-ρ^2 - 3 / 16) * σ +
                    ρ^3 +
                    (3 * ρ) / 2
                ) *
                σ *
                κ^3 - 135 * σ^2 * (-19 / 90 * ρ * σ * τ + ρ^2 + 1 / 4) * κ^2 +
                (-33 / 8 * σ^4 * τ + 315 / 4 * ρ * σ^3) * κ - (105 * σ^4) / 8
            ) *
            θ^2 *
            shk / 2 + (-3 / 2 + ρ * σ * τ) * τ^2 * ρ^2 * θ^2 * κ^7 -
            (
                (
                    ρ^3 * σ^3 * τ^2 + 3 * τ * ρ^2 * (τ * θ^2 - 1 / 2) * σ^2 -
                    9 * (ρ^2 + 2 / 3) * τ * ρ * θ^2 * σ + 18 * (ρ^2 + 1 / 4) * θ^2
                ) *
                τ *
                κ^6
            ) / 2 +
            (
                (3 * ρ^2 * σ^4 * τ^3) / 4 -
                3 * (-τ * θ^2 / 2 + ρ^2 + 1) * τ^2 * ρ * σ^3 / 2 -
                9 * (τ * (ρ^2 + 1 / 8) * θ^2 - ρ^2 / 12 - 1 / 48) * τ * σ^2 +
                15 * (ρ^2 + 17 / 10) * τ * ρ * θ^2 * σ +
                9 * (ρ^2 + 1 / 6) * θ^2
            ) * κ^5 -
            (
                3 *
                σ *
                (
                    ρ * σ^4 * τ^3 / 4 - 3 * (-τ * θ^2 / 18 + ρ^2 + 1 / 4) * τ^2 * σ^3 / 2 +
                    ρ * τ * (-(7 * τ * θ^2) / 2 + ρ^2 + 1) * σ^2 +
                    (τ * ((55 * ρ^2) / 2 + 57 / 8) * θ^2 - (3 * ρ^2) / 2 - 1 / 4) * σ +
                    10 * ρ * θ^2 * (ρ^2 + 1)
                ) *
                κ^4
            ) / 2 -
            3 *
            (
                -σ^4 * τ^3 / 48 + (5 * ρ * τ^2 * σ^3) / 16 -
                5 * (-τ * θ^2 / 2 + ρ^2 + 3 / 20) * τ * σ^2 / 8 +
                ρ * (-(141 * τ * θ^2) / 16 + ρ^2 + 1) * σ - 15 * (ρ^2 + 1 / 5) * θ^2 / 2
            ) *
            σ^2 *
            κ^3 +
            (
                15 *
                σ^3 *
                (
                    σ^3 * τ^2 / 40 - (3 * ρ * τ * σ^2) / 20 +
                    (-(49 * τ * θ^2) / 40 + ρ^2 + 1 / 5) * σ - (14 * ρ * θ^2) / 5
                ) *
                κ^2
            ) / 4 - (3 * σ^4 * (-5 / 128 * σ^2 * τ + ρ * σ - θ^2) * κ) / 2 + (3 * σ^6) / 16
        ) / κ^8
    c4 = c4H3 / 2 * v0^2 + c4H4 * v0 + c4H5
    return c1, c2, c4
end


text:
using CodeInfoTools
import Ghost: Tape, Variable, V, Input, Constant, Call, mkcall
import Core: CodeInfo, SSAValue, SlotNumber

mutable struct Frame
    ir2tape::Dict{Any, Variable}
end

function Frame(ci::CodeInfo, v_fargs...)
    ir2tape = Dict(SlotNumber(i) => v_fargs[i] for i=1:length(v_fargs) if v_fargs[i] isa V)
    return Frame(ir2tape)
end

function Base.show(io::IO, frame::Frame)
    s = "Frame(\n"
    for (sv, v) in sort(frame.ir2tape, by=sv -> sv.id)
        s *= "  $sv => $v\n"
    end
    s *= ")"
    print(io, s)
end

function resolve_tape_vars(frame::Frame, sv_fargs...)
    v_fargs = []
    for sv in sv_fargs
        if sv isa SlotNumber || sv isa SSAValue
            push!(v_fargs, frame.ir2tape[sv])
        elseif sv isa GlobalRef
            push!(v_fargs, getfield(sv.mod, sv.name))
        else
            # treat as constant value
            push!(v_fargs, sv)
        end
    end
    return v_fargs
end

struct Fragment
    ci::CodeInfo
    from::Int
    to::Int
end

struct IfSpec
    cond::SSAValue
    true_block::Fragment
    false_block::Fragment
end

mutable struct Tracer
    tape::Tape
    frames::Vector{Frame}  # TODO: remove, we don't need frames field here
    primitives   # TODO: make it more Ghost-like instead
end

"""
    push_call!(tape::Tape, fn, args...; kwargs)

Shortcut for `push!(tape, mkcall(fn, args..))` also handling
keyword arguments and respecting `ONNXCtx.exec` setting.
"""
function push_call!(tape::Tape, fn, args...; kwargs...)
    kwargs = NamedTuple(kwargs)
    if !isempty(kwargs)
        args = (kwargs, fn, args...)
        fn = Core.kwfunc(fn)
    end
    # op = tape.c.exec ? mkcall(fn, args...) : mkcall(fn, args...; val=nothing)
    op = mkcall(fn, args...)
    return push!(tape, op)
end

# TODO: imagine we have complete IfSpecs and LoopSpecs, so what? how are we going to use them?

function detect_ifs(ci::CodeInfo)
    ifs = []
    for (i, st) in enumerate(ci.code)
        if st isa Core.GotoIfNot
            j = i + 1
            while j < length(ci.code) && (ci.code[j] isa Core.GotoIfNot || ci.code[j] isa Core.ReturnNode)
                j += 1
            end
            if_spec = IfSpec(st.cond, Fragment(ci, st.cond.id, j), nothing)
        end
    end

end

function trace!(t::Tracer, ci::CodeInfo, v_fargs...)
    frame = Frame(ci, v_fargs...)
    i = 1
    while i <= length(ci.code)
        st = ci.code[i]
        @show i, st
        if Meta.isexpr(st, :call) || (Meta.isexpr(st, :(=)) && Meta.isexpr(st.args[2], :call))
            # function call
            sv = SSAValue(i)
            ex = Meta.isexpr(st, :(=)) ? st.args[2] : st
            vs = resolve_tape_vars(frame, ex.args...)
            fvals = [v isa V ? t.tape[v].val : v for v in vs]
            v = if fvals[1] in t.primitives
                try  # todo: turn into a macro
                    push_call!(t.tape, vs...)
                catch
                    global STATE = (t, ci, v_fargs, frame, i, st)
                    rethrow()
                end
            else
                trace!(t, code_lowered(fvals[1], fvals[2:end])[1], vs...)
            end
            frame.ir2tape[sv] = v
            if Meta.isexpr(st, :(=))
                # update mapping for slot
                slot = st.args[1]
                frame.ir2tape[slot] = v
            end
            i += 1
        elseif Meta.isexpr(st, :(=))
            # constant
            sv = st.args[1]
            v_rhs = resolve_tape_vars(frame, st.args[2])[1]
            v = if v_rhs isa V
                val = t.tape[v_rhs].val
                push!(t.tape, Constant(val))
            else
                # there's no assignment operator, so we record a dummy call to identity
                # would it be safe to simply update mapping ir2tape[sv] = v_rhs?
                push_call!(t.tape, identity, v_rhs)
            end
            frame.ir2tape[sv] = v
            i += 1
        elseif st isa Core.GotoIfNot
            # conditional jump
            v_cond = frame.ir2tape[st.cond]
            # if not cond, set i to destination, otherwise step forward
            i = !t.tape[v_cond].val ? st.dest : i + 1
        elseif st isa Core.GotoNode
            # unconditional jump
            i = st.label
        elseif st isa Core.ReturnNode
            # return statement
            sv = st.val
            return frame.ir2tape[sv]
        else
            error("Unexpected statement type in CodeInfo: $st")
        end
    end
    # for now assume the returned value is always the last recorded
    # we will update this when control flow instructions are implemented
    return V(t.tape[V(end)])
end

"""
Same as trace!(), but tries to parse control flow constructs as well
"""
function trace_cf!(t::Tracer, ci::CodeInfo, v_fargs...)
    # TODO
end

function trace(f, args...)
    ci = code_lowered(f, args)[1]
    t = Tracer(Tape(), [], Set([+, -, *, /, <, <=, >, >=]))
    v_fn = push!(t.tape, Input(f))
    v_args = [push!(t.tape, Input(a)) for a in args]
    v_fargs = [v_fn, v_args...]
    rv = trace!(t, ci, v_fargs...)
    t.tape.result = rv
    return t.tape
end

# v_xxx - Variable xxx
# sv_xxx - SSAValue xxx

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

function main()
    f = linear
    args = [1.0, 2.0, 3.0]
    trace(f, args...)

    f = while1
    args = (2.0,)
    tp = trace(f, args...)

    ci = @code_lowered cond1(2.0)
end

mul(x, y) = x * y
linear(x, a, b) = mul(a, x) + b

pow(x, n) = for i=1:n x = x * x end

function const1(x)
    y = 1
    return x + y
end

function cond1(x)
    y = 2x
    if x > 0
        y = 3x
    end
    return y
end

function while1(x)
    y = 10
    while x > 5
        y -= 1
        x /= 2
    end
    return y
end

function while_continue(x)
    y = 3x
    while y > 0
        if y < x
            continue
        end
        y -= 1
    end
    return y
end

function while_break(x)
    y = 3x
    while y > 0
        if y < x
            break
        end
        y -= 1
    end
    return y
end


text:
import Pkg
Pkg.activate(@__DIR__)

using RipsererBenchmarks

time_start = time_ns()

#RipserComparison.run_and_save_all()
#RipserComparison.grind_all()
RipsererBenchmarks.EireneComparison.run_and_save_all()
Involuted.run_and_save_all()

elapsed = (time_ns() - time_start)/1e9
@info "Time $(ProgressMeter.durationstring(elapsed))"


text:
test_suite_clone_path = joinpath(homedir(), "sbml-test-suite")
if !isdir(test_suite_clone_path)
    run(`git clone https://github.com/sbmlteam/sbml-test-suite $(test_suite_clone_path)`)
end
path = joinpath(test_suite_clone_path, "cases")
@test isdir(path)

function extract_xmls(path)
    dirs = filter(isdir, readdir(path; join=true))
    all_dirs = mapreduce(x -> filter(isdir, readdir(x; join=true)), vcat, dirs)
    all_files = vcat(readdir.(all_dirs; join=true)...)
    filter!(x -> occursin(".xml", x), all_files)
end

# theres probably an idiomatic way to do this
function call_save_err(f, x)
    n = length(x)
    arr = Vector(undef, n)
    for i in 1:n
        try
            arr[i] = f(x[i])
        catch e
            arr[i] = e
        end
    end
    arr
end

# theres probably an idiomatic way to do this
function filtermap(f, x)
    arr = call_save_err(f, x)
    arr[.!(@. typeof(arr) <: Exception)]
end

all_files = extract_xmls(path)
@test all_files isa Vector{String}

n = length(all_files)
@test n > 0

# speed up CI
TEST_ON = 50
all_files = all_files[1:TEST_ON]

models = getmodel.(all_files);
ps = getparameters.(models);
u0s = call_save_err(getinitialconditions, models);
rxs = call_save_err(getreactions, models);
systems = call_save_err(sbml2odesystem, all_files);
probs = call_save_err(x->sbml2odeproblem(x, (0., 1.)), all_files);

errored_mask = @. typeof(systems) <: Exception
errs = systems[errored_mask, :]
good = systems[.!(errored_mask), :]

good_probs = probs[.!(@. typeof(probs) <: Exception)] # eventually benchmark solve on these

@time mtks = filtermap(modelingtoolkitize, good_probs);
eqs = unique(equations.(mtks))
@test length(eqs) > 0

@show eqs


text:
"""
$(SIGNATURES)

Solve the adiabatic master equation defined by `A` for a total evolution time `tf`.

...
# Arguments
- `A::Annealing`: the `Annealing`/`Evolution` object.
- `tf::Real`: the total annealing time.
- `tspan = (0, tf)`: time interval to solve the dynamics.
- `ω_hint=[]`: specify a grid to precompute the ``S`` function in the Lamb shift term. Skip the precomputation if empty.
- `lambshift::Bool=true`: whether to include the Lamb shift in the simulation.
- `lambshift_S=nothing`: provide a custom routine to calculate the `S` function in the Lamb shift term.
- `lvl::Int=size(A.H, 1)`: number of levels to keep. The higher levels are ignored to speed up the computation.
- `vectorize::Bool = false`: whether to vectorize the density matrix.
- `one_sided=false`: whether to solve the one-sided AME.
- `kwargs` : other keyword arguments supported by `DifferentialEquations`.
...
"""
function solve_ame(
    A::Annealing,
    tf::Real;
    tspan=(0.0, tf),
    ω_hint=[],
    lambshift::Bool=true,
    lambshift_S=nothing,
    lvl::Int=size(A.H, 1),
    vectorize::Bool=false,
    one_sided=false,
    kwargs...,
)
    u0 = build_u0(A.u0, :m, vectorize=vectorize)
    if one_sided == false
        L = OpenQuantumBase.davies_from_interactions(A.interactions, ω_hint, lambshift, lambshift_S)
    else
        L = build_onesided_ame(A.interactions, ω_hint, lambshift)
    end
    if vectorize
        error("Vectorization is not yet supported for adiabatic master equation.")
    end
    f = DiffEqLiouvillian(A.H, L, [], lvl)
    p = ODEParams(f, float(tf), A.annealing_parameter)
    prob = ODEProblem(f, u0, tspan, p)
    solve(prob; alg_hints=[:nonstiff], kwargs...)
end

function build_ensemble_ame(
    A::Annealing,
    tf::Real,
    output_func,
    reduction;
    tspan=(0.0, tf),
    ω_hint=[],
    lambshift::Bool=true,
    lambshift_S=nothing,
    lvl::Int=size(A.H, 1),
    initializer=initializer,
    kwargs...,
)
    u0 = build_u0(A.u0, :v)

    flist = OpenQuantumBase.fluctuator_from_interactions(A.interactions)
    dlist = OpenQuantumBase.davies_from_interactions(A.interactions, ω_hint, lambshift, lambshift_S)
    L = DiffEqLiouvillian(A.H, dlist, flist, lvl)
    cb = build_jump_callback(dlist, flist, initializer)
    p = ODEParams(L, tf, A.annealing_parameter)
    update_func = function (cache, u, p, t)
        update_cache!(cache, p.L, p, t)
    end
    cache = get_cache(A.H)
    diff_op = DiffEqArrayOperator(cache, update_func=update_func)
    jac_cache = similar(cache)
    jac_op = DiffEqArrayOperator(jac_cache, update_func=update_func)
    ff = ODEFunction(diff_op; jac_prototype=jac_op)

    prob = ODEProblem{true}(ff, u0, tspan, p, callback=cb)
    ensemble_prob = EnsembleProblem(
        prob;
        output_func=output_func,
        reduction=reduction,
        kwargs...
    )
    ensemble_prob
end

function build_jump_callback(dlist, flist, initializer)
    initializer =
        initializer == DEFAULT_INITIALIZER ? (x, y) -> rand([-1, 1], x, y) :
        initializer
    if isempty(flist)
        cb = LindbladJumpCallback()
    elseif isempty(dlist)
        error("No interactions support AME. Use other ensemble instead.")
    else
        cb = CallbackSet(
            LindbladJumpCallback(),
            [FluctuatorCallback(f, initializer) for f in flist]...,
        )
    end
    cb
end

max_stars_repo_name: huangjingzheng/OpenQuantumTools

max_issues_repo_name: huangjingzheng/OpenQuantumTools

max_forks_repo_name: huangjingzheng/OpenQuantumTools


text:
"""
    A = JopRoughness(sp, dim, w)

`A`  can be used to regularize an optimization problem, applying a penalty to
models that are non-smooth.  This is similar to a finite difference operator,
but here no care is taken to ensure that `A` computes a derivative.  `dim` is
the dimension that one wishes to smooth and `w` is the half-width of the smoother.
The width of the window determines how strong the penalty is for being non-smooth.

For example, the form of `A` for `w=1` is,
```julia
A=
[
(1/3-2/3) 1/3       0         0         ;
1/3       (1/3-3/3) 1/3       0         ;
0         1/3       (1/3-3/3) 1/3       ;
0         0         1/3       (1/3-2/3)
]
```
"""
JopRoughness(sp::JetSpace, dim::Int; w::Int=1) = JopLn(dom = sp, rng = sp, df! = JopRoughness_df!, s = (dim=dim, w=w))

export JopRoughness

function JopRoughness_df!(d, m; dim, w, kwargs...)
    Rpre = CartesianIndices(size(m)[1:dim-1])
    Rpost = CartesianIndices(size(m)[dim+1:end])
    n = size(m, dim)
    roughness_helper(d, m, w, n, Rpre, Rpost)
end

function roughness_helper(d::AbstractArray{T}, m::AbstractArray{T}, w, n, Rpre, Rpost) where {T}
    d .= 0
    _w = 2*w + 1
    for Ipost in Rpost
        for i = 1:n
            f = -one(T)
            for j = 1:_w
                k = i + j - 1 - w
                if k ∈ 1:n
                    for Ipre in Rpre
                        d[Ipre, i, Ipost] += m[Ipre, k, Ipost]
                    end
                else
                    f += one(T)/_w
                end
            end
            for Ipre in Rpre
                d[Ipre, i, Ipost] /= _w
                d[Ipre, i, Ipost] += m[Ipre, i, Ipost]*f
            end
        end
    end
    d
end