text:
using Test
using NoisyQuantumComputerSimulator
sim = NoisyQuantumComputerSimulator

⊗ = kron

matrix_I = Matrix{Complex{Float64}}([1 0; 0 1])
matrix_H = Matrix{Complex{Float64}}(1 / √2 * [1 1; 1 -1])
matrix_X = Matrix{Complex{Float64}}([0 1; 1 0])
matrix_S = Matrix{Complex{Float64}}([1 0; 0 complex(0, 1)])

# Testing Circuit.

c = Circuit(1)
@test c.density_matrix == [1 0; 0 0]
c = Circuit(2)
@test c.density_matrix == [1 0 0 0; 0 0 0 0; 0 0 0 0; 0 0 0 0]
c = Circuit(3)
@test c.density_matrix == [1 0 0 0 0 0 0 0;
                           0 0 0 0 0 0 0 0;
                           0 0 0 0 0 0 0 0;
                           0 0 0 0 0 0 0 0;
                           0 0 0 0 0 0 0 0;
                           0 0 0 0 0 0 0 0;
                           0 0 0 0 0 0 0 0;
                           0 0 0 0 0 0 0 0]

# Testing reset_state!
c = Circuit(1, H(0))
@test exec(c) ≈ [.5 .5; .5 .5]
reset_state!(c)
@test c.density_matrix == [1 0; 0 0]

# Testing CNOT.

BELL = 1 / √2 * [1 0 1 0; 0 1 0 1; 0 1 0 -1; 1 0 -1 0] # (CX)(H⊗I)
BELL₀ = BELL * [1; 0; 0; 0]
BELL₁ = BELL * [0; 1; 0; 0]
BELL₂ = BELL * [0; 0; 1; 0]
BELL₃ = BELL * [0; 0; 0; 1]

@test exec(Circuit(2, H(1), CNOT(1, 0))) == BELL₀ * BELL₀'
@test exec(Circuit(2, X(0), H(1), CNOT(1, 0))) == BELL₁ * BELL₁'
@test exec(Circuit(2, X(1), H(1), CNOT(1, 0))) == BELL₂ * BELL₂'
@test exec(Circuit(2, X(0), X(1), H(1), CNOT(1, 0))) == BELL₃ * BELL₃'
@test exec(Circuit(4, H(2), CNOT(2, 1))) ≈ ([1; 0] ⊗ BELL₀ ⊗ [1; 0]) * ([1 0] ⊗ BELL₀' ⊗ [1 0])

# Testing CCNOT.

CCBELL = [.5  .5   0   0  .5  .5   0   0;
          .5 -.5   0   0  .5 -.5   0   0;
           0   0  .5  .5   0   0  .5  .5;
           0   0  .5 -.5   0   0  .5 -.5;
          .5  .5   0   0 -.5 -.5   0   0;
           0   0  .5 -.5   0   0 -.5  .5;
           0   0  .5  .5   0   0 -.5 -.5;
          .5 -.5   0   0 -.5  .5   0   0]
CCBELL₀ = CCBELL * [1; 0; 0; 0; 0; 0; 0; 0]
CCBELL₁ = CCBELL * [0; 1; 0; 0; 0; 0; 0; 0]
CCBELL₂ = CCBELL * [0; 0; 1; 0; 0; 0; 0; 0]
CCBELL₃ = CCBELL * [0; 0; 0; 1; 0; 0; 0; 0]
CCBELL₄ = CCBELL * [0; 0; 0; 0; 1; 0; 0; 0]
CCBELL₅ = CCBELL * [0; 0; 0; 0; 0; 1; 0; 0]
CCBELL₆ = CCBELL * [0; 0; 0; 0; 0; 0; 1; 0]
CCBELL₇ = CCBELL * [0; 0; 0; 0; 0; 0; 0; 1]

@test exec(Circuit(3, H(0), H(2), CCNOT(0, 2, 1))) ≈ CCBELL₀ * CCBELL₀'
@test exec(Circuit(3, X(0), H(0), H(2), CCNOT(0, 2, 1))) ≈ CCBELL₁ * CCBELL₁'
@test exec(Circuit(3, X(1), H(0), H(2), CCNOT(0, 2, 1))) ≈ CCBELL₂ * CCBELL₂'
@test exec(Circuit(3, X(0), X(1), H(0), H(2), CCNOT(0, 2, 1))) ≈ CCBELL₃ * CCBELL₃'
@test exec(Circuit(3, X(2), H(0), H(2), CCNOT(0, 2, 1))) ≈ CCBELL₄ * CCBELL₄'
@test exec(Circuit(3, X(0), X(2), H(0), H(2), CCNOT(0, 2, 1))) ≈ CCBELL₅ * CCBELL₅'
@test exec(Circuit(3, X(1), X(2), H(0), H(2), CCNOT(0, 2, 1))) ≈ CCBELL₆ * CCBELL₆'
@test exec(Circuit(3, X(0), X(1), X(2), H(0), H(2), Gate(matrix_X, 1, [0,2]))) ≈ CCBELL₇ * CCBELL₇'

# Testing X.
@test exec(Circuit(1, X(0))) == [0 0; 0 1]
@test exec(Circuit(1, H(0), X(0))) ≈ [.5 .5; .5 .5]

# Testing Y.
@test exec(Circuit(1, Y(0))) == [0 0; 0 1]
@test exec(Circuit(1, X(0), Y(0))) == [1 0; 0 0]

# Testing Z.
@test exec(Circuit(1, Z(0))) == [1 0; 0 0]
@test exec(Circuit(1, X(0), Z(0))) == [0 0; 0 1]
@test exec(Circuit(1, H(0), Z(0))) ≈ [.5 -.5; -.5 .5]
@test exec(Circuit(1, X(0), H(0), Z(0))) ≈ [.5 .5; .5 .5]

# Testing H.
@test exec(Circuit(1, H(0))) ≈ [.5 .5; .5 .5]

# Testing PHASE.
@test exec(Circuit(1, H(0), PHASE(pi / 2, 0))) ≈ [.5 complex(0, -.5); complex(0, .5) .5]

# Testing S.
@test exec(Circuit(1, H(0), S(0))) ≈ [.5 -.5im; .5im .5]
@test exec(Circuit(1, X(0), H(0), S(0))) ≈ [.5 .5im; -.5im .5]

# Testing T.
@test exec(Circuit(1, H(0), T(0))) ≈ [.5 (1 - 1im) / (2 * sqrt(2)); (1 + 1im) / (2 * sqrt(2)) .5]

# Testing CZ.
@test exec(Circuit(2, H(0), H(1), CZ(1, 0))) ≈ [.5; .5; .5; -.5] * [.5 .5 .5 -.5]

# Testing RX.
@test exec(Circuit(1, RX(pi, 0))) ≈ [0 0; 0 1]
@test exec(Circuit(1, RX(pi / 2, 0))) ≈ [.5 .5im; -.5im .5]

# Testing RY.
@test exec(Circuit(1, RY(pi, 0))) ≈ [0 0; 0 1]
@test exec(Circuit(1, RY(pi / 2, 0))) ≈ [.5 .5; .5 .5]

# Testing RZ.
@test exec(Circuit(1, RZ(pi, 0))) ≈ [1 0; 0 0]
@test exec(Circuit(1, H(0), RZ(pi, 0))) ≈ [.5 -.5; -.5 .5]
@test exec(Circuit(1, H(0), RZ(pi / 2, 0))) ≈ [.5 -.5im; .5im .5]

# Testing SWAP.
@test exec(Circuit(2, X(0), SWAP(0, 1))) == [0; 0; 1; 0] * [0 0 1 0]
@test exec(Circuit(4, X(1), X(0), X(3), CCSWAP(0, 3, 1, 2))) == [0; 1] ⊗ [0; 0; 1; 0] ⊗ [0; 1] * ([0; 1] ⊗ [0; 0; 1; 0] ⊗ [0; 1])'
# Control-SWAP does nothing because one of the control qubits is 0.
@test exec(Circuit(4, X(1), X(3), CCSWAP(0, 3, 1, 2))) == [0; 1] ⊗ [1; 0] ⊗ [0; 1] ⊗ [1; 0] * ([0; 1] ⊗ [1; 0] ⊗ [0; 1] ⊗ [1; 0])'
@test exec(Circuit(2, X(0), H(0), SWAP(0, 1))) == [1 / √2; 0; -1 / √2; 0] * [1 / √2 0 -1 / √2 0]
@test exec(Circuit(2, X(0), H(0), SWAP(0, 1))) == [1 / √2; 0; -1 / √2; 0] * [1 / √2 0 -1 / √2 0]
@test exec(Circuit(3, X(0), CSWAP(2, 0, 1))) == [0; 1; 0; 0; 0; 0; 0; 0] * [0 1 0 0 0 0 0 0]

# Testing parametric RX.
c = Circuit(1, RX("theta", 0))
@test exec(c, Dict("theta" => Float64(pi))) ≈ [0 0; 0 1]
reset_state!(c)
@test exec(c, Dict("theta" => Float64(pi / 2))) ≈ 1 / 2 * [1 1im; -1im 1]

# Testing parametric RY.
c = Circuit(1, RY("theta", 0))
@test exec(c, Dict("theta" => Float64(pi))) ≈ [0 0; 0 1]
reset_state!(c)
@test exec(c, Dict("theta" => Float64(pi / 2))) ≈ [.5 .5; .5 .5]

# Testing parametric RZ.
c = Circuit(1, H(0), RZ("theta", 0))
@test exec(c, Dict("theta" => Float64(pi))) ≈ [.5 -.5; -.5 .5]
reset_state!(c)
@test exec(c, Dict("theta" => Float64(pi / 2))) ≈ [.5 -.5im; .5im .5]

# Testing custom gate without parameters.
@test exec(Circuit(1, Gate(matrix_X, 0), X(0))) == [1 0; 0 0]
@test exec(Circuit(1, Gate(matrix_X, 0), H(0), X(0))) ≈ [.5 -.5; -.5 .5]

# Testing parametric custom gate.
matrix_x_func = function (params)
    angle = params["alpha"]
    return Matrix{Complex{Float64}}([cos(angle / 2) complex(0, -1) * sin(angle / 2); complex(0, -1) * sin(angle / 2) cos(angle / 2)])
end
c = Circuit(1, Gate(matrix_x_func, 0))
@test exec(c, Dict("alpha" => Float64(pi))) ≈ [0 0; 0 1]
reset_state!(c)
@test exec(c, Dict("alpha" => Float64(pi / 2))) ≈ 1 / 2 * [1 1im; -1im 1]

# Testing adding control bits to a gate.
@test exec(Circuit(2, CNOT(1, 0))) == exec(Circuit(2, controlled(1, X(0))))
@test exec(Circuit(3, CCNOT(1, 2, 0))) == exec(Circuit(3, controlled(1, controlled(2, X(0)))))
@test exec(Circuit(3, CCNOT(1, 2, 0))) == exec(Circuit(3, controlled(1, controlled(2, X(0)))))
@test exec(Circuit(3, CCNOT(1, 2, 0))) == exec(Circuit(3, controlled([1,2], X(0))))

# Testing += sintax for adding gates.
c = Circuit(2)
c += X(0)
c += X(1)
exec(c)
@test c.density_matrix == exec(Circuit(2, X(0), X(1)))

# Testing amplitude damping. Amplitude damping causes |1> to decay to |0> with some probability.
# When we start in |0> state, damping kraus with any |1> to |0> decay probability never gets applied,
# residual kraus always gets applied, so the state doesn't change.
@test exec(Circuit(1, noisify(ID(0), damping_kraus_map(.1)))) == [1 0; 0 0]
# Starting in |1> state, setting |1> to |0> decay probability to 0.
# Residual kraus should be applied and the state should remain |1>.
@test exec(Circuit(1, damp_amplitude(X(0), 0.0))) == [0 0; 0 1]
# Starting in |1> state, setting |1> to |0> decay probability to 1.
# Damping kraus should be applied and the state should decay to |0>.
@test exec(Circuit(1, damp_amplitude(X(0), 1.0))) == [1 0; 0 0]

# Starting in |+> state, setting |1> to |0> decay probability to 0.
# Residual kraus should be applied and the state should remain |1>.
@test exec(Circuit(1, damp_amplitude(H(0), 0.0))) ≈ [.5 .5; .5 .5]

# Starting in |+> state, setting |1> to |0> decay probability to 1.
# Damping kraus should be applied and the state should decay to |0>.
@test exec(Circuit(1, damp_amplitude(H(0), 1.0))) ≈ [1 0; 0 0]

# Testing private functions.

CX = [1 0 0 0;
      0 1 0 0;
      0 0 0 1;
      0 0 1 0]

CCX = [1 0 0 0 0 0 0 0;
       0 1 0 0 0 0 0 0;
       0 0 1 0 0 0 0 0;
       0 0 0 1 0 0 0 0;
       0 0 0 0 1 0 0 0;
       0 0 0 0 0 1 0 0;
       0 0 0 0 0 0 0 1;
       0 0 0 0 0 0 1 0]

CCCX = [1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1;
        0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0]

XC = [1 0 0 0;
      0 0 0 1;
      0 0 1 0;
      0 1 0 0]

XCC = [1 0 0 0 0 0 0 0;
       0 1 0 0 0 0 0 0;
       0 0 1 0 0 0 0 0;
       0 0 0 0 0 0 0 1;
       0 0 0 0 1 0 0 0;
       0 0 0 0 0 1 0 0;
       0 0 0 0 0 0 1 0;
       0 0 0 1 0 0 0 0]

CXC = [1 0 0 0 0 0 0 0;
       0 1 0 0 0 0 0 0;
       0 0 1 0 0 0 0 0;
       0 0 0 1 0 0 0 0;
       0 0 0 0 1 0 0 0;
       0 0 0 0 0 0 0 1;
       0 0 0 0 0 0 1 0;
       0 0 0 0 0 1 0 0]

CIXC = [1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1;
        0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0]

CXIC = [1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0;
        0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1;
        0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0;
        0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0;
        0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0]

τ = 1 / √2
CHIC = [1 0 0 0 0 0 0 0 0 0 0 0 0  0 0  0;
        0 1 0 0 0 0 0 0 0 0 0 0 0  0 0  0;
        0 0 1 0 0 0 0 0 0 0 0 0 0  0 0  0;
        0 0 0 1 0 0 0 0 0 0 0 0 0  0 0  0;
        0 0 0 0 1 0 0 0 0 0 0 0 0  0 0  0;
        0 0 0 0 0 1 0 0 0 0 0 0 0  0 0  0;
        0 0 0 0 0 0 1 0 0 0 0 0 0  0 0  0;
        0 0 0 0 0 0 0 1 0 0 0 0 0  0 0  0;
        0 0 0 0 0 0 0 0 1 0 0 0 0  0 0  0;
        0 0 0 0 0 0 0 0 0 τ 0 0 0  τ 0  0;
        0 0 0 0 0 0 0 0 0 0 1 0 0  0 0  0;
        0 0 0 0 0 0 0 0 0 0 0 τ 0  0 0  τ;
        0 0 0 0 0 0 0 0 0 0 0 0 1  0 0  0;
        0 0 0 0 0 0 0 0 0 τ 0 0 0 -τ 0  0;
        0 0 0 0 0 0 0 0 0 0 0 0 0  0 1  0;
        0 0 0 0 0 0 0 0 0 0 0 τ 0  0 0 -τ]

HH = Matrix{Complex{Float64}}([.5 .5 .5 .5; .5 -.5 .5 -.5; .5 .5 -.5 -.5; .5 -.5 -.5 .5])

CHH = [1 0 0 0  0   0   0   0;
       0 1 0 0  0   0   0   0;
       0 0 1 0  0   0   0   0;
       0 0 0 1  0   0   0   0;
       0 0 0 0 .5  .5  .5  .5;
       0 0 0 0 .5 -.5  .5 -.5;
       0 0 0 0 .5  .5 -.5 -.5;
       0 0 0 0 .5 -.5 -.5  .5]

@test sim.FullRegisterGate.expandGateToFullRegister(1, matrix_H, 0) == matrix_H
@test sim.FullRegisterGate.expandGateToFullRegister(2, HH, 0) == HH
@test sim.FullRegisterGate.expandGateToFullRegister(3, HH, 0, [2]) == CHH
@test sim.FullRegisterGate.expandGateToFullRegister(6, HH, 2, [4]) == matrix_I ⊗ CHH ⊗ matrix_I ⊗ matrix_I
@test sim.FullRegisterGate.expandGateToFullRegister(6, matrix_H, 3, [4,1]) == matrix_I ⊗ CHIC ⊗ matrix_I
@test sim.FullRegisterGate.expandGateToFullRegister(2, matrix_X, 0, [1]) == CX
@test sim.FullRegisterGate.expandGateToFullRegister(3, matrix_X, 0, [1,2]) == CCX
@test sim.FullRegisterGate.expandGateToFullRegister(4, matrix_X, 0, [2,1,3]) == CCCX
@test sim.FullRegisterGate.expandGateToFullRegister(3, matrix_X, 0, [1]) == matrix_I ⊗ CX
@test sim.FullRegisterGate.expandGateToFullRegister(3, matrix_X, 1, [2]) == CX ⊗ matrix_I
@test sim.FullRegisterGate.expandGateToFullRegister(4, matrix_X, 0, [1]) == matrix_I ⊗ matrix_I ⊗ CX
@test sim.FullRegisterGate.expandGateToFullRegister(5, matrix_X, 0, [1,2]) == matrix_I ⊗ matrix_I ⊗ CCX
@test sim.FullRegisterGate.expandGateToFullRegister(2, matrix_X, 1, [0]) == XC
@test sim.FullRegisterGate.expandGateToFullRegister(3, matrix_X, 1, [0]) == matrix_I ⊗ XC
@test sim.FullRegisterGate.expandGateToFullRegister(3, matrix_X, 2, [0,1]) == XCC
@test sim.FullRegisterGate.expandGateToFullRegister(4, matrix_X, 2, [1,0]) == matrix_I ⊗ XCC
@test sim.FullRegisterGate.expandGateToFullRegister(3, matrix_X, 1, [0,2]) == CXC
@test sim.FullRegisterGate.expandGateToFullRegister(4, matrix_X, 2, [1,3]) == CXC ⊗ matrix_I
@test sim.FullRegisterGate.expandGateToFullRegister(5, matrix_X, 2, [3,1]) == matrix_I ⊗ CXC ⊗ matrix_I
@test sim.FullRegisterGate.expandGateToFullRegister(4, matrix_X, 1, [0,3]) == CIXC
@test sim.FullRegisterGate.expandGateToFullRegister(3, matrix_X, 2, [1]) == XC ⊗ matrix_I
@test sim.FullRegisterGate.expandGateToFullRegister(4, matrix_X, 2, [3,0]) == CXIC
@test sim.FullRegisterGate.expandGateToFullRegister(7, matrix_X, 5, [3,6]) == CXIC ⊗ matrix_I ⊗ matrix_I ⊗ matrix_I
@test sim.FullRegisterGate.expandGateToFullRegister(6, matrix_X, 3, [1,4]) == matrix_I ⊗ CXIC ⊗ matrix_I
@test sim.FullRegisterGate.expandGateToFullRegister(1, matrix_S, 0) == matrix_S
# Testing asymmetric gates to ensure the gate doesn't get transposed by mistake.
@test sim.FullRegisterGate.expandGateToFullRegister(1, damping_residual_kraus(), 0) == damping_residual_kraus()
@test sim.FullRegisterGate.expandGateToFullRegister(3, damping_residual_kraus(), 1) == matrix_I ⊗ damping_residual_kraus() ⊗ matrix_I
@test sim.FullRegisterGate.expandGateToFullRegister(4, damping_residual_kraus(), 2) == matrix_I ⊗ damping_residual_kraus() ⊗ matrix_I ⊗ matrix_I
@test sim.FullRegisterGate.expandGateToFullRegister(1, damping_kraus(), 0) == damping_kraus()
@test sim.FullRegisterGate.expandGateToFullRegister(3, damping_kraus(), 1) == matrix_I ⊗ damping_kraus() ⊗ matrix_I
@test sim.FullRegisterGate.expandGateToFullRegister(4, damping_kraus(), 1) == matrix_I ⊗ matrix_I ⊗ damping_kraus() ⊗ matrix_I

@test sim.FullRegisterGate.gate_full_register_bitmask(1, 0) == 0b1
@test sim.FullRegisterGate.gate_full_register_bitmask(1, 1) == 0b10
@test sim.FullRegisterGate.gate_full_register_bitmask(1, 2) == 0b100
@test sim.FullRegisterGate.gate_full_register_bitmask(2, 0) == 0b11
@test sim.FullRegisterGate.gate_full_register_bitmask(2, 1) == 0b110
@test sim.FullRegisterGate.gate_full_register_bitmask(2, 2) == 0b1100
@test sim.FullRegisterGate.gate_full_register_bitmask(3, 0) == 0b111
@test sim.FullRegisterGate.gate_full_register_bitmask(3, 1) == 0b1110
@test sim.FullRegisterGate.gate_full_register_bitmask(3, 2) == 0b11100

@test sim.FullRegisterGate.control_full_register_bitmask(Array{Int64,1}([])) == 0
@test sim.FullRegisterGate.control_full_register_bitmask([0]) == 1
@test sim.FullRegisterGate.control_full_register_bitmask([0, 1]) == 0b11
@test sim.FullRegisterGate.control_full_register_bitmask([1, 5]) == 0b100010
@test sim.FullRegisterGate.control_full_register_bitmask([1, 7, 5]) == 0b10100010

@test sim.FullRegisterGate.all_control_bits_set(0b11, 0b10) == true
@test sim.FullRegisterGate.all_control_bits_set(0b1100, 0b1010) == false
@test sim.FullRegisterGate.all_control_bits_set(0b0000, 0b1010) == false
@test sim.FullRegisterGate.all_control_bits_set(0b1111, 0b0101) == true
@test sim.FullRegisterGate.all_control_bits_set(0b0, 0b0) == true
@test sim.FullRegisterGate.all_control_bits_set(0, 0) == true
@test sim.FullRegisterGate.all_control_bits_set(0xffffffffffffffff, 0xffffffffffffffff) == true
@test sim.FullRegisterGate.all_control_bits_set(0x8fffffffffffffff, 0xffffffffffffffff) == false

@test damping_kraus_map(.1) == [[1 0; 0 √.9], [0 √.1; 0 0]]

@test sim.value_index_by_probabilities(.0, [.3, .4, .2, .1]) == 1
@test sim.value_index_by_probabilities(.35, [.3, .4, .2, .1]) == 2
@test sim.value_index_by_probabilities(.8, [.3, .4, .2, .1]) == 3
@test sim.value_index_by_probabilities(.95, [.3, .4, .2, .1]) == 4
@test sim.value_index_by_probabilities(.9999, [.3, .4, .2, .1]) == 4
@test sim.value_index_by_probabilities(1.0, [.3, .4, .2, .1]) == 4


text:
using FractionalDiffEq
using Plots

α = 0.7
β = 1.8
κ = 1
T = 2
m = 21
n = 148

result = solve(α, β, κ, T, m, n, FPDEMatrixDiscrete())

XX, YY = meshgrid(0.05^2/6 .*(0:147), 0:0.05:1)

plotlyjs()

plot(XX, YY, result, st=:surface)

ext: jl


text:
export downsample

"""

  downsample(C::SeisChannel,fs::Real)

Downsample SeisChannel sampling rate to frequency `fs`.

For best results, lowpass filter data to `fs` before downsampling.
Implements the weighted average slopes interpolation scheme proposed in
[Wiggins1976] for evenly sampled data from
obspy.signal.interpolation.weighted_average_slopes.
"""
function downsample(C::SeisChannel,fs::Float64)
    dt = float(fs)

    if dt <= 0.
        throw(DomainError(dt, "Sampling rate must be nonnegative"))
    end

    dt = 1. / fs
    if fs == C.fs
        return C
    end

    old_start = C.t[1,2] * 1e-6
    old_dt = 1. / C.fs
    starttime, endtime = SeisIO.t_win(C.t,C.fs)
    starttime, endtime = starttime *1e-6, endtime * 1e-6
    npts = Int(floor((endtime-starttime)/dt)) + 1
    C.x = weighted_average_slopes(C.x,old_start,old_dt,starttime,dt,npts)
    C.fs = fs
    C.t[2,1] = npts
    return C
end

"""
    downsample(S::SeisData, fs::Float64)

Downsample SeisData sampling rate to frequency `fs`.
"""
function downsample(S::SeisData,fs::Float64)
    @inbounds for i = 1:S.n
        S[i]= downsample(S[i],fs)
    end
    return S
end

"""

  weighted_average_slopes(A::AbstractArray,old_start::Real,old_dt::Real,
                                 new_start::Real,new_dt::Real,new_npts::Int)

"""
function weighted_average_slopes(A::AbstractArray,
                                 old_start::Float64,
                                 old_dt::Float64,
                                 new_start::Float64,
                                 new_dt::Float64,
                                 new_npts::Int)
    old_end, new_end = validate_parameters(A, old_start, old_dt, new_start,
                                           new_dt, new_npts)
    new_time_array = Array(range(new_start, stop=new_end, length=new_npts))
    m = diff(A)
    m ./= old_dt
    w = abs.(m)
    w .= 1. ./ clamp.(w, eps(Float64), maximum(w))

    slope = zeros(eltype(A),length(A))
    slope[1] = m[1]
    for ii = 2:length(m)
        slope[ii] = (w[ii-1] * m[ii-1] + w[ii] * m[ii]) / (w[ii-1] + w[ii])
    end
    slope[end] = m[end]

    # If m_i and m_{i+1} have opposite signs then set the slope to zero.
    # This forces the curve to have extrema at the sample points and not
    # in-between.
    sign_change = find_nonzero(diff(sign.(m)))
    for ii = 1:length(sign_change)
        slope[sign_change[ii] + 1] = 0.
    end

    # Create interpolated value using hermite interpolation. In this case
    # it is directly applicable as the first derivatives are known.
    return hermite_interpolation(A,slope,new_time_array,old_dt,old_start)
end

"""

  hermite_interpolation(y_in, slope, x_out, h, x_start)

Hermite interpolation when zeroth and first derivatives are given for each time step.

# Arguements
- y_in: The data values to be interpolated.
- slope: The desired slope at each data point.
- x_out: The point at which to interpolate.
- h: The sample interval for y_in.
- x_start: Time of the first sample in y_in.
"""
function hermite_interpolation(y_in::AbstractArray,
                               slope::AbstractArray,
                               x_out::AbstractArray,
                               h::Float64,x_start::Float64)
    len_out = length(x_out)
    len_in = length(y_in)
    y_out = zeros(len_out)
    y_out[1] = y_in[1]

    for idx = 2:len_out
        i = (x_out[idx] - x_start) / h
        i_0 = convert(Int64,round(i))
        i_1= i_0 + 1

        if i == float(i_0)
            y_out[idx] = y_in[i_0]
            continue
        end

        t = i - float(i_0)
        a_0 = y_in[i_0]
        a_1 = y_in[i_1]
        b_minus_1 = h * slope[i_0]
        b_plus_1 = h * slope[i_1]
        b_0 = a_1 - a_0
        c_0 = b_0 - b_minus_1
        c_1 = b_plus_1 - b_0
        d_0 = c_1 - c_0

        y_out[idx] = a_0 + (b_0 + (c_0 + d_0 * t) * (t - 1.)) * t
    end
    return y_out
end

"""

  validate_parameters(data, old_start, old_dt, new_start, new_dt, new_npts)

Validates the parameters for various interpolation functions.

Returns the old and the new end.

"""
function validate_parameters(A::AbstractArray,
                             old_start::Float64,
                             old_dt::Float64,
                             new_start::Float64,
                             new_dt::Float64,
                             new_npts::Int)
    if new_dt <= 0.
        throw(DomainError(new_dt, "New sampling rate must be nonnegative"))
    end

    if ndims(A) != 1 || size(A) == (0,)
        throw(ArgumentError("Not a 1d array."))
    end

    old_end = old_start + old_dt * (length(A) - 1)
    new_end = new_start + new_dt * (new_npts - 1)

    if old_start > new_start || old_end < new_end
        throw(ArgumentError("The new array must be fully contained in the old array. No extrapolation can be performed."))
    end

    return old_end, new_end
end

function find_nonzero(c::Array{<:Union{Float32,Float64},1})
    a = similar(c, Int)
    count = 1
    @inbounds for i in eachindex(c)
        a[count] = i
        count += (c[i] != zero(eltype(c)))
    end
    return resize!(a, count-1)
end


text:
#function blkd(A :: Array,B :: Array)
#	# block diagonal matrix of A and B
#	al = size(A,1)
#	ac = size(A,2)
#	bl = size(B,1)
#	bc = size(B,2)
#	C = zeros(al+bl, ac + bc)
#	C[1:al, 1:ac] = A
#	C[al+1:end, ac+1:end] = B
#	C
#end
function eig(p :: Phs)
	# computes the eigenvalues and eigenvectors of a linear PHS
	#  (which haves a Q matrix!)
	#
	#	constrained (G != 0) cases can be addressed using
	#  generalized eigenvalues
	#
	#
	#  TODO: eigenvalues around a linearized condition:
	#       just need to compute the Hamiltonian Hessian at a given point!
	#       For e.g.: hess = ForwardDiff.hessian(p.Hamiltonian)
	#                 Q = hess(x)
	#               --> this can be useful for validating our nonlinear system!
	#

	if isdefined(p, :Q)
		if isdefined(p,:G)
			nconst = size(p.G,2)
			I = eye(size(p.J,1))
			z = zeros(nconst,nconst)
			E = blkdiag(I,z)
			A = [p.J*p.Q p.G;
			     p.G'*p.Q p.G_D]
			e = eigfact(A,E)
			ind = sortperm(imag(e.values))
			return e.values[ind], e.vectors[ind,ind]
		else
			a,v = eig(p.J*p.Q)
			ind = sortperm(imag(a))
			return a[ind], v[:,ind]
		end
	else
		fprintln("Undefined Q matrix (nonlinear system?)!")
		return
	end
end

function eigdamp(p :: Phs)
	# computes the frequency and damping of PHS with damping
	#  (which haves a Q matrix!)
	#
	#	constrained (G != 0) cases can be addressed using
	#  generalized eigenvalues
	#

	if isdefined(p, :Q)
		if isdefined(p,:G)
			nconst = size(p.G,2)
			I = eye(size(p.J,1))
			z = zeros(nconst,nconst)
			E = blkdiag(I,z)
			A = [(p.J-p.R)*p.Q p.G;
			     p.G'*p.Q p.G_D]
			e = eigfact(A,E)
			ind = sortperm(imag(e.values))
			return e.values[ind], e.vectors[ind,ind]
		else
			a,v = eig((p.J-p.R)*p.Q)
			ind = sortperm(imag(a))
			return a[ind], v[:,ind]
		end
	else
		fprintln("Undefined Q matrix (nonlinear system?)!")
		return
	end
end

function damp(p :: Phs)
	a,v = eigdamp(p)
	return [abs(a) -cos(angle(a))]
end
function frequencies(a :: Array)
		imag(a[imag(a).>=0])
end

function frequencies(ph :: Phs)
	a,v = eig(ph)
	return frequencies(a)
end

max_stars_repo_name: flavioluiz/PortHamiltonian

max_issues_repo_name: flavioluiz/PortHamiltonian

max_forks_repo_name: flavioluiz/PortHamiltonian


text:
using StrideArraysCore, ThreadingUtilities, Aqua
using Test

function closeopensum(x)
  s = zero(eltype(x))
  @inbounds @simd for i ∈ StrideArraysCore.CloseOpen(length(x))
    s += x[i+1]
  end
  s
end
function closeopensumfastmath(x)
  s = zero(eltype(x))
  @inbounds @fastmath for i ∈ StrideArraysCore.CloseOpen(length(x))
    s += x[i+1]
  end
  s
end
function cartesianindexsum(A)
  s = zero(eltype(A))
  @inbounds @simd for I ∈ CartesianIndices(A)
    s += A[I]
  end
  s
end

@testset "StrideArraysCore.jl" begin

  Aqua.test_all(StrideArraysCore)

  @testset "StrideArrays Basic" begin
    Acomplex = StrideArray{Complex{Float64}}(undef, (StaticInt(4), StaticInt(5)))
    Acomplex .= rand.(Complex{Float64})
    @test StrideArray(Acomplex) ===  reinterpret(reshape, Complex{Float64}, reinterpret(reshape, Float64, Acomplex))
    @test StrideArray(Acomplex) ===  reinterpret(reshape, Complex{Float64}, reinterpret(reshape, Complex{UInt64}, Acomplex))
    @test StrideArraysCore.size(Acomplex) === StrideArraysCore.size(StrideArray(Acomplex)) === (StaticInt(4),StaticInt(5))
    @test StrideArraysCore.size(reinterpret(reshape, Float64, Acomplex)) === (StaticInt(2),StaticInt(4),StaticInt(5))
    A = rand(100, 100);
    B = copy(A);
    C = StrideArraysCore.PtrArray(A);
    @test similar(C) isa StrideArraysCore.StrideArray
    let D = similar(C,Float32)
      @test D isa StrideArraysCore.StrideArray
      @test eltype(D) === Float32
    end
    GC.@preserve A begin
      @test closeopensum(C) == closeopensum(A)
      @test closeopensumfastmath(C) == closeopensumfastmath(A)
      @test sum(C) == sum(A)
      @test closeopensum(C) ≈ closeopensumfastmath(C) ≈ sum(C)
      @test A == B
      C .*= 3;
      @test A == 3 .* B
      @test C == A
      @test C == 3 .* B

      D = copy(A);
      Cslice = view(C, 23:48, 17:89)
      Cslice .= 2;
      @test D != C
      D[23:48,17:89] .= 2;
      @test D == C
      @test C === view(C, :, :)
      @test @inferred(size(view(C, StaticInt(1):StaticInt(8), :), 1)) === 8
      @test @inferred(size(view(C, StaticInt(1):StaticInt(8), :), StaticInt(1))) === StaticInt(8)
      @test @inferred(StrideArraysCore.size(view(C, StaticInt(1):StaticInt(8), :), 1)) === 8
      @test @inferred(StrideArraysCore.size(view(C, StaticInt(1):StaticInt(8), :), StaticInt(1))) === StaticInt(8)

      @test C  isa PtrArray
      @test C' isa PtrArray
      @test @inferred(permutedims(C, Val((2,1)))) isa PtrArray
      @test @inferred(adjoint(C)) == D'
      @test @inferred(axes(C)) == axes(D)
      @test @inferred(axes(C')) == axes(D')
      @test @inferred(eachindex(view(C, :, 2:6))) == 1:(5*size(C,1))
      @test @inferred(eachindex(view(C', 2:6, :)')) == 1:(5*size(C,1))
      @test @inferred(eachindex(view(C, 2:6, :))) == CartesianIndices((5, size(C,2)))
      @test @inferred(axes(StrideArraysCore.zview(C, StaticInt(2):StaticInt(6), :), StaticInt(1))) === StrideArraysCore.CloseOpen(StaticInt(5))
      @test @inferred(axes(StrideArraysCore.zview(C, StaticInt(2):StaticInt(6), :), StaticInt(1))) === StrideArraysCore.CloseOpen(StaticInt(5))
      @test @inferred(length(axes(StrideArraysCore.zview(C, StaticInt(2):StaticInt(6), :), StaticInt(1)))) === StaticInt(5)
      @test @inferred(length(axes(StrideArraysCore.zview(C, StaticInt(2):StaticInt(6), :), 1))) === 5
      ax1, ax2 = axes(StrideArraysCore.zview(C, StaticInt(2):StaticInt(6), :))
      @test StrideArraysCore.ArrayInterface.known_length(ax1) == 5
      @test StrideArraysCore.ArrayInterface.known_length(ax2) === nothing
      @test StrideArraysCore.ArrayInterface.known_first(ax1) == 0
      @test StrideArraysCore.ArrayInterface.known_first(ax2) == 0
      @test StrideArraysCore.ArrayInterface.known_step(ax1) == 1
      @test StrideArraysCore.ArrayInterface.known_step(ax2) == 1
      @test StrideArraysCore.ArrayInterface.known_last(ax1) == 4
      @test StrideArraysCore.ArrayInterface.known_last(ax2) === nothing
    end
    W = rand(2,3,4);
    X = PtrArray(W);
    GC.@preserve W begin
      @test W == X
      @test permutedims(W, (1,2,3)) == @inferred(permutedims(X, Val((1,2,3))))
      @test permutedims(W, (1,3,2)) == @inferred(permutedims(X, Val((1,3,2))))
      @test permutedims(W, (2,1,3)) == @inferred(permutedims(X, Val((2,1,3))))
      @test permutedims(W, (2,3,1)) == @inferred(permutedims(X, Val((2,3,1))))
      @test permutedims(W, (3,1,2)) == @inferred(permutedims(X, Val((3,1,2))))
      @test permutedims(W, (3,2,1)) == @inferred(permutedims(X, Val((3,2,1))))
      @test_throws BoundsError X[length(X) + 1]
      @test_throws BoundsError X[-4]
      @test_throws BoundsError X[2,5,3]
      @test cartesianindexsum(W) ≈ cartesianindexsum(X)
      @test iszero(@allocated cartesianindexsum(X))
    end
    @test X === PtrArray(pointer(X), size(X))
    y = rand(77);
    GC.@preserve y begin
      z = PtrArray(y);
      @test y == z
      @test y' == z'
      @test PtrArray(y') === z'
      @test PtrArray(transpose(y)) === transpose(PtrArray(y))
      @test pointer(y) === pointer(z)
      @test_throws BoundsError z[-8]
      @test_throws BoundsError z[88]
    end
  end
  @testset "ThreadingUtilities" begin
    xu = zeros(UInt, 100);
    x = rand(100); y = rand(100); z = rand(100);
    t = (x,y,z)
    pt, gt = StrideArraysCore.object_and_preserve(t)
    greet = ["Hello", "world", "hang", "in", "there"]
    pg, gg = StrideArraysCore.object_and_preserve(greet)
    ph, gh = StrideArraysCore.object_and_preserve(greet[1])
    pht, ght = StrideArraysCore.object_and_preserve((3,greet[1]))
    GC.@preserve xu gt gg gh ght begin
      ThreadingUtilities.store!(pointer(xu), pt, 0)
      @test @inferred(ThreadingUtilities.load(pointer(xu), typeof(pt), 0)) == (6sizeof(UInt), t)
      offset = sizeof(UInt)
      for a ∈ t
        _p, g = StrideArraysCore.object_and_preserve(a)
        ThreadingUtilities.store!(pointer(xu), _p, offset)
        offset, p = @inferred(ThreadingUtilities.load(pointer(xu), typeof(_p), offset))
        @test g === a
        @test p isa PtrArray
        @test pointer(p) === pointer(a)
        @test p !== a
        @test p == a
      end
      ThreadingUtilities.store!(pointer(xu), pg, offset)
      @test @inferred(ThreadingUtilities.load(pointer(xu), typeof(pg), offset)) === (offset+sizeof(UInt), greet)
      ThreadingUtilities.store!(pointer(xu), ph, offset)
      @test @inferred(ThreadingUtilities.load(pointer(xu), typeof(ph), offset)) === (offset+sizeof(UInt), greet[1])
      ThreadingUtilities.store!(pointer(xu), pht, offset)
      @test @inferred(ThreadingUtilities.load(pointer(xu), typeof(pht), offset)) === (offset+2sizeof(UInt), (3,greet[1]))
    end

    greetsa = StrideArray(greet);
    @test greetsa == greet
    greetsa[1] = "howdy"
    greetsa[5] = "there!"
    @test greet[1] == "howdy"
    @test greet[5] == "there!"
    A = rand(ComplexF64, 4, 5);
    Asa = StrideArray(A);
    @test A == Asa
    Asa[3,4] = 1234.5 + 678.9im;
    @test A[3,4] == 1234.5 + 678.9im;
  end
  @testset "StrideArrays Initialization" begin
    A = StrideArray{Float64}(undef, (3, 5));
    @test StrideArraysCore.size(A) === (3, 5)
    @test StrideArraysCore.strides(A) === (StaticInt(1),3)
    @test StrideArraysCore.static_length(A) === 15

    B = StrideArray(undef, (StaticInt(3), 5));
    @test StrideArraysCore.strides(B) === (StaticInt(1), StaticInt(3))
    @test StrideArraysCore.size(B) === (StaticInt(3), 5);
    @test StrideArraysCore.static_length(B) === 15

    D = StrideArray(undef, (StaticInt(3), StaticInt(5)));
    @test StrideArraysCore.strides(D) === (StaticInt(1), StaticInt(3))
    @test StrideArraysCore.size(D) === (StaticInt(3), StaticInt(5));
    @test StrideArraysCore.static_length(D) === StaticInt(15)
    @test D isa StrideArraysCore.StaticStrideArray
    for C ∈ [A,B,D]
      @test strides(C) === (1, 3)
      @test size(C) === (3, 5);
      @test StrideArraysCore.offsets(C) === (StaticInt(1),StaticInt(1))
      @test StrideArraysCore.offsets(StrideArraysCore.zeroindex(C)) === (StaticInt(0),StaticInt(0))
      C[2,3] = 4
      StrideArraysCore.zeroindex(C)[2,3] = -10
      @test C[2,3] === StrideArraysCore.zeroindex(C)[1,2] === 4.0
      @test C[3,4] === StrideArraysCore.zeroindex(C)[2,3] === -10.0
    end
    @test all(iszero, StrideArray(zero, static(4), static(8)))
    @test all(iszero, StrideArray(zero, 1000, 2000))
    @test all(isone, StrideArray(one, static(4), static(8)))
    @test all(isone, StrideArray(one, 100, 200))
  end
  @testset "views" begin
    A = StrideArray{Float64}(undef, (100,100)) .= rand.()
    vA = view(A, 3:40, 2:50)
    vAslice = view(A, :, 2:50)
    vaslice = view(A, 2:50, 4)
    x = StrideArray{Float64}(undef, (100,)) .= rand.()
    vxslice = view(x, 2:50)
    for (i,n) in enumerate(2:50)
      @test vaslice[i] == A[n,4]
      @test vxslice[i] == x[n]
      for (j,m) in enumerate(3:40)
        @test A[m, n] == vA[j,i]
      end
      for j in 1:50
        @test A[j,n] == vAslice[j,i]
      end
    end
  end
  # @testset "reinterpret" begin

  # end
end


text:
using Test, HoeffdingD
ab =[
    -1.0	2.03
    -1.0	2.11
    -1.0	2.03
    -0.8	1.58
    -0.8	1.79
    -0.8	1.58
    -0.6	1.35
    -0.6	1.46
    -0.6	1.37
    -0.4	1.28
    -0.4	1.15
    -0.4	1.19
    -0.2	0.93
    -0.2	0.99
    -0.2	0.99
    0.0 1.03
    0.0 0.98
    0.0 1.02
    0.2 1.00
    0.2 1.06
    0.2 1.00
    0.4 1.31
    0.4 1.33
    0.4 0.99
    0.6 1.28
    0.6 1.42
    0.6 1.50
    0.8 1.68
    0.8 1.64
    0.8 1.68
    1.0 2.07
    1.0 1.94
    1.0 2.14 ]

a = ab[:, 1]
b = ab[:, 2]

@testset "HoeffdingD" begin
    #results values are from SAS PROC Corr
    @test isapprox(0.1108736664, HoeffdingD.hoeffdingd(a, b))
    @test isapprox(0.1108736664, HoeffdingD.hoeffdingd(b, a))
    @test isapprox(0.9665923259, HoeffdingD.hoeffdingd(b, b))
    @test isapprox(0.8450535317, HoeffdingD.hoeffdingd(a, a))

    @test isapprox(0.6875, HoeffdingD.hoeffdingd(
            [9, 8,  7,  missing,    missing,    2,          5, 10, 2],
            [6, 5,  4,  missing,    3,          missing,    -5, 1000, 4]))

    @test isapprox([0.845053531701891 0.110873666447568; 0.110873666447568 0.9665923258586983], HoeffdingD.hoeffdingd(ab))

    linear_f(x) = @. 3x + 2
    quad_f(x) = @. x^2 + 1
    x_test = -1:0.05:1
    y_linear = linear_f.(x_test)
    y_quad = quad_f.(x_test)
    xyy = hcat(x_test, y_linear, y_quad)
    @test isapprox( [1.0 1.0 0.2183712793468891; 1.0 1.0 0.2183712793468891; 0.2183712793468891 0.2183712793468891 0.9369001852153328] , HoeffdingD.hoeffdingd(xyy))

end

@testset "HoeffdingD α" begin
    @test_throws DomainError HoeffdingD.hoeffdingd(ab, -1)
    @test [true false; false true] == last.(HoeffdingD.hoeffdingd(ab, .001))
    @test [true true; true true] == last.(HoeffdingD.hoeffdingd(ab, .05))
end

using Distributions

@testset "HoeffdingD independence" begin
    x = collect(1:200)
    y_uniform = rand(Distributions.Uniform(2,5), length(x))
    @test isapprox(0, last(HoeffdingD.hoeffdingd(x, y_uniform, 0.95)))
    @test isapprox(0, last(HoeffdingD.hoeffdingd(x, y_uniform, 0.5)))
    @test isapprox(0, last(HoeffdingD.hoeffdingd(x, y_uniform, 0.1)))
    @test isapprox(0, last(HoeffdingD.hoeffdingd(x, y_uniform, 0.05)))
    @test isapprox(0, last(HoeffdingD.hoeffdingd(x, y_uniform, 0.005)))
    @test isapprox(0, last(HoeffdingD.hoeffdingd(x, y_uniform, 0.001)))
end

# using BenchmarkTools

# @btime HoeffdingD.hoeffdingd(a, b)
# @btime HoeffdingD.hoeffdingd(b, a)
# @btime HoeffdingD.hoeffdingd(ab)

# @btime HoeffdingD.hoeffdingd(a, b, 0.05)
# @btime HoeffdingD.hoeffdingd(b, a, 0.05)
# @btime HoeffdingD.hoeffdingd(ab, 0.05)

# @btime HoeffdingD.hoeffdingd( [9, 8,  7.,  missing,    missing,    2,          5, 10, 2],
#             [6, 5,  4.,  missing,    3,          missing,    -5, 1000, 4])

max_stars_repo_name: ericqu/HoeffdingD

max_issues_repo_name: ericqu/HoeffdingD

max_forks_repo_name: ericqu/HoeffdingD


text:
# code for Video 5.2:
julia> using RDatasets
julia> iris = dataset("datasets", "iris")
julia> features = convert(Array, iris[:, 1:4]) # measurements of length / width
# 1) SepalLength, 2) SepalWidth, 3) PetalLength, 4) PetalWidth
julia> labels = convert(Array, iris[:, 5])     # species

julia> Pkg.add("DecisionTree")
julia> using DecisionTree

### Pruned Tree classifier:
# train full-tree classifier
julia> model = build_tree(labels, features)
Decision Tree
Leaves: 9
Depth:  5

# prune tree: merge leaves having >= 90% combined purity (default: 100%)
julia> model = prune_tree(model, 0.9)
Decision Tree
Leaves: 8
Depth:  5

# pretty print of the tree, to a depth of 5 nodes
julia> print_tree(model, 5)
Feature 3, Threshold 3.0 # PetalLength
L-> setosa : 50/50                                      (1)
R-> Feature 4, Threshold 1.8 # PetalWidth
    L-> Feature 3, Threshold 5.0 # PetalLength
        L-> versicolor : 47/48                          (2)
        R-> Feature 4, Threshold 1.6 # PetalWidth
            L-> virginica : 3/3                         (3)
            R-> Feature 1, Threshold 7 # SepalLength
                L-> versicolor : 2/2                    (4)
                R-> virginica : 1/1                     (5)
    R-> Feature 3, Threshold 4.9 # PetalLength
        L-> Feature 1, Threshold 6.0 # SepalLength
            L-> versicolor : 1/1                        (6)
            R-> virginica : 2/2                         (7)
        R-> virginica : 43/43                           (8)

julia> print_tree(model, 3)
Feature 3, Threshold 3.0 # PetalLength
L-> setosa : 50/50                                      (1)
R-> Feature 4, Threshold 1.8 # PetalWidth
    L-> Feature 3, Threshold 5.0 # PetalLength
        L-> versicolor : 47/48                          (2)
        R->
    R-> Feature 3, Threshold 4.9 # PetalLength
        L->
        R-> virginica : 43/43                           (3)

# apply learned model
## given feature measurements, predict the label (species)
julia> apply_tree(model, [5.9,3.0,5.1,1.9])
"virginica"
julia> apply_tree(model, [1.9,2.0,2.1,1.9])
"setosa"

# run n-fold cross validation for pruned tree,
# using 90% purity threshold purning, and 3 CV folds
accuracy = nfoldCV_tree(labels, features, 0.9, 3)

Fold 1
Classes:  Any["setosa","versicolor","virginica"]
Matrix:   3x3 Array{Int64,2}:
 19   0   0
  0  16   2
  0   0  13

Accuracy: 0.96
Kappa:    0.9396135265700483

Fold 2
Classes:  Any["setosa","versicolor","virginica"]
Matrix:   3x3 Array{Int64,2}:
 16   0   0
  1  14   2
  0   2  15

Accuracy: 0.9
Kappa:    0.8500299940011996

Fold 3
Classes:  Any["setosa","versicolor","virginica"]
Matrix:   3x3 Array{Int64,2}:
 15   0   0
  0  14   1
  0   1  19

Accuracy: 0.96
Kappa:    0.9393939393939393

Mean Accuracy: 0.94
3-element Array{Float64,1}:
 0.96
 0.9
 0.96

julia> using Base.Test
julia> @test mean(accuracy) > 0.8

### Adaptive-Boosted Decision Stumps Classifier
# train adaptive-boosted stumps, using 10 iterations
julia> model, coeffs = build_adaboost_stumps(labels, features, 10)
(Ensemble of Decision Trees
Trees:      10
Avg Leaves: 2.0
Avg Depth:  1.0,[0.3465735902799731,0.3711532082590139,0.38544475784944743,0.4093
352070410545,0.39211215411428185,0.43271872450399823,0.41523729251553393,0.413565
85168239374,0.42054477554575537,0.4125278240602734])

julia> apply_adaboost_stumps(model, coeffs, [5.9,3.0,5.1,1.9])
"versicolor"

julia> accuracy = nfoldCV_stumps(labels, features, 10, 3)

Fold 1
Classes:  Any["setosa","versicolor","virginica"]
Matrix:   3x3 Array{Int64,2}:
 19   0   0
  1  13   0
  0   1  16

Accuracy: 0.96
Kappa:    0.9394673123486682

Fold 2
Classes:  Any["setosa","versicolor","virginica"]
Matrix:   3x3 Array{Int64,2}:
 15  0   0
  0  3  18
  0  0  14

Accuracy: 0.64
Kappa:    0.489795918367347

Fold 3
Classes:  Any["setosa","versicolor","virginica"]
Matrix:   3x3 Array{Int64,2}:
 16   0   0
  0  14   1
  0   1  18

Accuracy: 0.96
Kappa:    0.9396863691194209

Mean Accuracy: 0.8533333333333334
3-element Array{Float64,1}:
 0.96
 0.64
 0.96

 ### Random Forest classifier
# train random forest classifier
# using 2 random features, 10 trees, and 0.5 portion of samples per tree (optional)
julia> model = build_forest(labels, features, 2, 10, 0.5)
Ensemble of Decision Trees
Trees:      10
Avg Leaves: 6.6
Avg Depth:  4.2

# apply learned model
julia> apply_forest(model, [5.9,3.0,5.1,1.9])
"virginica"

# run n-fold cross validation for forests
# using 2 random features, 10 trees, 3 folds and 0.5 of samples per tree (optional)
julia> accuracy = nfoldCV_forest(labels, features, 2, 10, 3, 0.5)
Fold 1
Classes:  Any["setosa","versicolor","virginica"]
Matrix:   3x3 Array{Int64,2}:
 17   0   0
  0  15   1
  0   1  16

Accuracy: 0.96
Kappa:    0.9399759903961584

Fold 2
Classes:  Any["setosa","versicolor","virginica"]
Matrix:   3x3 Array{Int64,2}:
 18   0   0
  0  14   2
  0   0  16

Accuracy: 0.96
Kappa:    0.9399038461538461

Fold 3
Classes:  Any["setosa","versicolor","virginica"]
Matrix:   3x3 Array{Int64,2}:
 15   0   0
  0  18   0
  0   1  16

Accuracy: 0.98
Kappa:    0.9698976520168573

Mean Accuracy: 0.9666666666666667
3-element Array{Float64,1}:
 0.96
 0.96
 0.98

max_stars_repo_path: Section 5/5_2.jl

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_path: Section 5/5_2.jl

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_path: Section 5/5_2.jl

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video


text:
using FactCheck

inc(x) = x + 1

facts("Succeeding examples") do

    @fact 1 --> 1

    context("You can define contexts") do
        @fact 1 --> 1
    end

    @fact_throws error("neat")

    context("group facts") do
        @fact 1 --> not(2)
        @fact 2 --> not(isodd)
    end

    context("generating facts") do
        for i=1:5
            @fact i --> i
        end
    end

end

facts("Failing examples") do

    @fact 1 --> 2

    context("strings are strings") do
        @fact "bar" --> "barr"
        @fact "baz" --> "bazz"
    end

    context("some numbers are even") do
        @fact 3 --> iseven
    end

    context("incrementing") do
        x = 10
        @fact inc(inc(inc(0))) --> 2
    end

    context("nonsense") do
        x = 5
        y = 10
        @fact x --> y
    end

    context("throws an error") do
        @fact error("foo") --> 1
    end

end

facts("Assertion helpers") do

    @fact 1 --> not(iseven)
    @fact iseven --> exactly(iseven)
    @fact 2.499999 --> roughly(2.5)
    @fact 1 --> less_than(2)
    @fact 1 --> less_than_or_equal(1)
    @fact 2 --> greater_than(1)
    @fact 2 --> greater_than_or_equal(2)

end


text:
# test_le.sac was generated in SAC 101.6a with "fg seismogram; write test_le.sac"
sac_file    = path*"/SampleFiles/SAC/test_le.sac"
sac_be_file = path*"/SampleFiles/SAC/test_be.sac"
sac_pz_file = path*"/SampleFiles/SAC/test_sac.pz"
sac_pz_wc   = path*"/SampleFiles/SAC/test_sac.*"
uw_file     = path*"/SampleFiles/UW/00012502123W"
sac_pat     = path*"/SampleFiles/SAC/*.sac"
sac_pz_out1 = path*"/local_sac_1.pz"
sac_pz_out2 = path*"/local_sac_2.pz"
sac_pz_out3 = path*"/local_sac_3.pz"
sac_pz_out4 = path*"/local_sac_4.pz"
sac_pz_out5 = path*"/local_sac_5.pz"
f_stub      = "1981.088.10.38.23.460"
f_out       = f_stub * "..CDV...R.SAC"
f_out_new   = f_stub * ".VU.CDV..NUL.R.SAC"

printstyled("  SAC\n", color=:light_green)
printstyled("    read\n", color=:light_green)
@test_throws ErrorException verified_read_data("sac", uw_file)

SAC1 = verified_read_data("sac", sac_file)[1]
@test ≈(SAC1.fs, 100.0)
@test ≈(length(SAC1.x), 1000)

# SAC with mmap
printstyled("      mmap\n", color=:light_green)
SACm = read_data("sac", sac_file, memmap=true)[1]
@test SAC1 == SACm

SAC2 = verified_read_data("sac", sac_file, full=true)[1]
@test ≈(1/SAC1.fs, SAC2.misc["delta"])
@test ≈(length(SAC1.x), SAC2.misc["npts"])

printstyled("      wildcard\n", color=:light_green)
SAC = verified_read_data("sac", sac_pat, full=true)

printstyled("      bigendian\n", color=:light_green)
SAC3 = verified_read_data("sac", sac_be_file, full=true)[1]
@test ≈(1/SAC3.fs, SAC3.misc["delta"])
@test ≈(length(SAC3.x), SAC3.misc["npts"])

redirect_stdout(out) do
  sachdr(sac_be_file)
end

printstyled("    write\n", color=:light_green)
S = SeisData(SAC2)
writesac(S) # change 2019-07-15 to cover writesac on GphysData
@test safe_isfile(f_out)
@test any([occursin(string("wrote to file ", f_out), S.notes[1][i]) for i in 1:length(S.notes[1])])

printstyled("      reproducibility\n", color=:light_green)
SAC4 = verified_read_data("sac", f_out, full=true)[1]
for f in SeisIO.datafields
  (f in [:src, :notes, :misc]) && continue
  @test isequal(getfield(SAC2, f), getfield(SAC4, f))
end
fn = f_out[3:end]
writesac(SAC4, fname=fn)
@test safe_isfile(fn)
SAC5 = verified_read_data("sac", fn, full=true)[1]
for f in SeisIO.datafields
  (f in [:src, :notes, :misc]) && continue
  @test isequal(getfield(SAC2, f), getfield(SAC4, f))
  @test isequal(getfield(SAC4, f), getfield(SAC5, f))
end
safe_rm(fn)

printstyled("      logging\n", color=:light_green)
@test any([occursin("write", n) for n in S.notes[1]])
redirect_stdout(out) do
  show_writes(S, 1)
  show_writes(S[1])
  show_writes(S)
end

SAC1.id = "VU.CDV..NUL"
SAC1.name = "VU.CDV..NUL"
writesac(SAC1)
@test safe_isfile(f_out_new)
writesac(SAC1, fname="POTATO.SAC")
@test safe_isfile("POTATO.SAC")

# testing custom naming formats
writesac(SAC1, fname="test_write_1.sac")
@test safe_isfile("test_write_1.sac")
safe_rm("test_write_1.sac")

# testing that appending ".sac" to a file string works
writesac(SAC1, fname="test_write_2", v=1)
@test safe_isfile("test_write_2.sac")
safe_rm("test_write_2.sac")

redirect_stdout(out) do
  writesac(SAC1, v=1)
end
@test safe_isfile(f_out_new)
safe_rm(f_out_new)

printstyled("      skip if fs==0.0\n", color=:light_green)
SAC1.id = "VU.FS0..NUL"
SAC1.fs = 0.0
writesac(SAC1)
@test safe_isfile(f_stub*"VU.FS0..NUL.R.SAC") == false

# SACPZ
printstyled("    SACPZ\n", color=:light_green)
printstyled("      read\n", color=:light_green)
S = read_meta("sacpz", sac_pz_wc)
S = read_meta("sacpz", sac_pz_file)
writesacpz(sac_pz_out1, S)
T = read_meta("sacpz", sac_pz_out1)
for f in (:n, :id, :name, :loc, :fs, :gain, :units)
  @test isequal(getfield(S, f), getfield(T, f))
end
for f in fieldnames(typeof(S.resp[1]))
  for i = 1:S.n
    @test isapprox(getfield(S.resp[i], f), getfield(T.resp[i], f))
  end
end
U = deepcopy(S)
for i = 1:U.n
  U.resp[i] = MultiStageResp(3)
  U.resp[i].stage[1] = S.resp[i]
end

printstyled("      write\n", color=:light_green)
writesacpz(sac_pz_out2, U)
T = read_meta("sacpz", sac_pz_out2)
for f in (:n, :id, :name, :loc, :fs, :gain, :units)
  @test isequal(getfield(S, f), getfield(T, f))
end
for f in fieldnames(typeof(S.resp[1]))
  for i = 1:S.n
    @test isapprox(getfield(S.resp[i], f), getfield(T.resp[i], f))
  end
end
U[1] = SeisChannel(id = "UW.HOOD..ENE")
writesacpz(sac_pz_out3, U)
read_meta!(S, "sacpz", sac_pz_out3)
for f in (:n, :id, :name, :loc, :fs, :gain, :units)
  @test isequal(getfield(S, f), getfield(T, f))
end
S = breaking_seis()[1:3]
S.resp[1].resp = rand(ComplexF64, 12, 2)
S.resp[3].stage[2] = nothing
writesacpz(sac_pz_out4, S)

printstyled("        GphysChannel\n", color=:light_green)
writesacpz(sac_pz_out5, S[1])

safe_rm(sac_pz_out1)
safe_rm(sac_pz_out2)
safe_rm(sac_pz_out3)
safe_rm(sac_pz_out4)
safe_rm(sac_pz_out5)


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export determinize_fsa
"""
  `determinize_fsa(fsa)`

Returns a deterministic finite state acceptor. `fsa` must be an acceptor ([`is_acceptor`](@ref)).

Requires the WFST to be defined in a left distributive and weakly divisible semiring.

See [M. Mohri, F. Pereira, Fernando, M. Riley, Michael "Speech Recognition with Weighted Finite-State Transducers" in Springer Handb. Speech Process. 2008](https://cs.nyu.edu/~mohri/pub/hbka.pdf) for details.
"""
function determinize_fsa(fst::WFST{W,A,D,I,O}) where {W,A,D,I,O}
  if !isleft(W)
    throw(ErrorException("Semiring $W is not left distributive, fst cannot be determinized."))
  end
  if !isweaklydivisible(W)
    throw(ErrorException("Semiring $W is not weakly divisible, fst cannot be determinized."))
  end
  dfst = WFST(get_isym(fst), get_osym(fst); W = W)
  i = get_initial(fst)
  w = get_initial_weight(fst,i)
  initial!(dfst,get_initial(fst))

  S = Queue{Set{Tuple{D,W}}}()
  qq = Set([(i,w)])
  cnt = 1
  qq2state = Dict( qq => cnt ) # new states
  enqueue!(S,qq)
  while !isempty(S)
    pp = dequeue!(S)
    lab2weight = Dict{D,W}() # set in line 8 of Fig.10 in [1]
    lab_nextstate2weight = Dict{D,Dict{D,W}}() # set in line 9 of Fig.10 in [1]

    for (p,v) in pp
      for arc in fst[p]
        lab, w, n = get_ilabel(arc), get_weight(arc), get_nextstate(arc)
        vw = v*w
        if get_ilabel(arc) in keys(lab2weight)
          lab2weight[lab] +=  vw
        else
          lab2weight[lab]  =  vw
        end
        if lab in keys(lab_nextstate2weight)
          if n in keys(lab_nextstate2weight[lab])
            lab_nextstate2weight[lab][n] += vw
          else
            lab_nextstate2weight[lab][n] = vw
          end
        else
          lab_nextstate2weight[lab] = Dict(n => vw)
        end
      end
    end

    for lab in keys(lab2weight)
      qq = Set{Tuple{D,W}}()
      wp = lab2weight[lab]
      for n in keys(lab_nextstate2weight[lab])
        vw = lab_nextstate2weight[lab][n]
        push!(qq,(n, one(W) / wp * vw))
      end
      if !(qq in keys(qq2state))
        cnt += 1
        qq2state[qq] = cnt
        src, dest = qq2state[pp], cnt
        if max(src,dest) - length(dfst) > 0
          add_states!(dfst, max(src,dest) - length(dfst) )  # expand states
        end
        push!(dfst[src],Arc{W,D}(lab,lab,wp,dest))
        rho = zero(W)
        any_final = false
        for (n,v) in qq
          if isfinal(fst,n)
            any_final = true
            rho += v*get_final_weight(fst,n)
          end
        end
        if any_final final!(dfst,cnt,rho) end
        enqueue!(S,qq)
      else
        src, dest = qq2state[pp], qq2state[qq]
        push!(dfst[src],Arc{W,D}(lab,lab,wp,dest))
      end
    end
  end
  return dfst
end


text:
#  github: part of RxNav.jl

const PRESCRIBABLE = [false]  # defaults to larger database

prescribable() = PRESCRIBABLE[1]

""" set whether to use Prescribeable RxNorm database or entire database including items no longer available """
prescribable(tf) = begin PRESCRIBABLE[1] = tf end

""" get the base url for RxNorm, default at loading is the full RxNorm set, not the Prescribable RxNorm set. """
baseurl() = prescribable() ? "prescribable" : "baseurl"

""" RESTuri is a Dict of RxNav REST urls keyed by strings. """
const RESTuri = Dict(
    "baseurl" => "https://rxnav.nlm.nih.gov/REST/",
    "prescribable" => "https://rxnav.nlm.nih.gov/REST/Prescribe/",
    "rcui" => "https://rxnav.nlm.nih.gov/REST/rxcui?name=",
    "drugs" => "https://rxnav.nlm.nih.gov/REST/drugs?name=",
    "interaction" => "https://rxnav.nlm.nih.gov/REST/interaction/interaction?rxcui=",
    "interactionlist" => "https://rxnav.nlm.nih.gov/REST/interaction/list?rxcuis=",
)

"""
    getdoc(urlkey, urltail)

get the document found by the string formed by: (the key indexed by urlkey) * (the urltail)
"""
function getdoc(urlkey, urltail)
    req = HTTP.request("GET", RESTuri[urlkey] * urltail)
    return parsexml(String(req.body)).root
end

"""
    is_in_rxcui_format(s)

Tests whether the string is in the format for an RxNav rxcui identifier.
Currently rxcui identifiers are composed of only digits 0 through 9,
though there is nothing in the schema that says these must be only digits.
If that changes in future RxNav updates, the parsing here may also change.
"""
is_in_rxcui_format(s) = all(c -> c in "0123456789", collect(s))

"""
    morearg(name, value)

Return further arguments to the REST query of the form "&name=value"
or if values is a vector, of the form "&name=val1+val2+val3"
Starts with &, not ?, so there must be a prior argument in the url string.
"""
morearg(name::String, value::String) = "&$name=$(HTTP.URIs.escapeuri(value))"
morearg(name::String, values::Vector{String}) = "&$name=" * join(values, "+")

"""
    morearg(pairs)

Return further arguments to the REST query of the form "&name=value&name2=value2",
etc, where pairs is a vector of [name, value] pairs. Starts with &, not ?,
so there must be a prior argument in the url string.
"""
morearg(pairs::Vector) = reduce(*, morearg(a[1], a[2]) for a in pairs)
morearg(d::Dict) = morearg(collect(d))


text:
function TabulatedLofOffset(arg0::Frame, arg1::LOFType, arg2::List, arg3::jint, arg4::AngularDerivativesFilter)
    return TabulatedLofOffset((Frame, LOFType, List, jint, AngularDerivativesFilter), arg0, arg1, arg2, arg3, arg4)
end

function get_attitude(obj::TabulatedLofOffset, arg0::FieldPVCoordinatesProvider, arg1::FieldAbsoluteDate, arg2::Frame)
    return jcall(obj, "getAttitude", FieldAttitude, (FieldPVCoordinatesProvider, FieldAbsoluteDate, Frame), arg0, arg1, arg2)
end

function get_attitude(obj::TabulatedLofOffset, arg0::PVCoordinatesProvider, arg1::AbsoluteDate, arg2::Frame)
    return jcall(obj, "getAttitude", Attitude, (PVCoordinatesProvider, AbsoluteDate, Frame), arg0, arg1, arg2)
end

function get_table(obj::TabulatedLofOffset)
    return jcall(obj, "getTable", List, ())
end


text:
function _get_all_atc_nodes(graph::PharmGraph)
    all_atc_nodes_set = Vector{PharmClass}(undef, 0)

    function f(class)
        system = class.system::String
        if startswith(system, "ATC")
            push!(all_atc_nodes_set, class)
        end
        return nothing
    end

    _for_each_node(f, graph)

    all_atc_nodes_vector = collect(all_atc_nodes_set)
    unique!(all_atc_nodes_vector)
    sort!(all_atc_nodes_vector)
    return all_atc_nodes_vector
end

function add_atc_hierarchy_edges!(graph::PharmGraph)
    all_atc_levels = _get_all_atc_nodes(graph)

    all_atc1 = Vector{PharmClass}(undef, 0)
    all_atc2 = Vector{PharmClass}(undef, 0)
    all_atc3 = Vector{PharmClass}(undef, 0)
    all_atc4 = Vector{PharmClass}(undef, 0)
    all_atc5 = Vector{PharmClass}(undef, 0)

    for x in all_atc_levels
        system = x.system
        if system == "ATC1"
            push!(all_atc1, x)
        elseif system == "ATC2"
            push!(all_atc2, x)
        elseif system == "ATC3"
            push!(all_atc3, x)
        elseif system == "ATC4"
            push!(all_atc4, x)
        elseif system == "ATC5"
            push!(all_atc5, x)
        end
    end

    unique!(all_atc1)
    unique!(all_atc2)
    unique!(all_atc3)
    unique!(all_atc4)
    unique!(all_atc5)

    sort!(all_atc1)
    sort!(all_atc2)
    sort!(all_atc3)
    sort!(all_atc4)
    sort!(all_atc5)

    atc = [all_atc1, all_atc2, all_atc3, all_atc4, all_atc5]

    for i in 1:4
        atc_i        = atc[i]
        atc_iplusone = atc[i + 1]
        for less_specific_node in atc_i
            for more_specific_node in atc_iplusone
                if startswith(more_specific_node.value, less_specific_node.value)
                    PharmaceuticalClassification._add_edge!(
                        graph,
                        less_specific_node => more_specific_node,
                    )
                end
            end
        end
    end

    return nothing
end


text:
# In v1.0 we should possibly just overload getproperty
# Base.getproperty(w::WrapVar, s::Symbol) = getfield(w.v, s)
struct MaximizationWrapper{T}
    res::T
end
res(r::MaximizationWrapper) = r.res

# ==============================================================================
# Univariate warppers
# ==============================================================================
function maximize(f, lb::Real, ub::Real, method::AbstractOptimizer; kwargs...)
    fmax = x->-f(x)
    MaximizationWrapper(optimize(fmax, lb, ub, method; kwargs...))
end

function maximize(f, lb::Real, ub::Real; kwargs...)
    fmax = x->-f(x)
    MaximizationWrapper(optimize(fmax, lb, ub; kwargs...))
end

# ==============================================================================
# Multivariate warppers
# ==============================================================================
function maximize(f, x0::AbstractArray; kwargs...)
    fmax = x->-f(x)
    MaximizationWrapper(optimize(fmax, x0; kwargs...))
end
function maximize(f, x0::AbstractArray, method::AbstractOptimizer, options = Optim.Options(); kwargs...)
    fmax = x->-f(x)
    MaximizationWrapper(optimize(fmax, x0, method, options; kwargs...))
end
function maximize(f, g, x0::AbstractArray, method::AbstractOptimizer, options = Optim.Options(); kwargs...)
    fmax = x->-f(x)
    gmax = (G,x)->(g(G,x); G.=-G)
    MaximizationWrapper(optimize(fmax, gmax, x0, method, options; kwargs...))
end

function maximize(f, g, h, x0::AbstractArray, method::AbstractOptimizer, options = Optim.Options(); kwargs...)
    fmax = x->-f(x)
    gmax = (G,x)->(g(G,x); G.=-G)
    hmax = (H,x)->(h(G,x); H.=-H)
    MaximizationWrapper(optimize(fmax, gmax, hmax, x0, method, options; kwargs...))
end

minimum(r::MaximizationWrapper) = throw(MethodError())
maximizer(r::Union{UnivariateOptimizationResults,MultivariateOptimizationResults}) = throw(MethodError())
maximizer(r::MaximizationWrapper) = minimizer(res(r))
maximum(r::Union{UnivariateOptimizationResults,MultivariateOptimizationResults}) = throw(MethodError())
maximum(r::MaximizationWrapper) = -minimum(res(r))
Base.summary(r::MaximizationWrapper) = summary(res(r))

for api_method in (:lower_bound, :upper_bound, :rel_tol, :abs_tol, :iterations, :initial_state, :converged, :x_tol, :x_converged,
               :x_abschange, :g_tol, :g_converged, :g_residual, :f_tol, :f_converged,
               :f_increased, :f_relchange, :iteration_limit_reached, :f_calls,
               :g_calls, :h_calls)
   @eval $api_method(r::MaximizationWrapper) = $api_method(res(r))
end

function Base.show(io::IO, r::MaximizationWrapper{<:UnivariateOptimizationResults})
    @printf io "Results of Maximization Algorithm\n"
    @printf io " * Algorithm: %s\n" summary(r)
    @printf io " * Search Interval: [%f, %f]\n" lower_bound(r) upper_bound(r)
    @printf io " * Maximizer: %e\n" maximizer(r)
    @printf io " * Maximum: %e\n" maximum(r)
    @printf io " * Iterations: %d\n" iterations(r)
    @printf io " * Convergence: max(|x - x_upper|, |x - x_lower|) <= 2*(%.1e*|x|+%.1e): %s\n" rel_tol(r) abs_tol(r) converged(r)
    @printf io " * Objective Function Calls: %d" f_calls(r)
    return
end

function Base.show(io::IO, r::MaximizationWrapper{<:MultivariateOptimizationResults})
    take = Iterators.take

    @printf io "Results of Optimization Algorithm\n"
    @printf io " * Algorithm: %s\n" summary(r.res)
    if length(join(initial_state(r), ",")) < 40
        @printf io " * Starting Point: [%s]\n" join(initial_state(r), ",")
    else
        @printf io " * Starting Point: [%s, ...]\n" join(take(initial_state(r),
2), ",")
    end
    if length(join(maximizer(r), ",")) < 40
        @printf io " * Maximizer: [%s]\n" join(maximizer(r), ",")
    else
        @printf io " * Maximizer: [%s, ...]\n" join(take(maximizer(r), 2), ",")
    end
    @printf io " * Maximum: %e\n" maximum(r)
    @printf io " * Iterations: %d\n" iterations(r)
    @printf io " * Convergence: %s\n" converged(r)
    if isa(r.res.method, NelderMead)
        @printf io "   *  √(Σ(yᵢ-ȳ)²)/n < %.1e: %s\n" g_tol(r) g_converged(r)
    else
        @printf io "   * |x - x'| ≤ %.1e: %s \n" x_tol(r) x_converged(r)
        @printf io "     |x - x'| = %.2e \n"  x_abschange(r)
        @printf io "   * |f(x) - f(x')| ≤ %.1e |f(x)|: %s\n" f_tol(r) f_converged(r)
        @printf io "     |f(x) - f(x')| = %.2e |f(x)|\n" f_relchange(r)
        @printf io "   * |g(x)| ≤ %.1e: %s \n" g_tol(r) g_converged(r)
        @printf io "     |g(x)| = %.2e \n"  g_residual(r)
        @printf io "   * Stopped by an decreasing objective: %s\n" (f_increased(r) && !iteration_limit_reached(r))
    end
    @printf io "   * Reached Maximum Number of Iterations: %s\n" iteration_limit_reached(r)
    @printf io " * Objective Calls: %d" f_calls(r)
    if !(isa(r.res.method, NelderMead) || isa(r.res.method, SimulatedAnnealing))
        @printf io "\n * Gradient Calls: %d" g_calls(r)
    end
    if isa(r.res.method, Newton) || isa(r.res.method, NewtonTrustRegion)
        @printf io "\n * Hessian Calls: %d" h_calls(r)
    end
    return
end


text:
### Compare LombScargle.jl with Astropy

using LombScargle, Test, PyCall, Random
ast = pyimport("astropy.timeseries")

Random.seed!(1)
ntimes = 401

@testset "Un-evenly spaced data" begin
    t = range(0.01, stop=10pi, length=ntimes)
    t += step(t)*rand(ntimes)
    for f in (x -> sinpi(x), x -> sin(x) + 1.5*cospi(4*x) + 3)
        s = f.(t)
        # "psd" normalization in LombScargle slightly differs from that of
        # Astropy in a few points and the test would fail if we includ it.
        @testset "$fitmean, $nrm, $fast, $center" for fitmean in (true, false),
            nrm in ("standard", "model", "log"), fast in ((true, "fast"), (false, "cython")),
            center in (true, false)
            f_jl, p_jl = freqpower(lombscargle(t, s, fit_mean = fitmean,
                                               center_data = center,
                                               normalization=Symbol(nrm),
                                               padding_factors=[2],
                                               maximum_frequency=20, fast = fast[1]))
            f_py, p_py =
                ast.LombScargle(t, s,
                                fit_mean = fitmean,
                                center_data = center).autopower(method=fast[2],
                                                                normalization=nrm,
                                                                maximum_frequency=20)
            @test f_jl ≈ f_py
            @test p_jl ≈ p_py
        end
    end
end

t = range(0.01, stop=10pi, length=ntimes)
errors = rand(0.1:1e-3:4.0, ntimes)

@testset "Evenly spaced data" begin # Use both heteroskedastic and homoskedastic uncertainties.
    for f in (x -> sinpi(x), x -> sin(x) + 1.5*cospi(4*x) + 3), err in (ones(ntimes), errors)
        s = f.(t)
        @testset "$fitmean, $nrm, $fast, $center" for fitmean in (true, false),
            nrm in ("standard", "model", "log", "psd"), fast in ((true, "fast"), (false, "cython")),
            center in (true, false)
            f_jl, p_jl = freqpower(lombscargle(t, s, err,
                                               fast = fast[1],
                                               fit_mean = fitmean,
                                               center_data = center,
                                               normalization=Symbol(nrm),
                                               padding_factors=[2],
                                               maximum_frequency=10,
                                               samples_per_peak=10))
            f_py, p_py =
                ast.LombScargle(t, s, dy = err,
                                fit_mean = fitmean,
                                center_data = center).autopower(method=fast[2],
                                                                normalization=nrm,
                                                                maximum_frequency=10,
                                                                samples_per_peak=10)
            @test f_jl ≈ f_py
            @test p_jl ≈ p_py
        end
    end
end

@testset "heteroskedastic uncertainties" begin #  with non-fast method
    for f in (x -> sinpi(x), x -> sin(x) + 1.5*cospi(4*x) + 3)
        s = f.(t)
        @testset "$fitmean, $nrm, $fast, $center" for fitmean in (true, false),
            nrm in ("standard", "model", "psd"), fast in ((true, "fast"), (false, "cython")),
            center in (true, false)
            f_jl, p_jl = freqpower(lombscargle(t, s, errors,
                                               fit_mean = fitmean,
                                               center_data = center,
                                               fast = fast[1],
                                               normalization = Symbol(nrm),
                                               padding_factors=[2],
                                               maximum_frequency=20))
            f_py, p_py =
                ast.LombScargle(t, s, dy = errors,
                                fit_mean = fitmean,
                                center_data = center).autopower(method=fast[2],
                                                                normalization = nrm,
                                                                maximum_frequency=20)
            @test f_jl ≈ f_py
            @test p_jl ≈ p_py
        end
    end
end

@testset "Model functions" begin
    for f in (x -> sinpi(x), x -> sin(x) + 1.5*cospi(4*x) + 3)
        s = f.(t)
        @testset "$fm, $cd" for fm in (true, false), cd in (true, false)
            m_jl = LombScargle.model(t, s, 1/2pi, fit_mean=fm, center_data=cd)
            m_py = ast.LombScargle(t, s, center_data=cd, fit_mean=fm).model(t, 1/2pi)
            @test m_jl ≈ m_py
        end
    end
end

@info("LombScargle.jl is consistent with Astropy")


text:
using FinancialToolbox, DualNumbers, HyperDualNumbers, FinancialFFT, FinancialMonteCarlo

A = 600.0;
N = 18;

method_cm = CarrMadanMethod(A, N);
S0 = 100.0
K = 100.0;
r = hyper(0.02, 1.0, 1.0, 0.0);

T = 1.0;
d = 0.01
sigma = 0.2
lam = 15.0;
mu1 = 0.03;
sigma1 = 0.02;
zero_rate = ZeroRate(r);

Model = VarianceGammaProcess(sigma, 0.03, 0.01, Underlying(S0, d))

EUData = EuropeanOption(T, K)
Nsim = 100000;
Nstep = 30;
function test_hyper(x, y, toll)
    @test abs(x.value - y.value) < toll
    @test abs(x.epsilon1 - y.epsilon1) < toll
    @test abs(x.epsilon2 - y.epsilon2) < toll
    @test abs(x.epsilon12 - y.epsilon12) < toll
    nothing
end

result_cm = pricer(Model, zero_rate, method_cm, EUData);
method_lewis = LewisMethod(700.0, 200000);
analytic_result = pricer(Model, zero_rate, method_lewis, EUData);
method_cm_lewis = CarrMadanLewisMethod(A, N);
result_cm_lewis = pricer(Model, zero_rate, method_cm_lewis, EUData);
toll = 1e-3;
#test_hyper(analytic_result, result_cm, toll)
@testset "hyper test" begin
    test_hyper(analytic_result, result_cm_lewis, toll)
end


text:
module Memoize

export @memoize, @dumpf, @showf

include("./Parser.jl")

macro memoize(f::Expr)
    esc(f_expr(f))
end

macro dumpf(f::Expr)
    dump(f_expr(f), maxdepth=32)
end

macro showf(f::Expr)
    show(f_expr(f))
end

end


text:
export @mul

#= TODO:

# add @check! version

# anonymous functions, and scalar output, for consistency?

=#

"""
    @mul A[i,j] := B[i,k,k′] * C[k,k′,j]
    @mul A[i,j] := sum(k,k′) B[i,k,k′] * C[k,k′,j]

Matrix multiplication macro. This expects two tensors on the right,
whose shared index (or indices) will be sumed over.
You may also explicitly supply these, as for `@reduce`, as a check.

    @mul A[_,i][j] := B[i\\k] * C[j,2,k]
    @mul A[_,i][j] := sum(k) B[i\\k] * C[j,2,k]

Each tensor factor will be processed in the same way as for `@cast` / `@reduce`,
allowing for slicing and reshaping, permuting indices, and fixing their values.
(The same tuple of options understood by `@cast` may be given after the expression, too.)
But the right hand side must be simply a product of two tensors, nothing else.
Thus it is both more limited, and more general, than `@tensor` / `@einsum`.

The same operations can be written `@reduce A[i,j] := sum(k) B[i,k] * C[k,j]`,
but this broadcasts out to a 3-tensor before reducing,
rather than just calling `A = B * C`.

The in-place form `@mul A[i,j] = B[i,k] * C[k,j]` calls instead `mul!(A,B,C)`.

    @mul A[i,j,n] := B[i,k,n] * C[k,j,n]
    @mul A[i,j,n] := sum(k) B[i,k,n] * C[k,j,n]

Batched matrix multiplication. Right now done by slicing each tensor and mapping `*` or `mul!`
over the slices. But it should be easy to hook up to Batched.jl or something.
"""
macro mul(exs...)
    where = (mod=__module__, src=__source__, str=unparse("@mul", exs...))
    _mul(exs...; icheck=false, where=where)
end

#=
"""
    @mul! A[i,j] := B[i,k] * C[k,j]

Variant of `@mul` which effectively runs `@check!()` on each tensor.
"""
macro mul!(exs...)
    where = (mod=__module__, src=__source__, str=unparse("@mul", exs...))
    _mul(exs...; icheck=true, where=where)
end
=#

function _mul(exone, extwo=nothing, exthree=nothing; icheck=false, where=where, recurse=false)
    flags = Set{Symbol}()

    #===== parse basic expression =====#

    if @capture(exone, left_ = mid_ * right_ )
        push!(flags, :inplace)
    elseif @capture(exone, left_ := mid_ * right_ )
    elseif @capture(exone, left_ |= mid_ * right_ )
        push!(flags, :mustcopy) # this will only affect output slicing I think

    elseif @capture(exone, left_ = sum(sind__) ) &&  @capture(extwo, mid_ * right_ )
        push!(flags, :inplace)
        push!(flags, :explicit)
    elseif @capture(exone, left_ := sum(sind__) ) &&  @capture(extwo, mid_ * right_ )
        push!(flags, :explicit)
    elseif @capture(exone, left_ |= sum(sind__) ) &&  @capture(extwo, mid_ * right_ )
        push!(flags, :mustcopy)
        push!(flags, :explicit)

    else
        throw(MacroError("don't understand input, should be @mul A[...] := B[...] * C[...]", where))
    end

    if :explicit in flags
        options = exthree
    else
        options = extwo
        exthree == nothing || throw(MacroError("don't understand input, should be @mul A[...] := B[...] * C[...]  options", where))
    end

    #===== find indices to sum, batch, etc =====#

    store = SizeDict()

    # parse options both to look for keywords and sizes
    @capture(options, (optvec__,)) || (optvec = Any[options])
    optind, _,_,_ = parse!(store, nothing, [], optvec; allowranges=true, flags=flags)

    # get list of indices from each factor
    indZ, nameZ = firstpass!(store, left, where, :inplace in flags) #; allowrecursion=false)
    indA, nameA = firstpass!(store, mid, where)
    indB, nameB = firstpass!(store, right, where)

    if nameZ == nothing
        :inplace in flags && throw(MacroError("can't write in-place to an un-named tensor",where))
        nameZ == gensym(:Z)
    end
    nameA == nothing || nameB == nothing && throw(MacroError("can't work without tensor names",where))
    checkrepeats(indZ, " on left hand side", where)
    checkrepeats(indA, " in first factor", where)
    checkrepeats(indB, " in second factor", where)

    sumind = setdiff(intersect(indA, indB), indZ)
    batchind = intersect(indA, indB, indZ)

    if :explicit in flags
        sort(sumind) == sort(sind) || throw(MacroError("explicit sum over $(join(sind, ", ")) does not match implicit $(join(sumind, ", "))", where))
    end

    if length(batchind) > 0
        push!(flags, :bmm)
    end
    if length(sumind)>1 || length(batchind)>1 # then result of * (or batchmul) is wrong shape
        push!(flags, :outshape)
    end

    indAonly = setdiff(indA, sumind, batchind)
    indBonly = setdiff(indB, sumind, batchind)
    indZonly = setdiff(indZ, indAonly, indBonly, batchind)
    length(indZonly) == 0 || throw(MacroError("did not see index $(join(indZonly, ", ")) on the right", where))

    if length(indAonly) > 1 || length(indBonly) > 1 || length(batchind) > 1
        push!(flags, :midshape) # result of * will not be have canonical dimensions
    end

    canon = vcat(indZ, sumind)

    V && @info "@mul indices" Tuple(indAonly) Tuple(indBonly) Tuple(canon) Tuple(sumind) Tuple(batchind)

    # having got canonical list, we can check the options
    if count(i -> i != nothing, setdiff(optind, canon)) > 0
        str = join(something.(setdiff(optind, canon), "nothing"), ", ")
        m_error("attempting to ignore unrecognised options: $str", where)
    end

    #===== process each factor in detail =====#

    outex = MacroTools.@q begin end

    if :nolazy in flags
        m_error("@mul always ignores option nolazy, try |= to copy output", where)
    end
    push!(flags, :nolazy) # exA & exB should use permutedims

    exA = matexin!(outex, nameA, mid, (indAonly, sumind, batchind), flags, store, icheck, where)
    exB = matexin!(outex, nameB, right, (sumind, indBonly, batchind), flags, store, icheck, where)

    if :inplace in flags
        pop!(flags, :nolazy) # exZ must not use permutedims

        exZ = matexin!(outex, nameZ, left, (indAonly, indBonly, batchind), flags, store, icheck, where)

        @capture(left, _[__][__]) && throw(MacroError("can't write to sliced arrays in-place, for now", where))
    else
        # exZ will not be used for in-place mul!, but for out-of-place... can borrow readleft from @cast:
        _, outUZ, nameZ, checkZ = readleft(left, sumind, flags, store, icheck, where)

        # but unlike @cast/@reduce, the RHS does not arrive in canonical order...
        V && @info "more" Tuple(indZ) Tuple(vcat(indAonly,indBonly,batchind))

        p1 = sortperm(indZ)
        p2 = sortperm(vcat(indAonly,indBonly,batchind))
        permU = p2[invperm(p1)]
        if permU != 1:length(permU)
            push!(flags, :outperm)
            permUold = invperm(p1)[p2] # not one of my old tests triggered this!
            V && permU != permUold && @warn "permutations differ in $where.str" Tuple(p1) Tuple(p2) Tuple(permU) Tuple(permUold)
        end
    end

    packagecheck(flags, where)

    canonsize = sizeinfer(store, canon, where, true)

    if :inplace in flags

        #===== in-place output =====#

        if :batched in flags && :bmm in flags
            error("can't use batched yet!")
        elseif :bmm in flags
            if length(indAonly)==0 # vectors * something
                opex =  :( TensorCast.batchmul!(
                    TensorCast.orient($exZ,(*,:,:)), TensorCast.orient($exA,(*,:,:)), $exB))
                # TODO if z::Transpose then orient() may copy here, issue a warning?
            else
                opex =  :( TensorCast.batchmul!($exZ, $exA, $exB))
            end
        else
            if length(indAonly)==0 # vector * something
                opex = :( TensorCast.mul!(
                    TensorCast.orient($exZ,(*,:)), TensorCast.orient($exA,(*,:)), $exB))
            else
                opex = :( TensorCast.mul!($exZ, $exA, $exB))  # LinearAlgebra.mul! is available here
            end
        end

        if :strided in flags # this functions mostly to prevent a slowdown if "inputex" has used strided
            opex = :( Strided.@strided $opex )
        end

        push!(outex.args, opex )
        push!(outex.args, nameZ )

    else

        #===== out-of-place output =====#

        if :batched in flags && :bmm in flags
            error("can't use batched yet!")
        elseif :bmm in flags
            if length(indAonly)==0 # vectors * something
                newright = :( TensorCast.batchmul(TensorCast.orient($exA,(*,:,:)), $exB) )
            else
                newright = :( TensorCast.batchmul($exA, $exB) )
            end
        else
            if length(indAonly)==0 # vector * something
                if :scalar in flags
                    length(indBonly)==0 || m_error("can't make scalar output from this...")
                    newright = :( TensorCast.orient($exA,(*,:)) * $exB )
                else
                    newright = :( TensorCast.rvec( TensorCast.orient($exA,(*,:)) * $exB ))
                    # rvec handles vector * vector here, else :midshape would be cleaner
                end
            else
                newright = :( $exA * $exB )
            end
        end

        if :midshape in flags
            sizeM = Any[ Symbol(:sz_, i) for i in vcat(indAonly,indBonly,batchind) ]
            sizeMex = :(($(sizeM...) ,))
            newright = :( reshape( $newright, $sizeMex ) )
        end

        if :outperm in flags
            if :strided in flags
                ex = :( strided_permutedims($ex, $perm) )
            elseif permU == [2,1]
                newright = :( TensorCast.PermuteDims($newright) )
            else
                newright = :( PermutedDimsArray($newright, ($(permU...),) ))
            end
        end

        if :strided in flags # this functions mostly to prevent a slowdown if "inputex" has used strided
            newright = :( Strided.@strided $newright )
        end

        finalright = outputnew(newright, outUZ, nothing, canonsize, canon, flags, store, where)

        push!(outex.args, :( $nameZ =  $finalright ) )
    end

    #===== finalise (almost identical to @cast) =====#

    if :needsize in flags
        szcanon = Any[ Symbol(:sz_,i) for i in canon ]
        pushfirst!(outex.args, :( local ($(szcanon...),) = ($(canonsize...),) ) )
    end
    if :assert in flags || :(!) in flags || check_options.size
        for ch in store.checks
            pushfirst!(outex.args, ch)
        end
    end
    for tex in store.topex
        pushfirst!(outex.args, tex)
    end

    if recurse==true # only for @mul inside of something
        indZ = outUZ[end-1]
        return outex, indZ
    end

    if length(outex.args) == 1
        return esc(outex.args[1])
    else
        return esc(outex)
    end
end

#="""
    indA, A = firstpass(store, ex)

Aim is just to get the list of indices for each piece, to figure out which ones are contracted,
and to have canonical order (now can't be got from one side alone).
"""=#
function firstpass!(sdict, ex, where, canseeA=true) #; allowrecursion=true)

    if @capture(ex, A_[outer__][inner__]) ||  @capture(ex, [outer__][inner__])
    elseif @capture(ex, A_[outer__]{inner__}) || @capture(ex, [outer__]{inner__})
    elseif @capture(ex, A_[outer__]) || @capture(ex, [outer__])
        inner = []

    #=
    # allow @reduce/@mul inside RHS, here just believe the indices
    elseif allowrecursion && @capture(ex, @reduce(redex__)) || @capture(ex, @mul(redex__))
        inner = []
        @capture(redex[1], A_[outer__] = B__) ||
        @capture(redex[1], A_[outer__] := B__) ||
        @capture(redex[1], A_[outer__] |= B__) ||
        @capture(redex[1], [outer__] = B__) ||
        @capture(redex[1], [outer__] := B__) ||
        @capture(redex[1], [outer__] |= B__) ||
            throw(MacroError("recursion got confused about $(redex[1])", where))
    # But this needs you to copy lots more from "walker" to "matexin!"
    =#
    elseif @capture(ex, @reduce(redex__)) || @capture(ex, @mul(redex__))
        throw(MacroError("@mul does not allow recursion, sorry", where))
    else
        throw(MacroError("don't know what to do with $ex", where))
    end

    if canseeA

        # if we have f(x)[i,j] then we should evaluate f just once
        if !isa(A, Symbol)
            Atop = gensym(:A)
            push!(sdict.topex, :(local $Atop = $A ) )
            A = Atop
        end

        flat, getafix, putsz, negated = parse!(sdict, A, outer, inner)
    else
        flat, getafix, putsz, negated = parse!(sdict, nothing, outer, inner)
    end
    flat, A
end

#="""
    matexin!(outex, A, ex, (ind1,ind2,ind3), flags, store, icheck, where)

Produces the expression which will glue/permutedims etc A to form `ind1,ind2,ind3` (as `@cast`),
and then reshape A to be either a matrix (indices `⨷ind1, ⨷ind2`)
or for the batched case, a 3-tensor (indices `⨷ind1, ⨷ind2, ⨷ind3`).
Pushes this into outex & returns the symbol it is assigned to.
"""=#
function matexin!(outex, A, ex, (ind1,ind2,ind3), flags, store, icheck, where)

    # same function called by @cast walker -- produces reshapes to 2nd input, was canon
    Aval = inputex(A, ex, vcat(ind1,ind2,ind3), flags, SizeDict(), icheck, where)
    # @pretty @mul A[n][i] |= B[k,i,n] * C[k\n] # TODO fix this, inputex has assumptions...

    # one more processing step: reshape to matrix or batch-tensor
    if length(ind1)<=1 && length(ind2)<=1 && length(ind3) <= 1
    else
        prodsA = prodsizeex(ind1,ind2,ind3)
        Aval = :( reshape($Aval, ($(prodsA...),)) )
        push!(flags, :needsize) # resize needs sz_i
    end

    # push pre-processing step into outex if nontrivial, as walker did
    if isa(Aval, Symbol)
        ex = Aval
    else
        Asym = gensym(:A)
        push!(outex.args, :(local $Asym = $Aval) )
        ex =  Asym
    end

    # return what should appear in mul! etc.
    return ex
end

function prodsizeex(ind123...)
    out = Any[]
    for ind in ind123
        if length(ind) == 1
            push!(out, szwrap(ind[1]) )
        elseif length(ind) >= 2
            push!(out, szwrap(ind) )
        end
    end
    out
end

using LinearAlgebra

"""
    batchmul(A, B)      # @mul Z[i,j,n] := A[i,k,n] * B[k,j,n]
    batchmul!(Z, A, B)  # @mul Z[i,j,n] = A[i,k,n] * B[k,j,n]

Batched matrix multiplication; batch index is the last dimension.
"""
function batchmul!(Z::AbstractArray{TZ,N}, A::AbstractArray{TA,3}, B::AbstractArray{TB,N}) where {TA,TZ,TB,N}
    2 <= N <= 3 || throw(DimensionMismatch("expected ndims = 3,3,3 or else 2,3,2, got ndims = $N,3,$N"))
    size(Z,N) == size(A,3) == size(B,N) || throw(DimensionMismatch("arrays must agree on batch dimension"))
    codeB = ntuple(i-> i==N ? (*) : (:), N)
    sZ = sliceview(Z, codeB)
    sA = sliceview(A, (:,:,*))
    sB = sliceview(B, codeB)
    foreach(mul!, sZ, sA, sB)
    Z
end

@doc @doc(batchmul!)
function batchmul(A::AbstractArray{TA,3}, B::AbstractArray{TB,N}) where {TA,TB,N}
    2 <= N <= 3 || throw(DimensionMismatch("expected ndims = 3,3 or else 3,2, got ndims = 3,$N"))
    size(A,3) == size(B,N) || throw(DimensionMismatch("arrays must agree on batch dimension"))
    codeB = ntuple(i-> i==N ? (*) : (:), N)
    sA = sliceview(A, (:,:,*))
    sB = sliceview(B, codeB)
    sZ = map(*, sA, sB)
    red_glue(sZ, codeB)
end

# https://github.com/FluxML/Flux.jl/issues/544
# https://github.com/Roger-luo/BatchedRoutines.jl
# https://github.com/Roger-luo/Batched.jl


text:
using BinaryBuilder

name = "libmseed"
version = v"3.0.10"

# Collection of sources required to build libmseed
sources = [
    GitSource("https://github.com/iris-edu/libmseed",
              "d8839780299b06fb09c6f1961c293f93c3d0427e"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libmseed/

# Specify LIBDIR and LIB_SO_BASE since these are different on Windows than otherwise
make install PREFIX="${prefix}" LIBDIR="${libdir}" LIB_SO_BASE=libmseed.${dlext}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libmseed", :libmseed),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
"""
	gaussian_path(n::Int64 = 4)

Returns a 'nx1' dimensional array of Gaussian random walk.
"""
function gaussian_path(n::Int64 = 4)
    return vcat(0.0, cumsum(randn(rng, n-1, 1), dims = 1)) #n stages
end

"""
	running_maximum(n::Int64 = 4)

Returns a 'nx1' dimensional array of Running Maximum process.
"""
function running_maximum(n::Int64 = 4)
    rmatrix = vcat(0.0, cumsum(randn(rng, 3, 1), dims = 1))
    for i = 2 : 4
        rmatrix[i] = max.(rmatrix[i-1], rmatrix[i])
    end
    return rmatrix
end


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

# example from http://openfst.org/twiki/bin/view/FST/RmEpsilonDoc
W = TropicalWeight{Float64}
A = WFST(Dict("a"=>1,"p"=>2),W = W)
add_arc!(A,1,2,"<eps>","<eps>",1)
add_arc!(A,2,3,"a","<eps>",2)
add_arc!(A,2,3,"<eps>","p",3)
add_arc!(A,2,3,"<eps>","<eps>",4)
add_arc!(A,3,3,"<eps>","<eps>",5)
add_arc!(A,3,4,"<eps>","<eps>",6)
final!(A,4,7)
initial!(A,1)

# println(A)
B = rm_eps(A)
# println(B)
@test isinitial(B,1)
@test get_final_weight(B,1) == W(18)
@test get_initial_weight(B,1) == one(W)
@test isfinal(B,1)
@test isfinal(B,2)
@test get_final_weight(B,2) == W(13)
@test B[1][1] == FiniteStateTransducers.Arc{W,Int}(1,0,W(3),2)
@test B[1][2] == FiniteStateTransducers.Arc{W,Int}(0,2,W(4),2)
@test size(B) == (2,2)