text:
using ITensors, Test

@testset "Fermions" begin
  ITensors.enable_auto_fermion()

  @testset "parity_sign function" begin

    # Full permutations
    p1 = [1, 2, 3]
    @test ITensors.parity_sign(p1) == +1
    p2 = [2, 1, 3]
    @test ITensors.parity_sign(p2) == -1
    p3 = [2, 3, 1]
    @test ITensors.parity_sign(p3) == +1
    p4 = [3, 2, 1]
    @test ITensors.parity_sign(p4) == -1

    ## Partial permutations
    p5 = [2, 7]
    @test ITensors.parity_sign(p5) == +1
    p6 = [5, 3]
    @test ITensors.parity_sign(p6) == -1
    p7 = [1, 9, 3, 10]
    @test ITensors.parity_sign(p7) == -1
    p8 = [1, 12, 9, 3, 11]
    @test ITensors.parity_sign(p8) == +1
  end

  @testset "Fermionic QNs" begin
    q = QN("Nf", 1, -1)
    @test isfermionic(q[1])
    @test fparity(q) == 1

    q = q + q + q
    @test val(q, "Nf") == 3

    p = QN("P", 1, -2)
    @test fparity(p) == 1
    @test fparity(p + p) == 0
    @test fparity(p + p + p) == 1
  end

  @testset "Fermionic IndexVals" begin
    sn = Index([QN("Nf", 0, -1) => 1, QN("Nf", 1, -1) => 1], "sn")
    @test fparity(sn => 1) == 0
    @test fparity(sn => 2) == 1

    sp = Index([QN("Nfp", 0, -2) => 1, QN("Nfp", 1, -2) => 1], "sp")
    @test fparity(sp => 1) == 0
    @test fparity(sp => 2) == 1
  end

  @testset "Get and Set Elements" begin
    s = Index([QN("Nf", 0, -1) => 1, QN("Nf", 1, -1) => 1], "s")

    N = ITensor(s', dag(s))

    N[s' => 2, s => 2] = 1.0
    @test N[s' => 2, s => 2] ≈ +1.0
    @test N[s => 2, s' => 2] ≈ -1.0

    N[s => 2, s' => 2] = 1.0
    @test N[s' => 2, s => 2] ≈ -1.0
    @test N[s => 2, s' => 2] ≈ 1.0

    C = ITensor(s', dag(s))

    C[s' => 1, s => 2] = 1.0
    @test C[s' => 1, s => 2] ≈ 1.0
    @test C[s => 2, s' => 1] ≈ 1.0

    I = ITensor(s', dag(s))
    I[s' => 1, s => 1] = 1.0
    I[s' => 2, s => 2] = 1.0
    @test I[s' => 1, s => 1] ≈ 1.0
    @test I[s' => 2, s => 2] ≈ 1.0

    @test I[s => 1, s' => 1] ≈ 1.0
    @test I[s => 2, s' => 2] ≈ -1.0
  end

  @testset "Making operators different ways" begin
    s = Index([QN("Nf", 0, -1) => 1, QN("Nf", 1, -1) => 1], "s")

    N1 = ITensor(s', dag(s))
    N1[s' => 2, s => 2] = +1.0

    N2 = ITensor(dag(s), s')
    N2[s' => 2, s => 2] = +1.0
    @test norm(N1 - N2) ≈ 0.0

    N3 = ITensor(s', dag(s))
    N3[s => 2, s' => 2] = -1.0
    @test norm(N1 - N3) ≈ 0.0

    N4 = ITensor(dag(s), s')
    N4[s => 2, s' => 2] = -1.0
    @test norm(N1 - N4) ≈ 0.0
  end

  @testset "Permute and Add Fermionic ITensors" begin
    @testset "Permute Operators" begin
      s = Index([QN("Nf", 0, -1) => 1, QN("Nf", 1, -1) => 1], "s")

      N1 = ITensor(s', dag(s))
      N1[s' => 2, s => 2] = 1.0

      N2 = ITensor(dag(s), s')
      N2[s' => 2, s => 2] = 1.0

      pN1 = permute(N1, dag(s), s')
      @test pN1[s' => 2, s => 2] ≈ 1.0

      pN2 = permute(N2, s', dag(s))
      @test pN2[s' => 2, s => 2] ≈ 1.0

      #TODO add cases resulting in minus signs
    end

    @testset "Add Operators" begin
      s = Index([QN("Nf", 0, -1) => 1, QN("Nf", 1, -1) => 1], "sn")

      N1 = ITensor(s', dag(s))
      N1[s' => 2, s => 2] = 1.0

      N2 = ITensor(dag(s), s')
      N2[s' => 2, s => 2] = 1.0

      NN = N1 + N2
      @test NN[s' => 2, s => 2] ≈ 2.0

      NN = N1 + N1
      @test NN[s' => 2, s => 2] ≈ 2.0

      NN = N2 + N2
      @test NN[s' => 2, s => 2] ≈ 2.0
    end

    @testset "Wavefunction Tests" begin
      s = [Index([QN("N", 0, -2) => 2, QN("N", 1, -2) => 2], "s$n") for n in 1:4]

      psi0 = ITensor(s...)

      psi0[s[1] => 1, s[2] => 1, s[3] => 1, s[4] => 1] = 1111
      psi0[s[1] => 3, s[2] => 3, s[3] => 1, s[4] => 1] = 3311
      psi0[s[1] => 1, s[2] => 3, s[3] => 1, s[4] => 3] = 1313

      psi1 = permute(psi0, s[2], s[1], s[3], s[4])
      @test norm(psi1 - psi0) ≈ 0.0

      @test psi0[s[1] => 1, s[2] => 1, s[3] => 1, s[4] => 1] ≈ 1111
      @test psi1[s[1] => 1, s[2] => 1, s[3] => 1, s[4] => 1] ≈ 1111
      @test psi0[s[2] => 1, s[1] => 1, s[3] => 1, s[4] => 1] ≈ 1111
      @test psi1[s[2] => 1, s[1] => 1, s[3] => 1, s[4] => 1] ≈ 1111

      @test psi0[s[1] => 3, s[2] => 3, s[3] => 1, s[4] => 1] ≈ 3311
      @test psi1[s[1] => 3, s[2] => 3, s[3] => 1, s[4] => 1] ≈ 3311
      @test psi0[s[2] => 3, s[1] => 3, s[3] => 1, s[4] => 1] ≈ -3311
      @test psi1[s[2] => 3, s[1] => 3, s[3] => 1, s[4] => 1] ≈ -3311
      @test psi0[s[4] => 1, s[2] => 3, s[1] => 3, s[3] => 1] ≈ -3311
      @test psi1[s[4] => 1, s[2] => 3, s[1] => 3, s[3] => 1] ≈ -3311

      psi2 = permute(psi0, s[4], s[1], s[3], s[2])
      @test norm(psi2 - psi0) ≈ 0.0
      @test norm(psi2 - psi1) ≈ 0.0

      @test psi0[s[1] => 1, s[2] => 3, s[3] => 1, s[4] => 3] ≈ 1313
      @test psi1[s[1] => 1, s[2] => 3, s[3] => 1, s[4] => 3] ≈ 1313
      @test psi2[s[1] => 1, s[2] => 3, s[3] => 1, s[4] => 3] ≈ 1313
      @test psi0[s[4] => 3, s[1] => 1, s[3] => 1, s[2] => 3] ≈ -1313
      @test psi1[s[4] => 3, s[1] => 1, s[3] => 1, s[2] => 3] ≈ -1313
      @test psi2[s[4] => 3, s[1] => 1, s[3] => 1, s[2] => 3] ≈ -1313
    end
  end

  @testset "C Cdag operators" begin
    s = siteinds("Fermion", 3; conserve_qns=true)

    p110 = ITensor(s[1], s[2], s[3])
    p110[s[1] => 2, s[2] => 2, s[3] => 1] = 1.0

    p011 = ITensor(s[1], s[2], s[3])
    p011[s[1] => 1, s[2] => 2, s[3] => 2] = 1.0

    np011 = ITensor(s[1], s[2], s[3])
    np011[s[1] => 1, s[3] => 2, s[2] => 2] = 1.0

    dag_p011 = ITensor(dag(s[3]), dag(s[2]), dag(s[1]))
    dag_p011[s[3] => 2, s[2] => 2, s[1] => 1] = 1.0

    @test norm(dag(p011) - dag_p011) ≈ 0

    C1 = op(s, "C", 1)
    Cdag3 = op(s, "Cdag", 3)

    # Anti-commutator
    @test norm(Cdag3 * C1 + C1 * Cdag3) ≈ 0.0

    # Commutator
    @test norm(Cdag3 * C1 - C1 * Cdag3) ≈ 2.0

    let # <011|Cdag3*C1|110> = -1
      t1 = noprime(C1 * p110)
      t2 = noprime(Cdag3 * t1)
      @test scalar(dag_p011 * t2) ≈ -1.0
    end

    let # <011|C1*Cdag3|110> = +1
      t1 = noprime(Cdag3 * p110)
      t2 = noprime(C1 * t1)
      @test scalar(dag_p011 * t2) ≈ +1.0
    end

    let # <011|(Cdag3*C1)|110> = -1
      t = noprime((Cdag3 * C1) * p110)
      @test scalar(dag(p011) * t) ≈ -1.0
    end

    let # <011|(C1*Cdag3)|110> = +1
      t = noprime((C1 * Cdag3) * p110)
      @test scalar(dag(p011) * t) ≈ +1.0
    end

    #
    # Commuting B tensors
    #
    # These commute by carrying additional
    # g-indices (Grassman indices)
    #

    g = Index(QN("Nf", 1, -1) => 1; tags="g")

    Bdag3 = Cdag3 * setelt(dag(g) => 1)
    B1 = setelt(g => 1) * C1

    # Commutator
    @test norm(Bdag3 * B1 - B1 * Bdag3) ≈ 0.0

    # Anti-commutator
    @test norm(Bdag3 * B1 + B1 * Bdag3) ≈ 2.0

    let # <011|Cdag3*C1|110> = <011|Bdag3*B1|110> = -1
      t1 = noprime(B1 * p110)
      t2 = noprime(Bdag3 * t1)
      @test scalar(dag(p011) * t2) ≈ -1.0
    end

    let # <011|(Cdag3*C1)|110> = <011|(Bdag3*B1)|110> = -1
      t = noprime((Bdag3 * B1) * p110)
      @test scalar(dag(p011) * t) ≈ -1.0
    end

    let # <011|Cdag3*C1|110> = <011|B1*Bdag3|110> = -1
      t1 = noprime(Bdag3 * p110)
      t2 = noprime(B1 * t1)
      @test scalar(dag(p011) * t2) ≈ -1.0
    end

    let # <011|(Cdag3*C1)|110> = <011|(B1*Bdag3)|110> = -1
      t = noprime((B1 * Bdag3) * p110)
      @test scalar(dag(p011) * t) ≈ -1.0
    end

    #
    # Leave out middle fermion, test for cases <001|...|100>
    #
    p100 = ITensor(s[1], s[2], s[3])
    p100[s[1] => 2, s[2] => 1, s[3] => 1] = 1.0

    p001 = ITensor(s[1], s[2], s[3])
    p001[s[1] => 1, s[2] => 1, s[3] => 2] = 1.0

    let # <001|Cdag3*C1|100> = <001|Bdag3*B1|100> = +1
      t1 = noprime(B1 * p100)
      t2 = noprime(Bdag3 * t1)
      @test scalar(dag(p001) * t2) ≈ +1.0
    end

    let # <001|Cdag3*C1|100> = <001|(Bdag3*B1)|100> = +1
      t = noprime((Bdag3 * B1) * p100)
      @test scalar(dag(p001) * t) ≈ +1.0
    end

    let # <001|Cdag3*C1|100> = <001|B1*Bdag3|100> = +1
      t1 = noprime(Bdag3 * p100)
      t2 = noprime(B1 * t1)
      @test scalar(dag(p001) * t2) ≈ +1.0
    end

    let # <001|Cdag3*C1|100> = <001|(B1*Bdag3)|100> = +1
      t = noprime((B1 * Bdag3) * p100)
      @test scalar(dag(p001) * t) ≈ +1.0
    end
  end

  @testset "Combiner conjugation" begin
    s = siteinds("Fermion", 4; conserve_qns=true)
    C = combiner(s[1], s[2])
    @test NDTensors.isconj(storage(C)) == false

    dC = dag(C)
    @test NDTensors.isconj(storage(dC)) == true
  end

  @testset "Combine Uncombine Permute Test" begin
    s = siteinds("Fermion", 4; conserve_qns=true)

    @testset "Two Site Test" begin
      p11 = ITensor(s[1], s[2])
      p11[s[1] => 2, s[2] => 2] = 1.0

      C = combiner(s[1], s[2])

      dp11 = dag(p11)

      Cp11_A = C * p11
      dCp11_A = dag(Cp11_A)
      dp11_A = C * dCp11_A
      @test dp11_A ≈ dp11

      Cp11_B = p11 * C
      dCp11_B = dag(Cp11_B)
      dp11_B = C * dCp11_B
      @test dp11_B ≈ dp11
    end

    @testset "Longer two-site tests" begin
      s1, s2, s3, s4 = s
      C12 = combiner(s1, s2)
      C21 = combiner(s2, s1)
      C13 = combiner(s1, s3)
      C31 = combiner(s3, s1)

      T = randomITensor(QN("Nf", 3, -1), s1, s2, s3, s4)
      T .= abs.(T)

      #
      # 1a, 2a tests
      #

      c12 = combinedind(C12)
      c12T = C12 * T
      u12T = dag(C12) * c12T
      @test norm(u12T - T) < 1E-10

      c21 = combinedind(C21)
      c21T = C21 * T
      u21T = dag(C21) * c21T
      @test norm(u21T - T) < 1E-10

      c13 = combinedind(C13)
      c13T = C13 * T
      u13T = dag(C13) * c13T
      @test norm(u13T - T) < 1E-10

      c31 = combinedind(C31)
      c31T = C31 * T
      u31T = dag(C31) * c31T
      @test norm(u31T - T) < 1E-10

      #
      # 1b, 2b tests
      #

      dc12T = dag(C12) * dag(T)
      @test norm(dc12T - dag(c12T)) < 1E-10
      du12T = C12 * dc12T
      @test norm(du12T - dag(T)) < 1E-10

      dc21T = dag(C21) * dag(T)
      @test norm(dc21T - dag(c21T)) < 1E-10
      du21T = C21 * dc21T
      @test norm(du21T - dag(T)) < 1E-10

      dc13T = dag(C13) * dag(T)
      @test norm(dc13T - dag(c13T)) < 1E-10
      du13T = C13 * dc13T
      @test norm(du13T - dag(T)) < 1E-10

      dc31T = dag(C31) * dag(T)
      @test norm(dc31T - dag(c31T)) < 1E-10
      du31T = C31 * dc31T
      @test norm(du31T - dag(T)) < 1E-10
    end

    @testset "Three Site Test" begin
      p111 = ITensor(s[1], s[2], s[3])
      p111[s[1] => 2, s[2] => 2, s[3] => 2] = 1.0

      dp111 = dag(p111)

      C = combiner(s[1], s[3])
      Cp111 = C * p111
      dCp111 = dag(Cp111)
      dp111_U = C * dCp111
      @test dp111_U ≈ dp111
    end
  end

  @testset "Mixed Arrow Combiner Tests" begin
    @testset "One wrong-way arrow" begin
      q1 = QN("Nf", 1, -1)

      s0 = Index([q1 => 1]; tags="s0")
      s1 = Index([q1 => 1]; tags="s1")
      s2 = Index([q1 => 1]; tags="s2")
      s3 = Index([q1 => 1]; tags="s3")
      s4 = Index([q1 => 1]; tags="s4")

      A = randomITensor(QN("Nf", 0, -1), s0, s1, dag(s2), dag(s3))
      B = randomITensor(QN("Nf", 0, -1), s3, s2, dag(s1), dag(s4))
      A .= one.(A)
      B .= one.(B)
      @test norm(A) ≈ 1.0
      @test norm(B) ≈ 1.0

      Ru = A * B

      C = combiner(s3, s2, dag(s1))
      Bc = C * B
      Ac = A * dag(C)
      Rc = Ac * Bc

      @test norm(Ru - Rc) < 1E-8
    end

    @testset "Two wrong-way arrows" begin
      q1 = QN("Nf", 1, -1)

      s0 = Index([q1 => 1]; tags="s0")
      s1 = Index([q1 => 1]; tags="s1")
      s2 = Index([q1 => 1]; tags="s2")
      s3 = Index([q1 => 1]; tags="s3")
      s4 = Index([q1 => 1]; tags="s4")

      A = randomITensor(QN("Nf", 2, -1), s0, s1, s2, dag(s3))
      B = randomITensor(QN("Nf", -2, -1), s3, dag(s2), dag(s1), dag(s4))
      A .= one.(A)
      B .= one.(B)
      @test norm(A) ≈ 1.0
      @test norm(B) ≈ 1.0

      Ru = A * B

      C = combiner(s3, dag(s2), dag(s1))
      Bc = C * B
      Ac = A * dag(C)
      Rc = Ac * Bc

      @test norm(Ru - Rc) < 1E-8
    end
  end

  @testset "Permutedims Regression Test" begin
    s1 = Index([QN("N", 0, -1) => 1, QN("N", 1, -1) => 1], "s1")
    s2 = Index([QN("N", 0, -1) => 1, QN("N", 1, -1) => 1], "s2")
    i = Index([QN("N", 0, -1) => 1, QN("N", 1, -1) => 1, QN("N", 2, -1) => 1], "i")

    A = ITensor(QN("N", 4, -1), s1, s2, i)
    A[s1 => 2, s2 => 2, i => 3] = 223

    B = ITensor(QN("N", 4, -1), s1, i, s2)
    B[s1 => 2, i => 3, s2 => 2] = 223
    @test A ≈ B

    C = ITensor(QN("N", 4, -1), s1, i, s2)
    C[s2 => 2, i => 3, s1 => 2] = -223
    @test A ≈ C
  end

  @testset "Product MPS consistency checks" begin
    s = siteinds("Fermion", 3; conserve_qns=true)

    pA = productMPS(s, [2, 1, 2])
    TA = ITensor(s[1], s[2], s[3])
    TA[s[1] => 2, s[2] => 1, s[3] => 2] = 1.0
    A = pA[1] * pA[2] * pA[3]
    @test norm(A - TA) < 1E-8

    pB = productMPS(s, [1, 2, 2])
    TB = ITensor(s[1], s[2], s[3])
    TB[s[1] => 1, s[2] => 2, s[3] => 2] = 1.0
    B = pB[1] * pB[2] * pB[3]
    @test norm(B - TB) < 1E-8
  end

  @testset "Orthogonalize of Product MPS" begin
    N = 3

    sites = siteinds("Fermion", N; conserve_qns=true)

    state = [1 for n in 1:N]
    state[1] = 2
    state[2] = 2
    psi = productMPS(sites, state)
    psi_fluxes = [flux(psi[n]) for n in 1:N]

    psi_orig = copy(psi)
    orthogonalize!(psi, 1)
    @test inner(psi_orig, psi) ≈ 1.0
    @test inner(psi, psi_orig) ≈ 1.0
  end

  @testset "MPS inner regression test" begin
    sites = siteinds("Fermion", 3; conserve_qns=true)
    psi = productMPS(sites, [2, 2, 1])
    @test inner(psi, psi) ≈ 1.0
  end

  @testset "Fermionic AutoMPO Tests" begin
    @testset "Spinless Fermion Hamiltonian" begin
      N = 2
      sites = siteinds("Fermion", N; conserve_qns=true)
      t1 = 1.0
      ampo = AutoMPO()
      for b in 1:(N - 1)
        ampo += -t1, "Cdag", b, "C", b + 1
        ampo += -t1, "Cdag", b + 1, "C", b
      end
      H = MPO(ampo, sites)

      HH = H[1]
      for n in 2:N
        HH *= H[n]
      end
      HHc = dag(swapprime(HH, 0, 1))
      @test norm(HHc - HH) < 1E-8
    end

    @testset "Fermion Hamiltonian Matrix Elements" begin
      N = 10
      t1 = 0.654
      V1 = 1.23

      sites = siteinds("Fermion", N; conserve_qns=true)

      ampo = AutoMPO()
      for b in 1:(N - 1)
        ampo += -t1, "Cdag", b, "C", b + 1
        ampo += -t1, "Cdag", b + 1, "C", b
        ampo += V1, "N", b, "N", b + 1
      end
      H = MPO(ampo, sites)

      for j in 1:(N - 2)
        stateA = [1 for n in 1:N]
        stateA[j] = 2
        stateA[N] = 2 # to make MPS bosonic

        stateB = [1 for n in 1:N]
        stateB[j + 1] = 2
        stateB[N] = 2 # to make MPS bosonic

        psiA = productMPS(sites, stateA)
        psiB = productMPS(sites, stateB)

        @test inner(psiA', H, psiB) ≈ -t1
        @test inner(psiB', H, psiA) ≈ -t1
      end

      for j in 1:(N - 1)
        state = [1 for n in 1:N]
        state[j] = 2
        state[j + 1] = 2
        psi = productMPS(sites, state)
        @test inner(psi', H, psi) ≈ V1
      end
    end

    @testset "Fermion Second Neighbor Hopping" begin
      N = 4
      t1 = 1.79
      t2 = 0.427
      s = siteinds("Fermion", N; conserve_qns=true)
      a = AutoMPO()
      for n in 1:(N - 1)
        a += -t1, "Cdag", n, "C", n + 1
        a += -t1, "Cdag", n + 1, "C", n
      end
      for n in 1:(N - 2)
        a += -t2, "Cdag", n, "C", n + 2
        a += -t2, "Cdag", n + 2, "C", n
      end
      H = MPO(a, s)

      state1 = [1 for n in 1:N]
      state1[1] = 2
      state1[4] = 2
      psi1 = productMPS(s, state1)

      state2 = [1 for n in 1:N]
      state2[2] = 2
      state2[4] = 2
      psi2 = productMPS(s, state2)

      state3 = [1 for n in 1:N]
      state3[3] = 2
      state3[4] = 2
      psi3 = productMPS(s, state3)

      @test inner(psi1', H, psi2) ≈ -t1
      @test inner(psi2', H, psi1) ≈ -t1
      @test inner(psi2', H, psi3) ≈ -t1
      @test inner(psi3', H, psi2) ≈ -t1

      @test inner(psi1', H, psi3) ≈ -t2
      @test inner(psi3', H, psi1) ≈ -t2

      # Add stationary particle to site 2,
      # hopping over should change sign:
      state1[2] = 2
      psi1 = productMPS(s, state1)
      state3[2] = 2
      psi3 = productMPS(s, state3)
      @test inner(psi1', H, psi3) ≈ +t2
      @test inner(psi3', H, psi1) ≈ +t2
    end
  end

  @testset "Fermionic SVD" begin
    N = 4
    s = siteinds("Fermion", N; conserve_qns=true)

    A = randomITensor(QN("Nf", 2, -1), s[1], s[2], s[3], s[4])
    for n1 in 1:4, n2 in 1:4
      (n1 == n2) && continue
      U, S, V = svd(A, (s[n1], s[n2]))
      @test norm(U * S * V - A) < 1E-10
    end
    for n1 in 1:4, n2 in 1:4, n3 in 1:4
      (n1 == n2) && continue
      (n1 == n3) && continue
      (n2 == n3) && continue
      U, S, V = svd(A, (s[n1], s[n2], s[n3]))
      @test norm(U * S * V - A) < 1E-10
    end

    B = randomITensor(QN("Nf", 3, -1), s[1], s[2], s[3], s[4])
    for n1 in 1:4, n2 in 1:4
      (n1 == n2) && continue
      U, S, V = svd(B, (s[n1], s[n2]))
      @test norm(U * S * V - B) < 1E-10
    end
    for n1 in 1:4, n2 in 1:4, n3 in 1:4
      (n1 == n2) && continue
      (n1 == n3) && continue
      (n2 == n3) && continue
      U, S, V = svd(B, (s[n1], s[n2], s[n3]))
      @test norm(U * S * V - B) < 1E-10
    end
  end # Fermionic SVD tests

  @testset "Fermion Contraction with Combined Indices" begin
    N = 10
    s = siteinds("Fermion", N; conserve_qns=true)

    begin
      A = randomITensor(QN("Nf", 3, -1), s[1], s[2], s[3], s[4])
      B = randomITensor(QN("Nf", 2, -1), s[1], s[3], s[4])

      CC = combiner(s[1], s[3])

      cA = CC * A
      cB = CC * B

      R1 = dag(cA) * cB
      R2 = dag(A) * B

      @test norm(R1 - R2) < 1E-10
    end

    begin
      A = randomITensor(QN("Nf", 3, -1), s[1], s[2], s[3], s[4])
      B = randomITensor(QN("Nf", 2, -1), s[1], s[3], s[4])

      CC = combiner(s[1], s[3])

      cA = CC * A
      cdB = dag(CC) * dag(B)

      R1 = cA * cdB
      R2 = A * dag(B)

      @test norm(R1 - R2) < 1E-10
    end

    begin
      A = randomITensor(QN("Nf", 3, -1), s[1], s[2], s[3], s[4])
      B = randomITensor(QN("Nf", 2, -1), s[1], s[3], s[4])

      CC = combiner(s[1], s[4], s[3])

      cA = CC * A
      cdB = dag(CC) * dag(B)

      R1 = cA * cdB
      R2 = A * dag(B)

      @test norm(R1 - R2) < 1E-10
    end

    begin
      CC = combiner(s[3], s[4])
      c = combinedind(CC)

      A = randomITensor(QN("Nf", 3, -1), c, s[1], s[2])
      B = randomITensor(QN("Nf", 2, -1), s[1], c, s[5])

      uA = dag(CC) * A
      uB = dag(CC) * B

      R1 = dag(uA) * uB
      R2 = dag(A) * B

      @test norm(R1 - R2) < 1E-10
    end

    @testset "Combiner Regression Test" begin
      T = randomITensor(QN("Nf", 3, -1), s[1], s[2], s[3], s[4])

      C12 = combiner(s[1], s[2])
      c12 = combinedind(C12)
      c12T = C12 * T

      u12T = dag(C12) * c12T

      @test norm(u12T - T) < 1E-10
    end
  end # Fermion Contraction with Combined Indices

  @testset "DMRG Tests" begin
    @testset "Nearest Neighbor Fermions" begin
      N = 8
      t1 = 1.0
      V1 = 4.0

      s = siteinds("Fermion", N; conserve_qns=true)

      ampot = AutoMPO()
      ampoV = AutoMPO()
      for b in 1:(N - 1)
        ampot += -t1, "Cdag", b, "C", b + 1
        ampot += -t1, "Cdag", b + 1, "C", b
        ampoV += V1, "N", b, "N", b + 1
      end
      Ht = MPO(ampot, s)
      HV = MPO(ampoV, s)

      state = ["Emp" for n in 1:N]
      for i in 1:2:N
        state[i] = "Occ"
      end
      psi0 = productMPS(s, state)

      sweeps = Sweeps(3)
      maxdim!(sweeps, 20, 20, 40, 80, 200)
      cutoff!(sweeps, 1E-6)

      correct_energy = -2.859778

      energy, psi = dmrg([Ht, HV], psi0, sweeps; outputlevel=0)
      @test abs(energy - correct_energy) < 1E-4

      # Test using SVD within DMRG too:
      energy, psi = dmrg([Ht, HV], psi0, sweeps; outputlevel=0, which_decomp="svd")
      @test abs(energy - correct_energy) < 1E-4
    end

    @testset "Further Neighbor and Correlations" begin
      N = 8
      t1 = 1.0
      t2 = 0.2

      s = siteinds("Fermion", N; conserve_qns=true)

      ampot = AutoMPO()
      for b in 1:(N - 1)
        ampot += -t1, "Cdag", b, "C", b + 1
        ampot += -t1, "Cdag", b + 1, "C", b
      end
      for b in 1:(N - 2)
        ampot += -t2, "Cdag", b, "C", b + 2
        ampot += -t2, "Cdag", b + 2, "C", b
      end
      Ht = MPO(ampot, s)

      state = ["Emp" for n in 1:N]
      for i in 1:2:N
        state[i] = "Occ"
      end
      psi0 = productMPS(s, state)

      sweeps = Sweeps(3)
      maxdim!(sweeps, 20, 20, 40, 80, 200)
      cutoff!(sweeps, 1E-6)

      energy, psi = dmrg(Ht, psi0, sweeps; outputlevel=0)

      energy_inner = inner(psi', Ht, psi)

      C = correlation_matrix(psi, "Cdag", "C")
      C_energy =
        sum(j -> -2t1 * C[j, j + 1], 1:(N - 1)) + sum(j -> -2t2 * C[j, j + 2], 1:(N - 2))

      @test energy_inner ≈ energy
      @test C_energy ≈ energy
    end
  end

  @testset "MPS gate system" begin
    @testset "Fermion sites" begin
      N = 3

      s = siteinds("Fermion", N; conserve_qns=true)

      # Ground state |000⟩
      ψ000 = productMPS(s, "0")

      # Start state |011⟩
      ψ011 = productMPS(s, n -> n == 2 || n == 3 ? "1" : "0")

      # Reference state |110⟩
      ψ110 = productMPS(s, n -> n == 1 || n == 2 ? "1" : "0")

      function ITensors.op(::OpName"CdagC", ::SiteType, s1::Index, s2::Index)
        return op("Cdag", s1) * op("C", s2)
      end

      os = [("CdagC", 1, 3)]
      Os = ops(os, s)

      # Results in -|110⟩
      ψ1 = product(Os, ψ011; cutoff=1e-15)

      @test inner(ψ1, ψ110) == -1

      a = AutoMPO()
      a += "Cdag", 1, "C", 3
      H = MPO(a, s)

      # Results in -|110⟩
      ψ2 = noprime(contract(H, ψ011; cutoff=1e-15))

      @test inner(ψ2, ψ110) == -1
    end
  end

  @testset "Regression Tests" begin
    @testset "SVD DiagBlockSparse Regression Test" begin
      l1 = Index(QN("Nf", 0, -1) => 1, QN("Nf", 1, -1) => 1; tags="Link,l=1")
      s2 = Index(QN("Nf", 0, -1) => 1, QN("Nf", 1, -1) => 1; tags="Site,n=2")
      s3 = Index(QN("Nf", 0, -1) => 1, QN("Nf", 1, -1) => 1; tags="Site,n=3")
      l3 = Index(QN("Nf", 2, -1) => 1; tags="Link,l=3")

      phi = randomITensor(QN("Nf", 4, -1), l1, s2, s3, l3)

      U, S, V = svd(phi, (l1, s2))

      @test norm((U * S) * V - phi) < 1E-10
      @test norm(U * (S * V) - phi) < 1E-10
    end

    @testset "Eigen Positive Semi Def Regression Test" begin
      #
      # Test was failing without using combiners in
      # eigen which were conjugates of each other
      #
      cutoff = 1E-12
      N = 2
      s = siteinds("Fermion", N; conserve_qns=true)

      T = ITensor(QN("Nf", 0, -1), dag(s[1]), s[1]')
      T[2, 2] = 1

      F = eigen(T; ishermitian=true, cutoff=cutoff)
      D, U, spec = F
      Ut = F.Vt

      @test norm(dag(U) * D * Ut - T) < 1E-10
    end

    @testset "Factorize Eigen Regression Test" begin
      N = 3
      s = siteinds("Fermion", N; conserve_qns=true)
      A = ITensor(QN("Nf", 2, -1), s[1], s[2], s[3])
      A[s[1] => 1, s[2] => 2, s[3] => 2] = 1.0

      U, R = factorize(A, (s[1], s[2]); which_decomp="eigen", cutoff=1E-18, ortho="left")

      @test norm(U * R - A) < 1E-12
    end

    @testset "Contraction Regression Test" begin
      s = siteinds("Fermion", 3; conserve_qns=true)
      l = Index(QN("Nf", 1, -1) => 1; tags="l")

      q2 = QN("Nf", 2, -1)
      q0 = QN("Nf", 0, -1)

      T1 = ITensor(q2, s[1], s[2], l)
      T1[s[1] => 1, s[2] => 2, l => 1] = 1.0

      T2 = ITensor(q0, dag(l), s[3])
      T2[dag(l) => 1, s[3] => 2] = 1.0

      @test norm(T1 * T2 - T2 * T1) < 1E-10
    end

    @testset "SVD Regression Test" begin
      Pf0 = QN("Pf", 0, -2)
      Pf1 = QN("Pf", 1, -2)

      l22 = Index([Pf0 => 1, Pf1 => 1], "Link,dir=2,n=2")
      l23 = Index([Pf0 => 1, Pf1 => 1], "Link,dir=3,n=2")
      s1 = Index([Pf0 => 1, Pf1 => 1, Pf1 => 1, Pf0 => 1], "Site,n=1")
      l11 = Index([Pf0 => 1, Pf1 => 1], "Link,dir=1,n=1")

      T = randomITensor(dag(l22), dag(l23), s1, l11)

      U, S, V = svd(T, dag(l22), dag(l23), s1)

      @test norm(T - U * S * V) < 1E-10
    end

    @testset "AutoMPO Regression Test" begin
      N = 3
      s = siteinds("Fermion", N; conserve_qns=true)

      a = AutoMPO()
      a += "Cdag", 1, "C", 3
      @test_nowarn H = MPO(a, s)
    end
  end # Regression Tests

  ITensors.disable_auto_fermion()
end

nothing


text:
using CompScienceMeshes
using BEAST
using Test

m = meshrectangle(1.0, 1.0, 0.25)

X = raviartthomas(m, BEAST.Continuity{:none})
@test numfunctions(X) == 16*2*3
@test all(length.(X.fns) .== 1)

p = positions(X)

i = 12
c = X.fns[i][1].cellid
ch = chart(m, cells(m)[c])
ctr = cartesian(center(ch))
@test ctr ≈ p[i]

for (i,f) in enumerate(X.fns)
    @test f[1].coeff == 1.0
end

ext: jl


text:
module Radon
export radon!
# package code goes here

function radon{T<:AbstractFloat}(im::AbstractArray{T, 2})
    Nx, Ny = size(im)
    radonim = similar(im, min(Nx, Ny), max(Nx, Ny))
    radon!(radonim, im)
    radonim
end

function radon!{T<:AbstractFloat}(radonim::AbstractArray{T, 2}, im::AbstractArray{T, 2})
    Nφ, Ns = size(radonim)
    Nx, Ny = size(im)
    Ldiag = hypot(Nx, Ny)
    K = round(Int, Ldiag)

    dLz = Ldiag/(K-1)
    dLs = Ldiag/(Ns-1)
    dφ = π/Nφ
    for is = 1:Ns, iφ = 1:Nφ
        φ = (iφ-1)*dφ
        s = -Ldiag/2 + (is-1)*dLs
        tx = cos(φ)
        ty = sin(φ)

        radonim[iφ, is] = zero(T)
        for k=1:K
            z = -Ldiag/2 + (k-1)*dLz
            x = round(Int, Nx/2 + z*ty + s*tx)
            y = round(Int, Ny/2 - z*tx + s*ty)

            if 1 <= x <= Nx && 1 <= y <= Ny
                radonim[iφ, is] += im[x, y]*dLz
            end
        end
    end
    radonim
end

function backproject{T<:AbstractFloat}(radonim::AbstractArray{T, 2})
    Nφ, Ns = size(radonim)
    im = similar(radonim, Ns, Ns)
    backproject!(im, radonim)
    im
end

function backproject!{T<:AbstractFloat}(im::AbstractArray{T, 2}, radonim::AbstractArray{T, 2})
    Nφ, Ns = size(radonim)
    Nx, Ny = size(im)
    Ldiag = hypot(Nx, Ny)
    K = round(Int, Ldiag)

    dLz = Ldiag/(K-1)
    dLs = Ldiag/(Ns-1)
    dφ = π/Nφ

    im[:] = 0.
    for is = 1:Ns, iφ = 1:Nφ
        φ = (iφ-1)*dφ
        s = -Ldiag/2 + (is-1)*dLs
        tx = cos(φ)
        ty = sin(φ)

        for k=1:K
            z = -Ldiag/2 + (k-1)*dLz
            x = round(Int, Nx/2 + z*ty + s*tx)
            y = round(Int, Ny/2 - z*tx + s*ty)

            if 1 <= x <= Nx && 1 <= y <= Ny
                im[x, y] += radonim[iφ, is]*dLz
            end
        end
    end
    im
end

end # module


text:
using FountainCodes, Test, SparseArrays

function init(K)
    error("Deprecated")
    p = FountainCodes.RQ(K)
    d = FountainCodes.Decoder(p)
    C = [Vector{GF256}([i % 256]) for i in 1:p.L]
    precode!(C, p)
    return p, d, C
end

# test RQ parameter choice
@test FountainCodes.RQ_parameters(27) == (30, 566, 11, 10, 41)
@test FountainCodes.RQ_parameters(30) == (30, 566, 11, 10, 41)
@test FountainCodes.RQ_parameters(90) == (91, 66, 17, 10, 103)

# test degree distribution
@test FountainCodes.RQ_deg(0) == 1
@test FountainCodes.RQ_deg(5243) == 2
@test FountainCodes.RQ_deg(1048576-1) == 30

"""test RaptorQ tuple function"""
function test_RQ_tuple()
    c = RQ(10)
    for X in 1:10
        d, a, b, d1, a1, b1 = FountainCodes.RQ_tuple(X, c)
        if !(0 < d)
            error("d must be positive")
        end
        if !(1 <= a <= c.W-1)
            error("a must be between 1 and W-1 inclusive")
        end
        if !(0 <= b <= c.W-1)
            error("b must be between 0 and W-1 inclusive")
        end
        if !(d1 in [2, 3])
            error("d1 must be either 2 or 3")
        end
        if !(1 <= a1 <= c.P1-1)
            error("a1 must be between 1 and P1-1 inclusive")
        end
        if !(0 <= b1 <= c.P1-1)
            error("b1 must be between 0 and P1-1 inclusive")
        end
    end
    return true
end
@test test_RQ_tuple()

"""test LDPC constraints against known correct constraints"""
function test_ldpc_constraints()
    c = RQ(10)
    A = FountainCodes.ldpc_constraint_matrix(c)

    ri = 1
    Is = [1, 6, 7, 8, 11, 18, 19]
    Vs = ones(GF256, 7)
    correct = sparsevec(Is, Vs, 27)
    ans = A[:, ri]
    @test ans == correct

    ri = 2
    Is = [1, 2, 7, 9, 12, 19, 20]
    Vs = ones(GF256, 7)
    correct = sparsevec(Is, Vs, 27)
    ans = A[:, ri]
    @test ans == correct

    ri = 3
    Is = [1, 2, 3, 8, 10, 13, 20, 21]
    Vs = ones(GF256, 8)
    correct = sparsevec(Is, Vs, 27)
    ans = A[:, ri]
    @test ans == correct

    ri = 7
    Is = [5, 6, 7, 10, 17, 24, 25]
    Vs = ones(GF256, 7)
    correct = sparsevec(Is, Vs, 27)
    ans = A[:, ri]
    @test ans == correct
    return
end
test_ldpc_constraints()

"""test HDPC constraints against known correct constraints"""
function test_hdpc_constraints()
    c = RQ(10)
    A = FountainCodes.hdpc_constraint_matrix(c)

    i = 1
    Is = append!(collect(1:17), 18)
    Vs = GF256.([
        0xfa, 0xf3, 0xf7, 0xf5, 0xf4, 0xf4, 0xf4, 0x7a, 0x3d,
        0x90, 0x48, 0x24, 0x12, 0x09, 0x04, 0x02, 0x01, 0x01,
    ])
    correct = sparsevec(Is, Vs, 27)
    ans = A[:, i]
    @test ans == correct

    i = 10
    Is = append!(collect(1:17), 27)
    Vs = GF256.([
        0xeb, 0x75, 0x3a, 0x1d, 0x80, 0x40, 0x20, 0x10, 0x08,
        0x04, 0x02, 0x01, 0x8e, 0xc9, 0xea, 0x75, 0x3a ,0x01,
    ])
    correct = sparsevec(Is, Vs, 27)
    ans = A[:, i]
    @test ans == correct
    return
end
test_hdpc_constraints()

"""Test that decoding with ESIs 0 through Kp-1 succeeds"""
function test_precode(K::Integer)
    rng = MersenneTwister(123)
    code = RQ(K)
    K = code.Kp # padding to the closest supported number of source symbols
    src = rand_nonzero(rng, GF256, K)

    # compute the intermediate symbols
    A = constraint_matrix(code, 0:K-1)
    pre = zeros(GF256, size(A, 2))
    pre[end-K+1:end] .= src
    decoder = Decoder(A)

    # permanent inactivations
    for cpi in (code.L-code.P+1):code.L
        FountainCodes.mark_inactive!(decoder, A, cpi)
    end

    # decode
    int = decode(A, pre; decoder)

    # test that the first K coded symbols are equal to the src symbols,
    # i.e., that the code is systematic
    G = generator_matrix(code, 0:K-1)
    enc = G'*int
    @test enc[1:K] == src
end
for K in 4:100
    test_precode(K)
end
test_precode(7000)
test_precode(8192)

function test_decoding(K::Integer)
    rng = MersenneTwister(123)
    code = RQ(K)
    K = code.Kp # padding to the closest supported number of source symbols
    src = rand_nonzero(rng, GF256, K)
    Xs = code.Kp:(2*code.Kp-1)

    # compute the intermediate symbols
    A = constraint_matrix(code, 0:K-1)
    pre = zeros(GF256, size(A, 2))
    pre[end-K+1:end] .= src

    ## permanent inactivations
    decoder = Decoder(A)
    for cpi in (code.L-code.P+1):code.L
        FountainCodes.mark_inactive!(decoder, A, cpi)
    end
    int = decode(A, pre; decoder)

    # compute the encoded symbols
    G = generator_matrix(code, Xs)
    enc = G'*int

    # decode the intermediate symbols from the encoded symbols
    A = constraint_matrix(code, Xs)
    pre = zeros(GF256, size(A, 2))
    pre[end-length(enc)+1:end] .= enc

    ## permanent inactivations
    decoder = Decoder(A)
    for cpi in (code.L-code.P+1):code.L
        FountainCodes.mark_inactive!(decoder, A, cpi)
    end
    dec_int = decode(A, pre)

    # generate the source symbols from the decoded intermediate symbols
    G = generator_matrix(code, 0:K-1)
    dec_src = G'*dec_int
    @test dec_src == src
end
for K in 4:100
    test_decoding(K)
end

max_stars_repo_name: severinson/RaptorCodes

max_issues_repo_name: severinson/RaptorCodes

max_forks_repo_name: severinson/RaptorCodes


text:
# Create example data
x = LinRange(0, 30, 1000)
y = cos.(x) .* x
z = sin.(x) .* x
# Draw a solid line and another with star markers
# in one of every 10 points
plot3(x, y, z, x[1:10:end], y[1:10:end], z[1:10:end], "p")


text:
using PQBaseCamp
using DataFrames
using CSV

# compute array -
compute_model_array = Array{LinearReturnComputationModel,1}()

# ticker_array -
ticker_array = [
    "MRNA", "PFE", "JNJ", "SPY", "AAPL"
]

# build up a list of compute objects -
for ticker ∈ ticker_array

    # load up a sample data file -
    path_to_data_file = "$(pwd())/test/data/$(ticker).csv"
    df = CSV.read(path_to_data_file, DataFrame)

    # compute -
    model = LinearReturnComputationModel()
    model.ticker = ticker
    model.data = df
    model.map = :timestamp => :close

    # grab -
    push!(compute_model_array, model)
end

# compute the cov array -
# first: compute the dictionary of returns -
price_retrun_dictionary = Δ(compute_model_array; multiplier = 100.0)
covm = covariance(ticker_array, price_retrun_dictionary)

# compute the β_array -
β_array = β(ticker_array, price_retrun_dictionary)


text:
struct TraceSim <: MethodType
    α::Float64
    β::Float64
    γ::Float64

    sigmoid::Bool
    b::Float64

    idf::Bool
    reciprocal_func::Bool
    no_norm::Bool

    sum_match::Bool

    local_weights::Vector{Float64}
end

function TraceSim(α::Float64, β::Float64, γ::Float64, sigmoid::Bool, b::Float64, idf::Bool, reciprocal_func::Bool,  no_norm::Bool, sum_match::Bool=false)
    local_weights = Vector{Float64}(undef, 400)

    @simd for pos in 1:400
        @inbounds local_weights[pos] = reciprocal_func ? 1.0 / ((pos)^α) : exp(-α * (pos - 1.0))
    end

    @info "TraceSim: sum_match=$(sum_match)"
    return TraceSim(α, β, γ, sigmoid, b, idf, reciprocal_func, no_norm, sum_match, local_weights)
end

computeDocFreq(method::TraceSim, docFreq) = @. max(exp(-method.β * (docFreq.doc_freq / docFreq.nDocs) * 100.0), 1.0e-100)

@inline function frame_weight(method::TraceSim, pos, df)
    @inbounds lw = method.local_weights[pos]

    return lw * df
end

gap_score(method::TraceSim, pos, df) = -1.0 * frame_weight(method, pos, df)
@inline function match_score(method::TraceSim, pos1, pos2, df, qw, cw)
    diff = exp(-method.γ * abs(pos1 - pos2))

    if method.sum_match
        return (qw + cw) * diff
    else
        return max(qw, cw) * diff
    end
end

function compute_frame_weights(method::TraceSim, stacktrace::Vector{UInt32}, doc_freq::Vector{Float64})
    weights = Vector{Float64}(undef, length(stacktrace))

    @simd for pos in eachindex(stacktrace)
        @inbounds weights[pos] = frame_weight(method, pos, doc_freq[stacktrace[pos]])
    end

    return weights
end

@inline norm_sum(sim, q_weights::Vector{Float64}, c_weights::Vector{Float64})= sim/(sum(q_weights) + sum(c_weights))

function norm_jaccard(sim, query::Vector{UInt32}, candidate::Vector{UInt32}, q_weights::Vector{Float64}, c_weights::Vector{Float64},q_sorted_idxs, c_sorted_idxs)
    i = 1
    j = 1
    den = 0.0

    query_len = length(query) + 1
    cand_len = length(candidate) + 1

    func_id = 0x00000

    @inbounds while i < query_len || j < cand_len
        aux1 = 0.0

        if i == query_len || (j < cand_len && query[q_sorted_idxs[i]] > candidate[c_sorted_idxs[j]])
            # Query does not contain the candidate[c_sorted_idxs[j]]
            func_id = candidate[c_sorted_idxs[j]]

            while j < cand_len
                candidate[c_sorted_idxs[j]] != func_id && break

                aux1 += c_weights[c_sorted_idxs[j]]
                j += 1
            end

            den += aux1
        elseif j == cand_len || ( i < query_len && query[q_sorted_idxs[i]] < candidate[c_sorted_idxs[j]])
            # Candidate does not contain the query[q_sorted_idxs[i]]
            func_id = query[q_sorted_idxs[i]]

            while i < query_len
                query[q_sorted_idxs[i]] != func_id && break

                aux1 += q_weights[q_sorted_idxs[i]]
                i += 1
            end

            den += aux1
        else
            # Both sequences contain a function x. Get maximum weights between them

            # Sum of the weights of a function in candidate
            func_id = candidate[ c_sorted_idxs[j]]

            while j < cand_len
                candidate[ c_sorted_idxs[j]] != func_id && break
                aux1 += c_weights[ c_sorted_idxs[j]]
                j += 1
            end

            # Sum of the weights of a function in query
            aux2 = 0.0
            func_id = query[q_sorted_idxs[i]]

            while i < query_len
                query[q_sorted_idxs[i]] != func_id && break

                aux2 += q_weights[q_sorted_idxs[i]]
                i += 1
            end

            den += max(aux1, aux2)
        end
    end

    den == 0.0 && return -1.0

    return sim / den
end

@inline function similarity(method::TraceSim, query::Vector{UInt32}, candidate::Vector{UInt32}, doc_freq::Vector{Float64})
    q_weights = compute_frame_weights(method, query, doc_freq)

    return similarity(method, query, candidate, q_weights, sortperm(query), sortperm(candidate), doc_freq)
end

function similarity(method::TraceSim, query::Vector{UInt32}, candidate::Vector{UInt32}, q_weights::Vector{Float64}, q_sorted_idxs, c_sorted_idxs, doc_freq::Vector{Float64})
    c_weights = compute_frame_weights(method, candidate, doc_freq)

    query_len = length(q_weights)
    cand_len = length(c_weights)
    offset = cand_len + 1

    # previous_row = zeros(offset)
    row = zeros(offset)

    # Create initial values
    for i in 1:cand_len
        @inbounds row[i + 1] =  row[i] - c_weights[i]
    end

    @inbounds for (q_pos, q_func) in enumerate(query)
        # Set first column of the row
        diagonal = row[1]
        row[1] -= q_weights[q_pos]

        for (c_pos, c_func) in enumerate(candidate)
            col_idx = c_pos + 1

            # Align gap to query position
            above = row[col_idx] - q_weights[q_pos]
            # println("\tPrevious row= $(previous_row[col_idx]) - $(q_weights[q_pos])")

            # Align gap to candidate position
            left = row[col_idx - 1] - c_weights[c_pos]
            # println("\tPrevious col= $(cur_row[col_idx - 1]) - $(c_weights[c_pos])")

            # Value of diagonal
            # diagonal = previous_row[col_idx - 1]

            if q_func == c_func
                # Matching
                diagonal += match_score(method, q_pos, c_pos, 0.0, q_weights[q_pos], c_weights[c_pos])
            else
                # Mismatch
                diagonal -= q_weights[q_pos] + c_weights[c_pos]
                # println("\tDiagonal = $(previous_row[col_idx - 1]) - $( q_weights[q_pos] + c_weights[c_pos])")
            end

            diagonal_tmp = row[col_idx]

            # Replace max with IF ELSE considerably improves the performance =O
            # row[col_idx] = max(left, above, diagonal)
            m = left

            if above > m
                m = above
            end

            if diagonal > m
                m = diagonal
            end

            row[col_idx] = m

            diagonal = diagonal_tmp
            # println("$(q_pos), $(c_pos) max($(left), $(diagonal_value), $(above))")
        end
    end

    sim = row[offset]

    if !method.no_norm
        if method.sum_match
            return norm_sum(sim, q_weights, c_weights)
        else
            return norm_jaccard(sim, query, candidate, q_weights, c_weights, q_sorted_idxs, c_sorted_idxs)
        end
    end

    return sim
end

# @inline function compute_final_similarity!(partial_similarities::Vector{Float64}, method::TraceSim, query::Report, candidates::Vector{Report}, df_vec::Vector{Float64})
#     query_weight_sum = 0.0
#     q_weights = Vector{Float64}(undef, query.stats.doc_len)
#     q_frames = Vector{UInt32}(undef, query.stats.doc_len)

#     stack_idx = 1

#     for stack in query.stacks
#         for (pos, func_id) in enumerate(stack)
#             fw = frame_weight(method, pos, df_vec[func_id])

#             q_weights[stack_idx] = fw
#             q_frames[stack_idx] = func_id

#             query_weight_sum -= fw
#             stack_idx+=1
#         end
#     end

#     Threads.@threads for idx in 1:length(partial_similarities)
#         sim = partial_similarities[idx]

#         if sim == 0.0
#             partial_similarities[idx] = -1.0
#         else
#             weight_sum = query_weight_sum
#             @inbounds candidate = candidates[idx]

#             c_weights = Vector{Float64}(undef, candidate.stats.doc_len)
#             c_frames = Vector{UInt32}(undef, candidate.stats.doc_len)

#             c_stack_idx = 1

#             for stack in candidate.stacks
#                 for (pos, func_id) in enumerate(stack)
#                     @inbounds fw = frame_weight(method, pos, df_vec[func_id])
#                     c_weights[c_stack_idx] = fw
#                     c_frames[c_stack_idx] = func_id

#                     weight_sum -= fw
#                     c_stack_idx+=1
#                 end
#             end

#             sim += weight_sum

#             partial_similarities[idx] = normalize_sim(sim, q_frames, c_frames, q_weights, c_weights, sortperm(q_frames), sortperm(c_frames))
#         end
#     end

#     return partial_similarities
# end

# function similarity(method::TraceSim, query::Vector{UInt32}, candidate::Vector{UInt32}, doc_freq::Vector{Float64})
#     q_weights = compute_frame_weights(method, query, doc_freq)
#     c_weights = compute_frame_weights(method, candidate,doc_freq)

#     offset = length(candidate) + 1
#     M = zeros(2 * offset)

#     query_len = length(q_weights)
#     cand_len = length(c_weights)

#     ONE_AGO = 0
#     THIS_ROW = 1

#     # Create initial values
#     for i in 1:cand_len
#         M[THIS_ROW * offset + (i + 1)] =  M[THIS_ROW * offset + i] - c_weights[i]
#     end

#     # println(M[THIS_ROW * offset + 1 : THIS_ROW * offset + offset])
#     for (q_pos, q_func) in enumerate(query)
#         # # Copy THIS_ROW to ONE_AGO
#         for k in 1:offset
#             M[ONE_AGO * offset + k] = M[THIS_ROW * offset + k]
#         end

#         # Reset THIS_ROW
#         for k in 1:offset
#             M[THIS_ROW * offset + k] = 0.0
#         end

#         # Set first column of the row
#         M[THIS_ROW * offset + 1] = M[ONE_AGO * offset + 1] - q_weights[q_pos]

#         for (c_pos, c_func) in enumerate(candidate)
#             col_idx = c_pos +1

#             # Align gap to query position
#             previous_row = M[ONE_AGO * offset + col_idx] - q_weights[q_pos]
#             # println("\tPrevious row= $(M[ONE_AGO * offset + col_idx]) - $(q_weights[q_pos])")

#             # Align gap to candidate position
#             previous_col = M[THIS_ROW * offset + col_idx - 1] - c_weights[c_pos]
#             # println("\tPrevious col= $(M[THIS_ROW * offset + col_idx - 1]) - $(c_weights[c_pos])")

#             # Value of diagonal
#             diagonal_value = M[ONE_AGO * offset + col_idx - 1]

#             if q_func == c_func
#                 #Matching
#                 diagonal_value += max(q_weights[q_pos], c_weights[c_pos]) * exp(-method.γ * abs(q_pos - c_pos))
#             else
#                 #Mismatch
#                 diagonal_value -= q_weights[q_pos] + c_weights[c_pos]
#                 # println("\tDiagonal = $(M[ONE_AGO * offset + col_idx - 1]) - $( q_weights[q_pos] + c_weights[c_pos])")
#             end

#             M[THIS_ROW * offset + col_idx] = max(previous_row, previous_col, diagonal_value)
#             # println("$(q_pos), $(c_pos) max($(previous_row), $(diagonal_value), $(previous_col))")
#         end
#         # println(M[THIS_ROW * offset + 1: THIS_ROW * offset + offset])
#     end

#     sim = M[THIS_ROW * offset + offset]

#     if !method.no_norm
#         return normalize(sim, query, candidate, q_weights, c_weights)
#     end

#     return sim
# end

max_stars_repo_name: irving-muller/FaST

max_issues_repo_name: irving-muller/FaST

max_forks_repo_name: irving-muller/FaST


text:
import Sugar: LazyMethod, supports_overloading, expr_type
import Sugar: getsource!, dependencies!, istype, isfunction, getfuncargs, isintrinsic
import Sugar: typename, functionname, show_name, show_type, show_function

const GLMethod = LazyMethod{:GL}
const GEOMMethod = LazyMethod{:GEOM}
const GLMethods = Union{GLMethod, GEOMMethod}

@compat abstract type AbstractGLIO <: CIO end

immutable GLFunction{Args <: Tuple}
    program::GLProgram
    local_size::NTuple{3, Int}
end
type GLIO{T <: IO} <: AbstractGLIO
    io::T
    method::GLMethods
end

function (::Type{CIO})(io, m::GLMethod)
    GLIO(io, m)
end

supports_overloading(io::GLIO) = true

include("intrinsics.jl")
include("printing.jl")
include("rewriting.jl")

#to_gl_types{T}(::Type{T}) = T
to_gl_types{T}(::T) = to_gl_types(T)
to_gl_types{T}(::Type{T}) = T
to_gl_types{T}(::Type{Type{T}}) = Type{T}
to_gl_types(::Type{Int64}) = Int32
to_gl_types(::Type{Float64}) = Float32

function to_gl_types{T <: Texture}(arg::Type{T})
    return gli.GLTexture{eltype(arg), ndims(arg)}
end
function to_gl_types{T <: GLBuffer}(arg::Type{T})
    return gli.GLArray{eltype(arg), ndims(arg)}
end

function gl_convert{T}(x::T)
    # empty objects are empty and are only usable for dispatch
    isbits(x) && sizeof(x) == 0 && nfields(x) == 0 && return EmptyStruct()
    # same applies for types
    isa(x, Type) && return EmptyStruct()
    convert(to_gl_types(x), x)
end

const gl_compiled_functions = Dict{Any, GLFunction}()

function gl_empty_compile_cache!()
    empty!(gl_compiled_functions)
    return
end

function GLFunction{T}(f::Function, args::T, window)
    gltypes = to_gl_types.(args)
    get!(gl_compiled_functions, (f, gltypes)) do # TODO make this faster
        decl = GLFunction((f, gltypes))
        funcsource = getsource!(decl)
        # add compute program dependant infos
        io = GLIO(IOBuffer(), decl)
        deps = reverse(collect(dependencies!(decl, true)))
        types = filter(istype, deps)
        funcs = filter(isfunction, deps)
        println(io, "// dependant type declarations")
        for typ in types
            if !isintrinsic(typ)
                println(io, getsource!(typ))
            end
        end
        println(io, "// dependant function declarations")
        for func in funcs
            if !isintrinsic(func)
                println(io, getsource!(func))
            end
        end

        println(io, "// Main inner function")
        funcargs = getfuncargs(decl)
        declare_global(io, getfuncargs(decl))
        varnames = map(x-> string(global_identifier, x.args[1]), funcargs)
        print(io, "void main(){\n    ")
        show_name(io, f)
        print(io, "(", join(varnames, ", "), ");\n}")
        shader = GLAbstraction.compile_shader(take!(io.io), GL_COMPUTE_SHADER, Symbol(f))
        program = GLAbstraction.compile_program([shader], [])
        GLFunction{T}(program, local_size)
    end::GLFunction{T}
end

useprogram(p::GLFunction) = glUseProgram(p.program.id)

function (p::GLFunction{Args}){Args}(args::Args, size::NTuple{3})
    useprogram(p)
    for i = 1:length(args)
        bindlocation(args[i], i-1)
    end
    size = ntuple(Val{3}) do i
        ceil(Int, size[i] / p.local_size[i])
    end
    glDispatchCompute(size[1], size[2], size[3])
    return
end


text:
using Documenter, WaterLily

makedocs(
    modules = [WaterLily],
    format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
    authors = "Gabriel Weymouth",
    sitename = "WaterLily.jl",
    pages = Any["index.md"]
    # strict = true,
    # clean = true,
    # checkdocs = :exports,
)

deploydocs(
    repo = "github.com/gabrielweymouth/WaterLily.jl.git",
    push_preview = true
)


text:
export SimpleReplanningRequest

"""
    SimpleReplanningRequest

Intermediate representation of a `ProjectRequest` (useful for I/O)
"""
struct SimpleReplanningRequest
    def::SimpleProblemDef
    t_request::Int
    t_arrival::Int
end
function TOML.parse(def::SimpleReplanningRequest)
    toml_dict = TOML.parse(def.def)
    toml_dict["t_request"] = def.t_request
    toml_dict["t_arrival"] = def.t_arrival
    toml_dict
end

export read_simple_request
function read_simple_request(toml_dict)
    request = SimpleReplanningRequest(
        read_problem_def(toml_dict),
        toml_dict["t_request"],
        toml_dict["t_arrival"]
    )
end
read_simple_request(path::String) = read_simple_request(TOML.parsefile(path))

function ProjectRequest(def::SimpleReplanningRequest,prob_spec)
    ProjectRequest(
        construct_partial_project_schedule(
            def.def.project_spec,prob_spec),
        def.t_request,
        def.t_arrival
        )
    set_t0!(request.schedule,t_request)
    process_schedule!(request.schedule)
    return request.schedule
end

export SimpleRepeatedProblemDef

"""
    SimpleRepeatedProblemDef

Intermediate representation of a `RepeatedPC_TAPF` (useful for I/O)
"""
@with_kw struct SimpleRepeatedProblemDef
    requests::Vector{SimpleReplanningRequest} = Vector{SimpleReplanningRequest}()
    r0::Vector{Int} = Int[]
    env_id::String = ""
end

export write_simple_repeated_problem_def
# problem_path/
#   config.toml - # robot initial positions, env_id
#   stages/
#       stage1/
#           problem.toml
#       ...
#       stageN/
#           problem.toml
function write_simple_repeated_problem_def(path,def::SimpleRepeatedProblemDef)
    toml_dict = Dict("r0" => def.r0,"env_id" => def.env_id)
    mkpath(path)
    open(joinpath(path,"def.toml"),"w") do io
        TOML.print(io,toml_dict)
    end
    for (i,stage_def) in enumerate(def.requests)
        p = joinpath(path,@sprintf("stage%3.3i",i))
        mkpath(p)
        toml_dict = TOML.parse(stage_def)
        open(joinpath(p,"problem.toml"),"w") do io
            TOML.print(io,toml_dict)
        end
    end
    return true
end

export read_simple_repeated_problem_def
function read_simple_repeated_problem_def(path)
    toml_dict = TOML.parsefile(joinpath(path,"def.toml"))
    prob_def = SimpleRepeatedProblemDef(
        SimpleReplanningRequest[],
        toml_dict["r0"],
        toml_dict["env_id"],
    )
    for p in readdir(path;join=true,sort=true)
        isdir(p) ? nothing : continue
        # prob_dict = TOML.parsefile(joinpath(p,"problem.toml"))
        # push!(prob_def.requests,read_simple_request(prob_dict))
        push!(prob_def.requests,read_simple_request(joinpath(p,"problem.toml")))
    end
    prob_def
end
function RepeatedPC_TAPF(simple_def::SimpleRepeatedProblemDef,solver,loader::ReplanningProblemLoader)
    env         = loader.envs[simple_def.env_id]
    prob_spec   = loader.prob_specs[simple_def.env_id]
    sched = construct_partial_project_schedule(
        [ROBOT_AT(r,x) for (r,x) in enumerate(simple_def.r0)],prob_spec)
    search_env = construct_search_env(solver,sched,prob_spec,env)
    r_pc_tapf = RepeatedPC_TAPF(
        search_env,
        map(req->ProjectRequest(req,prob_spec), simple_def.requests)
    )
end

"""
    instantiate_random_repeated_pctapf_problem(env,config)

Instantiate a random `RepeatedPC_TAPF` problem based on the parameters of config.
"""
function random_pctapf_sequence(env::GridFactoryEnvironment,config;
        num_projects        = get(config,:num_projects,10),
        kwargs...)
    projects = map(i->random_pctapf_def(env,config;kwargs...), 1:num_projects)
end

function random_repeated_pctapf_def(env,config;
        arrival_interval    = get(config,:arrival_interval,40),
        warning_time        = get(config,:warning_time,0),
        N                   = get(config,:N,0),
        env_id              = get(config,:env_id,""),
        kwargs...)
    projects = random_pctapf_sequence(env,config;N=0,kwargs...)
    requests = map(i->SimpleReplanningRequest(
        projects[i],
        i*arrival_interval, # request
        i*arrival_interval + warning_time # arrival time
        ), 1:length(projects))
    r0,_,_ = get_random_problem_instantiation(N,0,[],[],get_free_zones(env))
    prob_def = SimpleRepeatedProblemDef(
        requests = requests,
        r0 = r0,
        env_id = env_id
        )
end

export write_problems!

padded_problem_name(number,name="problem",ext=".toml") = @sprintf("%s%4.4i%s",name,number,ext)
init_random_problem(loader::ReplanningProblemLoader,env,config) = random_repeated_pctapf_def(env,config)
write_problem(loader::ReplanningProblemLoader,problem_def,prob_file,args...) = write_simple_repeated_problem_def(prob_file,problem_def)

export compile_replanning_results!

function compile_replanning_results!(
        cache::ReplanningProfilerCache,solver,env,
        timer_results,prob,stage,request)
    push!(cache.schedules,deepcopy(get_schedule(env)))
    # store project ids
    for v in vertices(request.schedule)
        id = get_vtx_id(request.schedule,v)
        cache.project_ids[id] = stage
    end
    # store results
    stage_results = compile_results(solver,cache.features,prob,env,timer_results)
    if debug(solver) && !feasible_status(solver)
        merge!(stage_results, compile_results(solver,cache.final_features,prob,env,timer_results))
    end
    push!(cache.stage_results, stage_results)
    # store final results
    if stage == length(prob.requests)
        merge!(cache.final_results, compile_results(solver,cache.final_features,prob,env,timer_results))
    end
    return cache
end

export write_replanning_results

# function write_replanning_results(cache::ReplanningProfilerCache,results_path)
#     mkpath(results_path)
#     for (i,results) in enumerate(cache.stage_results)
#         path = joinpath(results_path,padded_problem_name(i,"stage"))
#         open(path,"w") do io
#             TOML.print(io,results)
#         end
#     end
#     open(joinpath(results_path,"final_results.toml"),"w") do io
#         TOML.print(io,cache.final_results)
#     end
#     return results_path
# end

function write_replanning_results(
    loader::ReplanningProblemLoader,
    planner::ReplannerWithBackup,
    results_path,
    primary_prefix="primary_planner",
    backup_prefix="backup_planner",
    )
    results = Dict()
    for (pre,cache) in [
            (backup_prefix,planner.backup_planner.cache),
            (primary_prefix,planner.primary_planner.cache),
        ]

        results[pre]  = cache.stage_results
        results[string(pre,"_final")] = cache.final_results
    end
    mkpath(results_path)
    open(joinpath(results_path,"results.toml"),"w") do io
        TOML.print(io,results)
    end
    # write_replanning_results(planner.primary_planner.cache,
    #     joinpath(results_path,primary_prefix))
    # write_replanning_results(planner.backup_planner.cache,
    #     joinpath(results_path,backup_prefix))
    return results_path
end

function load_replanning_results(loader::ReplanningProblemLoader,solver::ReplannerWithBackup,results_path,
    # primary_prefix="primary_planner",
    # backup_prefix="backup_planner",
    )
    results = TOML.parsefile(joinpath(results_path,"results.toml"))
    results["problem_name"] = CRCBS.get_problem_name(loader,results_path)
    return results
    # results_dict = Dict(
    #     p=>load_replanning_results(loader,joinpath(results_path,p)) for p in [
    #         primary_prefix,backup_prefix
    #     ]
    # )
end
function CRCBS.load_results(loader::ReplanningProblemLoader,results_path)
    dict = TOML.parsefile(joinpath(results_path,"results.toml"))
    results = to_symbol_dict(dict)
    results[:problem_name] = CRCBS.get_problem_name(loader,results_path)
    return results
end

export profile_replanner!

function is_problem_file(loader::ReplanningProblemLoader,path)
    if isdir(path) # each replanning problem is in a directory
        prob_name = splitdir(path)[end]
        if !isnothing(findfirst("problem",prob_name))
            return true
        end
    end
    return false
end
function profile_replanner!(
    loader::ReplanningProblemLoader,
    planner::ReplannerWithBackup,
    base_problem_dir,
    base_results_dir,
    )
    for f in readdir(base_problem_dir;join=true)
        if is_problem_file(loader,f)
            problem_name = splitdir(f)[end]
            outpath = joinpath(base_results_dir,problem_name)
            ispath(outpath) ? continue : nothing # skip over this if results already exist
            @log_info(-1,0,"PROFILING: testing on problem ",problem_name)
            simple_prob_def = read_simple_repeated_problem_def(f)
            prob = RepeatedPC_TAPF(simple_prob_def,planner.primary_planner.solver,loader)
            search_env, planner = profile_replanner!(planner,prob)
            @log_info(-1,0,"PROFILING: writing results for problem ",problem_name)
            write_replanning_results(loader,planner,outpath)
        end
    end
end

function profile_replanner!(planner::FullReplanner,prob::RepeatedAbstractPC_TAPF)
    profile_replanner!(planner.solver,planner.replanner,prob,planner.cache)
end
function profile_replanner!(solver,replan_model,prob::RepeatedAbstractPC_TAPF,
        cache = ReplanningProfilerCache()
    )
    env = get_env(prob)
    for (stage,request) in enumerate(prob.requests)
        @log_info(1,verbosity(solver),"REPLANNING: Stage ",stage)
        remap_object_ids!(request.schedule,get_schedule(env))
        base_env = replan!(solver,replan_model,env,request)
        reset_solver!(solver)
        # env, cost = solve!(solver,base_env)
        env, timer_results = profile_solver!(solver,construct_routing_problem(prob,base_env))
        compile_replanning_results!(cache,solver,env,timer_results,prob,stage,request)
        @log_info(1,verbosity(solver),"Stage ",stage," - ","route planner iterations: ",
            iterations(route_planner(solver)))
    end
    return env, cache
end

function profile_replanner!(planner::ReplannerWithBackup,prob::RepeatedAbstractPC_TAPF)
    reset_cache!(planner)
    plannerA = planner.primary_planner
    plannerB = planner.backup_planner
    hard_reset_solver!(planner.primary_planner.solver)
    hard_reset_solver!(planner.backup_planner.solver)
    env = get_env(prob)
    for (stage,request) in enumerate(prob.requests)
        remap_object_ids!(request.schedule,get_schedule(env))

        base_envB = replan!(plannerB,env,request)
        envB, resultsB = profile_solver!(plannerB.solver,construct_routing_problem(prob,base_envB))
        compile_replanning_results!(plannerB.cache,plannerB.solver,envB,
            resultsB,prob,stage,request)

        base_envA = replan!(plannerA,env,request)
        envA, resultsA = profile_solver!(plannerA.solver,construct_routing_problem(prob,base_envA))
        compile_replanning_results!(plannerA.cache,plannerA.solver,envA,
            resultsA,prob,stage,request)

        # if failed_status(plannerA) == false
        if feasible_status(plannerA)
            @log_info(-1,0,"REPLANNING: ",
            "Primary planner succeeded at stage $stage.")
            env = envA
        elseif feasible_status(plannerB)
            @log_info(-1,0,"REPLANNING: ",
            "Primary planner failed at stage $stage. Proceeding with backup plan.")
            env = envB
            validate(get_schedule(envB))
        else
            @log_info(-1,0,"REPLANNING:",
                "Both primary and backup planners failed at stage $stage.",
                " Returning early.")
            break
        end
    end
    return env, planner
end

function warmup(planner::ReplannerWithBackup,loader::ReplanningProblemLoader,
        base_dir            = joinpath("/tmp","warmup"),
        base_problem_dir    = joinpath(base_dir,"problem_instances"),
        base_results_dir    = joinpath(base_dir,"results"),
        config = Dict(
            :N => 10,
            :M => 5,
            :num_projects => 3,
            :env_id => sort(collect(keys(loader.envs)))[1],
            ),
    )
    Random.seed!(0)
    write_problems!(loader,config,base_problem_dir)
    env = get_env!(loader,config[:env_id])
    # turn off time constraints for warm up
    set_real_time_flag!(planner,false)
    profile_replanner!(loader,planner,base_problem_dir,base_results_dir)
    # turn real time constraints back on
    set_real_time_flag!(planner,true)
    run(`rm -rf $base_dir`) # remove dummy warmup folder
    reset_cache!(planner)
    return planner
end

export product_config_dicts

function product_config_dicts(configs...)
    map(dicts->merge(dicts...),Base.Iterators.product(configs...))[:]
end

export
    replanning_config_1,
    replanning_config_2,
    replanning_config_3,
    replanning_config_4,
    replanning_config_5

"""
    replanning_config_1

Returns a vector of config dictionaries, which can be used to generate random
problem instances for profiling.
"""
function replanning_config_1()
    base_configs = [
        Dict(
            :warning_time=>0,
            :commit_threshold=>10,
            :fallback_commit_threshold=>10,
            :num_trials => 4,
            :max_parents => 3,
            :depth_bias => 0.4,
            :dt_min => 0,
            :dt_max => 0,
            :dt_collect => 0,
            :dt_deliver => 0,
            # :task_sizes => (1=>1.0,2=>0.0,4=>0.0),
            :env_id=>"env_2",
            )
    ]
    stream_configs = [
        Dict(:N=>30, :M=>10, :num_projects=>10, :arrival_interval=>40, ),
        Dict(:N=>30, :M=>15, :num_projects=>10, :arrival_interval=>50, ),
        Dict(:N=>30, :M=>20, :num_projects=>10, :arrival_interval=>60, ),
        Dict(:N=>30, :M=>25, :num_projects=>10, :arrival_interval=>70, ),
        Dict(:N=>30, :M=>30, :num_projects=>10, :arrival_interval=>80, ),
    ]
    product_config_dicts(base_configs,stream_configs)
end
function replanning_config_2()
    base_configs = [
        Dict(
            :warning_time=>4,
            :commit_threshold=>4,
            :fallback_commit_threshold=>4,
            :num_trials => 1,
            :max_parents => 3,
            :depth_bias => 0.4,
            :dt_min => 0,
            :dt_max => 0,
            :dt_collect => 0,
            :dt_deliver => 0,
            # :task_sizes => (1=>1.0,2=>0.0,4=>0.0),
            :save_paths => false,
            :env_id=>"env_2",
            )
    ]
    stream_configs = [
        Dict(:N=>2, :M=>3, :num_projects=>4, :arrival_interval=>10, ),
    ]
    product_config_dicts(base_configs,stream_configs)
end
function replanning_config_3()
    base_configs = [
        Dict(
            :warning_time=>0,
            :commit_threshold=>10,
            :fallback_commit_threshold=>10,
            :num_trials => 4,
            :max_parents => 3,
            :depth_bias => 0.4,
            :dt_min => 0,
            :dt_max => 0,
            :dt_collect => 0,
            :dt_deliver => 0,
            # :task_sizes => (1=>1.0,2=>0.0,4=>0.0),
            :env_id=>"env_2",
            )
    ]
    robot_configs = [
        Dict(:N=>30)
    ]
    project_configs = [
        Dict(:num_projects=>10),
        Dict(:num_projects=>20),
        Dict(:num_projects=>30),
    ]
    stream_configs = [
        Dict(:M=>10, :arrival_interval=>20, ),
        Dict(:M=>10, :arrival_interval=>30, ),
        Dict(:M=>10, :arrival_interval=>40, ),

        Dict(:M=>15, :arrival_interval=>30, ),
        Dict(:M=>15, :arrival_interval=>40, ),
        Dict(:M=>15, :arrival_interval=>50, ),

        Dict(:M=>20, :arrival_interval=>40, ),
        Dict(:M=>20, :arrival_interval=>50, ),
        Dict(:M=>20, :arrival_interval=>60, ),

        Dict(:M=>25, :arrival_interval=>50, ),
        Dict(:M=>25, :arrival_interval=>60, ),
        Dict(:M=>25, :arrival_interval=>70, ),

        Dict(:M=>30, :arrival_interval=>60, ),
        Dict(:M=>30, :arrival_interval=>70, ),
        Dict(:M=>30, :arrival_interval=>80, ),
    ]
    product_config_dicts(base_configs,robot_configs,project_configs,stream_configs)
end
function replanning_config_4()
    base_configs = [
        Dict(
            :warning_time=>0,
            :commit_threshold=>10,
            :fallback_commit_threshold=>10,
            :num_trials => 8,
            :max_parents => 3,
            :depth_bias => 0.4,
            :dt_min => 0,
            :dt_max => 0,
            :dt_collect => 0,
            :dt_deliver => 0,
            # :task_sizes => (1=>1.0,2=>0.0,4=>0.0),
            :env_id=>"env_2",
            )
    ]
    robot_configs = [
        Dict(:N=>30)
    ]
    project_configs = [
        Dict(:num_projects=>30),
    ]
    stream_configs = [
        # Dict(:M=>10, :arrival_interval=>20, ),
        # Dict(:M=>10, :arrival_interval=>30, ),
        # Dict(:M=>10, :arrival_interval=>40, ),

        Dict(:M=>15, :arrival_interval=>20, ),
        # Dict(:M=>15, :arrival_interval=>30, ),
        # Dict(:M=>15, :arrival_interval=>40, ),
        # Dict(:M=>15, :arrival_interval=>50, ),

        Dict(:M=>20, :arrival_interval=>20, ),
        Dict(:M=>20, :arrival_interval=>30, ),
        # Dict(:M=>20, :arrival_interval=>40, ),
        # Dict(:M=>20, :arrival_interval=>50, ),
        # Dict(:M=>20, :arrival_interval=>60, ),

        Dict(:M=>25, :arrival_interval=>30, ),
        Dict(:M=>25, :arrival_interval=>40, ),
        # Dict(:M=>25, :arrival_interval=>50, ),
        # Dict(:M=>25, :arrival_interval=>60, ),
        # Dict(:M=>25, :arrival_interval=>70, ),

        Dict(:M=>30, :arrival_interval=>40, ),
        Dict(:M=>30, :arrival_interval=>50, ),
        # Dict(:M=>30, :arrival_interval=>60, ),
        # Dict(:M=>30, :arrival_interval=>70, ),
        # Dict(:M=>30, :arrival_interval=>80, ),
    ]
    product_config_dicts(base_configs,robot_configs,project_configs,stream_configs)
end
function replanning_config_5()
    base_configs = [
        Dict(
            :warning_time=>0,
            :commit_threshold=>10,
            :fallback_commit_threshold=>10,
            :num_trials => 8,
            :max_parents => 3,
            :depth_bias => 0.4,
            :dt_min => 0,
            :dt_max => 0,
            :dt_collect => 0,
            :dt_deliver => 0,
            :env_id=>"env_2",
            )
    ]
    robot_configs = [
        Dict(:N=>30)
    ]
    project_configs = [
        Dict(:num_projects=>30),
    ]
    stream_configs = [
        Dict(:M=>10, :arrival_interval=>20, ),
        Dict(:M=>10, :arrival_interval=>30, ),
        Dict(:M=>10, :arrival_interval=>40, ),

        Dict(:M=>15, :arrival_interval=>20, ),
        Dict(:M=>15, :arrival_interval=>30, ),
        Dict(:M=>15, :arrival_interval=>40, ),
        Dict(:M=>15, :arrival_interval=>50, ),

        Dict(:M=>20, :arrival_interval=>20, ),
        Dict(:M=>20, :arrival_interval=>30, ),
        Dict(:M=>20, :arrival_interval=>40, ),
        Dict(:M=>20, :arrival_interval=>50, ),
        Dict(:M=>20, :arrival_interval=>60, ),

        Dict(:M=>25, :arrival_interval=>20, ),
        Dict(:M=>25, :arrival_interval=>30, ),
        Dict(:M=>25, :arrival_interval=>40, ),
        Dict(:M=>25, :arrival_interval=>50, ),
        Dict(:M=>25, :arrival_interval=>60, ),
        Dict(:M=>25, :arrival_interval=>70, ),

        Dict(:M=>30, :arrival_interval=>20, ),
        Dict(:M=>30, :arrival_interval=>30, ),
        Dict(:M=>30, :arrival_interval=>40, ),
        Dict(:M=>30, :arrival_interval=>50, ),
        Dict(:M=>30, :arrival_interval=>60, ),
        Dict(:M=>30, :arrival_interval=>70, ),
        Dict(:M=>30, :arrival_interval=>80, ),
    ]
    product_config_dicts(base_configs,robot_configs,project_configs,stream_configs)
end

export setup_replanning_experiments

function setup_replanning_experiments(base_problem_dir,base_results_dir)
    loader = ReplanningProblemLoader()
    add_env!(loader,"env_2",init_env_2())
    # get probem config
    problem_configs = replanning_config_3()
    # write problems
    Random.seed!(0)
    # Problem generation and profiling
    reset_task_id_counter!()
    reset_operation_id_counter!()
    write_problems!(loader,problem_configs,base_problem_dir)
    # set up profiling features
    feats = [
        RunTime(),IterationCount(),LowLevelIterationCount(),
        TimeOutStatus(),IterationMaxOutStatus(),
        SolutionCost(),OptimalityGap(),OptimalFlag(),FeasibleFlag(),NumConflicts(),
        ObjectPathSummaries(),
        ]
    final_feats = [SolutionCost(),NumConflicts(),RobotPaths()]
    # set up planners
    planners = Vector{ReplannerWithBackup}()
    for (primary_replanner, backup_replanner) in [
            (MergeAndBalance(),MergeAndBalance()),
            (ReassignFreeRobots(),ReassignFreeRobots()),
            (DeferUntilCompletion(),DeferUntilCompletion()),
        ]
        # Primary planner
        path_finder = DefaultAStarSC()
        set_iteration_limit!(path_finder,5000)
        primary_route_planner = CBSSolver(ISPS(path_finder))
        set_iteration_limit!(primary_route_planner,5000)
        primary_planner = FullReplanner(
            solver = NBSSolver(path_planner=primary_route_planner),
            replanner = primary_replanner,
            cache = ReplanningProfilerCache(features=feats,final_features=final_feats)
            )
        # Backup planner
        backup_planner = FullReplanner(
            solver = NBSSolver(
                assignment_model = TaskGraphsMILPSolver(GreedyAssignment()),
                path_planner = PIBTPlanner{NTuple{3,Float64}}()
                ),
            replanner = backup_replanner,
            cache = ReplanningProfilerCache(features=feats,final_features=final_feats)
            )
        set_iteration_limit!(backup_planner,1)
        set_iteration_limit!(route_planner(backup_planner.solver),5000)
        # Full solver
        planner = ReplannerWithBackup(primary_planner,backup_planner)
        push!(planners,planner)
    end
    return loader, planners
end

"""
    add start_time, completion_time, and makespan for each stage
"""
function post_process_replanning_results!(results,config)
    M = config[:M]
    for k in [:primary_planner,:backup_planner]
        for (i,stage_dict) in enumerate(results[k])
            id_range = (1+(i-1)*M,i*M)
            stage_dict[:start_time] = i*config[:arrival_interval]
            completion_time = 0
            if stage_dict[:FeasibleFlag]
                for (id_string,object_summary) in stage_dict[:ObjectPathSummaries]
                    if id_range[1] <= object_summary[:object_id] <= id_range[2]
                        completion_time = max(completion_time,object_summary[:deposit_time])
                    end
                end
                # Make sure object path summaries are up to date
                results[:ObjectPathSummaries] = merge!(
                    get(results,:ObjectPathSummaries,Dict{Symbol,Any}()),
                    stage_dict[:ObjectPathSummaries]
                    )
                if i == length(results[k])
                    kp = Symbol(string(k,"_final"))
                    results[:RobotPaths] = results[kp][:RobotPaths]
                end
            else
                completion_time = typemax(Int)
            end
            stage_dict[:completion_time] = completion_time
            stage_dict[:makespan] = completion_time - stage_dict[:start_time]
        end
    end
    makespans = Int[]
    arrival_times = Int[]
    backup_flags = Bool[]
    primary_runtimes = Float64[]
    backup_runtimes = Float64[]
    runtime_gaps = Float64[]
    for (resultsA,resultsB) in zip(results[:primary_planner],results[:backup_planner])
        push!(primary_runtimes,resultsA[:RunTime])
        push!(backup_runtimes,resultsB[:RunTime])
        push!(runtime_gaps,(resultsA[:RunTime])/(resultsA[:RunTime]+resultsB[:RunTime]))
        if resultsA[:FeasibleFlag] == true
            push!(makespans,resultsA[:makespan])
            push!(arrival_times,resultsA[:start_time])
            push!(backup_flags,false)
        elseif resultsB[:FeasibleFlag] == true
            push!(makespans,resultsB[:makespan])
            push!(arrival_times,resultsB[:start_time])
            push!(backup_flags,true)
        else
            break
        end
    end
    results[:makespans] = makespans
    results[:arrival_times] = arrival_times
    results[:backup_flags] = backup_flags
    results[:runtime_gaps] = runtime_gaps
    results[:primary_runtimes] = primary_runtimes
    results[:backup_runtimes] = backup_runtimes
    results[:completion_times] = arrival_times .+ makespans
    return results
end

function construct_replanning_results_dataframe(loader,solver_config,feats)
    results_df = init_dataframe(feats)
    for results_file in readdir(solver_config.results_path;join=true)
        results = load_results(loader,results_file)
        prob_name = results[:problem_name]
        prob_file = joinpath(solver_config.problem_path,prob_name)
        config = CRCBS.load_config(loader,prob_file)
        TaskGraphs.post_process_replanning_results!(results,config)
        push!(results_df,results)
    end
    results_df
end


text:
module BenchmarkingLinearInterpolators

using Printf
using LinearInterpolators
using LazyAlgebra
using BenchmarkTools

mimeshow(io::IO, x) = show(io, MIME"text/plain"(), x)
mimeshow(x) = mimeshow(stdout, x)

prt(op::AbstractString, t) = begin
    @printf("%s:\n", op)
    mimeshow(t)
end

prt(op::AbstractString, t, nops::Integer) = begin
    # times in nanoseconds, speed in Gigaflops
    @printf("%s [%.6f Gflops]:\n", op, nops/minimum(t.times))
    mimeshow(t)
    println()
end

n1, n2 = 300, 200
m1, m2 =  43,  29
y = randn(n1, n2)
ker = CatmullRomSpline()
A1 = SparseUnidimensionalInterpolator(ker, 1, 1:n1, range(1,n1,length=m1))
A2 = SparseUnidimensionalInterpolator(ker, 2, 1:n2, range(1,n2,length=m2))
H = A1*A2

function LazyAlgebra.vcreate(::Type{LazyAlgebra.Direct},
                             ::Gram,
                             x::AbstractArray{T,N},
                             scratch::Bool) where {T<:AbstractFloat,N}
    return similar(x)
end

r1 = 2*length(A1.C)
r2 = 2*length(A2.C)
x = conjgrad(H'*H, H'*y; verb=true, ftol=1e-8, maxiter=50)
z = vcopy(y)
prt("A1*x",      @benchmark($A1*$x),        r1*m2)
prt("A2*x",      @benchmark($A2*$x),        r2*m1)
prt("A2*A1*x",   @benchmark($(A2*A1)*$x),   r1*m2 + r2*n1)
prt("A1*A2*x",   @benchmark($(A1*A2)*$x),   r2*m1 + r1*n2)
prt("A1'*y",     @benchmark($(A1')*$y),     r1*m2)
prt("A2'*y",     @benchmark($(A2')*$y),     r2*m1)
prt("A2'*A1'*y", @benchmark($(A2'*A1')*$y), r1*m2 + r2*n1)
prt("A1'*A2'*y", @benchmark($(A1'*A2')*$y), r2*m1 + r1*n2)
println("")

n = 100
pos = [1.1*i + 0.2*j for i in 1:n1, j in 1:n2]
pmin, pmax = extrema(pos)
A = SparseInterpolator(ker, pos, range(pmin, pmax, length=n))
x = randn(n)
y = randn(n1,n2)
nops = 2*length(A.C)

prt("A*x",  @benchmark($A*$x), nops)
prt("A'*x", @benchmark($(A')*$y), nops)

end # module


text:
# # Brusselator
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/models/Brusselator.ipynb)
#
#md # !!! note "Overview"
#md #     System type: polynomial continuous system\
#md #     State dimension: 2\
#md #     Application domain: Chemical kinetics
#
# ## Model description
#
# A chemical reaction is said to be *autocatalytic* if one of the reaction products is
# also a catalyst for the same or a coupled reaction, and such a reaction is called an autocatalytic reaction.
# We refer to the wikipedia article [Autocatalysis](https://en.wikipedia.org/wiki/Autocatalysis) for details.

# The Brusselator is a mathematical model for a class of autocatalytic reactions.
# The dynamics of the Brusselator is given by the two-dimensional ODE
#
# ```math
#   \left\{ \begin{array}{lcl} \dot{x} & = & A + x^2\cdot y - B\cdot x - x \\
#    \dot{y} & = & B\cdot x - x^2\cdot y \end{array} \right.
# ```

# The numerical values for the model's constants (in their respective units) are
# given in the following table.
#
# |Quantity|Value|
# |-----|-------|
# |A    | 1     |
# |B    | 1.5   |

using ReachabilityAnalysis #!jl

const A = 1.0
const B = 1.5
const B1 = B + 1

@taylorize function brusselator!(du, u, p, t)
    x, y = u
    x² = x * x
    aux = x² * y
    du[1] = A + aux - B1*x
    du[2] = B*x - aux
    return du
end

#md # !!! tip "Performance tip"
#md #     The auxiliary variables `B1`, `x²` and `aux` have been defined to make
#md #     better use of `@taylorize` and help to reduce allocations.

# ## Reachability settings
#
# The initial set is defined by ``x \in [0.8, 1]``, ``y \in [0, 0.2]``.
# These settings are taken from [^CAS13].

U₀ = (0.8 .. 1.0) × (0.0 .. 0.2); #!jl
prob = @ivp(u' = brusselator!(U), u(0) ∈ U₀, dim: 2); #!jl

# ## Results

# We use `TMJets` algorithm with sixth-order expansion in time and second order expansion
# in the spatial variables.

sol = solve(prob, T=18.0, alg=TMJets20(orderT=6, orderQ=2)); #!jl

#-

using Plots #!jl

plot(sol, vars=(1, 2), xlab="x", ylab="y", lw=0.2, color=:blue, lab="Flowpipe", legend=:bottomright) #!jl
plot!(U₀, color=:orange, lab="Uo") #!jl

# We observe that the system converges to the equilibrium point `(1.0, 1.5)`.

# Below we plot the flowpipes projected into the time domain.

plot(sol, vars=(0, 1), xlab="t", lw=0.2, color=:blue, lab="x(t)", legend=:bottomright)  #!jl
plot!(sol, vars=(0, 2), xlab="t", lw=0.2, color=:red, lab="y(t)")  #!jl

# ## Changing the initial volume

# The model was considered in [^GCLASG20] but using a different set of initial conditions. Let us parametrize
# the initial states as a ball centered at ``x = y = 1`` and radius ``r > 0``:

U0(r) = Singleton([1.0, 1.0]) ⊕ BallInf(zeros(2), r)

#-

# The parametric initial-value problem is defined accordingly.

bruss(r) = @ivp(u' = brusselator!(u), u(0) ∈ U0(r), dim: 2)

# First we solve for ``r = 0.01``:

sol_01 = solve(bruss(0.01), T=30.0, alg=TMJets20(orderT=6, orderQ=2))  #!jl

LazySets.set_ztol(Float64, 1e-15)  #!jl

plot(sol_01, vars=(1, 2), xlab="x", ylab="y", lw=0.2, color=:blue, lab="Flowpipe (r = 0.01)", legend=:bottomright)  #!jl

plot!(U0(0.01), color=:orange, lab="Uo", xlims=(0.6, 1.3))  #!jl

# We observe that the wrapping effect is controlled and the flowpipe doesn't blow up even for the large time horizon ``T = 30.0``.
# Next we plot the flowpipe zoomed to the last portion and compare ``r = 0.01`` with a set of larger initial states, ``r = 0.1``.

sol_1 = solve(bruss(0.1), T=30.0, alg=TMJets20(orderT=6, orderQ=2))

fig = plot(xlab="x", ylab="y", xlims=(0.9, 1.05), ylims=(1.43, 1.57), legend=:bottomright)  #!jl

plot(sol_1, vars=(1, 2), lw=0.2, color=:red, lab="r = 0.1", alpha=.4)  #!jl

plot!(sol_01, vars=(1, 2), lw=0.2, color=:blue, lab="r = 0.01")  #!jl

# The volume at time ``T = 9.0`` can be obtained by evaluating the flowpipe and computing the volume of the hyperrectangular overapproximation:

vol_01 = overapproximate(sol_01(9.0), Hyperrectangle) |> set |> volume  #!jl

#-

vol_1 = overapproximate(sol_1(9.0), Hyperrectangle) |> set |> volume  #!jl

# ## References

# [^CAS13]: X. Chen, E. Abraham, S. Sankaranarayanan. *Flow*: An Analyzer for Non-Linear Hybrid Systems.* In Proceedings of the 25th International Conference on Computer Aided Verification (CAV’13). Volume 8044 of LNCS, pages 258-263, Springer, 2013.

# [^GCLASG20]: S. Gruenbacher, J. Cyranka, M. Lechner, Md. Ariful Islam,, Scott A. Smolka, R. Grosu *Lagrangian Reachtubes: The Next Generation*. arXiv: 2012.07458


text:
module Stheno

    using Reexport

    using AbstractGPs
    using BlockArrays
    using ChainRulesCore
    using FillArrays
    @reexport using KernelFunctions
    using LinearAlgebra
    using Random
    using Zygote
    using ZygoteRules

    import Base.Broadcast: broadcasted

    using AbstractGPs: AbstractGP, FiniteGP, GP
    import AbstractGPs:
        mean,
        cov,
        var,
        mean_and_cov,
        mean_and_var,
        rand,
        logpdf,
        elbo,
        dtc,
        posterior,
        approx_posterior,
        marginals

    using MacroTools: @capture, combinedef, postwalk, splitdef

    const AV{T} = AbstractVector{T}

    # Various bits of utility that aren't inherently GP-related. Often very type-piratic.
    include(joinpath("util", "zygote_rules.jl"))
    include(joinpath("util", "covariance_matrices.jl"))
    include(joinpath("util", "block_arrays", "dense.jl"))
    include(joinpath("util", "block_arrays", "diagonal.jl"))
    include(joinpath("util", "abstract_data_set.jl"))
    include(joinpath("util", "proper_type_piracy.jl"))

    # Supertype for GPs.
    include("abstract_gp.jl")

    # Atomic GP objects.
    include(joinpath("gp", "gp.jl"))

    # Composite GPs, constructed via affine transformation of CompositeGPs and GPs.
    include(joinpath("composite", "composite_gp.jl"))
    include(joinpath("composite", "cross.jl"))
    include(joinpath("composite", "product.jl"))
    include(joinpath("composite", "addition.jl"))
    include(joinpath("composite", "compose.jl"))
    # include(joinpath("composite", "gradient.jl"))
    # include(joinpath("composite", "integrate.jl"))

    # Gaussian Process Probabilistic Programme object which implements the AbstractGPs API.
    include("gaussian_process_probabilistic_programme.jl")

    # Sparse GP hack to make pseudo-point approximations play nicely with Turing.jl.
    include("sparse_finite_gp.jl")

    include("deprecate.jl")

    export wrap, BlockData, GPC, GPPPInput, @gppp
    export ∘, select, stretch, periodic, shift
    export cov_diag, mean_and_cov_diag
end # module


text:
using .ARCSolver
using Test

@testset "ARCSolver.jl" begin
    # Write your tests here.
end


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Many functions to perform a preliminary analysis of an optical payload.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export minimum_swath_grss, minimum_half_FOV_grss, swath_width

"""
    function minimum_half_FOV_grss(h::Real, T::Real, i::Real, To::Integer)

Compute the minimum half FOV of a ground repeating Sun-synchronous (GRSS) orbit
to cover the entire Equator within the revisit interval.

# Args

* `h`: Orbit altitude in the Equator [m].
* `T`: Orbit period [s].
* `i`: Inclination [rad].
* `To`: Orbit cycle [days].

# Returns

The minimum half FOV [rad].

"""
function minimum_half_FOV_grss(h::Real, T::Real, i::Real, To::Integer)
    adjacent_track_angle_grss(h, T, i, To, 0.0)
end

"""
    function minimum_half_FOV_grss(h::Real, a::Real, e::Real, i::Real, To::Integer)

Compute the minimum half FOV of a ground repeating Sun-synchronous (GRSS) orbit
to cover the entire Equator within the revisit interval.

# Args

* `h`: Orbit altitude in the Equator [m].
* `a`: Semi-major axis [m].
* `e`: Eccentricity.
* `i`: Inclination [rad].
* `To`: Orbit cycle [days].

# Returns

The minimum half FOV [rad].

"""
function minimum_half_FOV_grss(h::Real, a::Real, e::Real, i::Real, To::Integer)
    adjacent_track_angle_grss(h, a, e, i, To, 0.0)
end

"""
    function minimum_swath_grss(T::Real, i::Real, To::Integer)

Compute the minimum swath of a ground repeating Sun-synchronous (GRSS) orbit to
cover the entire Equator within the revisit interval.

# Args

* `T`: Orbit period [s].
* `i`: Inclination [rad].
* `To`: Orbit cycle [days].

# Returns

The minimum swath [m].

"""
function minimum_swath_grss(T::Real, i::Real, To::Integer)
    adjacent_track_distance_grss(T, i, To, 0.0)
end

"""
    function minimum_swath_grss(a::Real, e::Real, i::Real, To::Integer)

Compute the minimum swath of a ground repeating Sun-synchronous (GRSS) orbit to
cover the entire Equator within the revisit interval.

# Args

* `a`: Semi-major axis [m].
* `e`: Eccentricity.
* `i`: Inclination [rad].
* `To`: Orbit cycle [days].

# Returns

The minimum swath [m].

"""
function minimum_swath_grss(a::Real, e::Real, i::Real, To::Integer)
    adjacent_track_distance_grss(a, e, i, To, 0.0)
end

"""
    function swath_width(h::real, HalfFOV::real)

Compute the swath width given the orbit altitude and the half FOV.

# Args

* `h`: Orbit altitude [m].
* `HalfFOV`: Half field of view [rad].

# Returns

The swath width [m].

"""
function swath_width(h::Real, HalfFOV::Real)
    gamma = pi - asin((R0+h)/R0*sin(HalfFOV))
    alpha = pi - gamma - HalfFOV
    S = R0*alpha
end


text:
# REINFORCE loss
function reinforce_loss(π, 𝒫, 𝒟; info = Dict())
    new_probs = logpdf(π, 𝒟[:s], 𝒟[:a])

    ignore() do
        info[:entropy] = mean(entropy(π, 𝒟[:s]))
        info[:kl] = mean(𝒟[:logprob] .- new_probs)
    end

    -mean(new_probs .* 𝒟[:return])
end

# Build a REINFORCE solver
function REINFORCE(;π,
                    a_opt::NamedTuple=(;),
                    log::NamedTuple=(;),
                    kwargs...)

    OnPolicySolver(;agent=PolicyParams(π),
                    log = LoggerParams(;dir = "log/reinforce", log...),
                    a_opt = TrainingParams(;loss = reinforce_loss, early_stopping = (infos) -> (infos[end][:kl] > 0.015), name = "actor_", a_opt...),
                    kwargs...)
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE

using FEMBase
using Base.Test

using FEMBase: group_by_element_type
using FEMBase: get_local_coordinates, inside
using FEMBase: get_basis, get_dbasis, get_integration_order
using FEMBase: get_reference_element_coordinates, get_reference_coordinates
using FEMBase: get_element_type, is_element_type, get_element_id, filter_by_element_type

mutable struct Dummy <: FieldProblem
end

@testset "add elements to problem" begin
    problem = Problem(Dummy, "test", 2)
    element = Element(Quad4, [1, 2, 3, 4])
    elements = [element]
    add_elements!(problem, elements)
    @test problem.elements[1] == element
end

@testset "no unknown field name or assemble!-function defined" begin
    el = Element(Quad4, [1, 2, 3, 4])
    pr = Problem(Dummy, "problem", 2)
    add_elements!(pr, [el])
    assemble!(pr, 0.0)
    @test true
end

@testset "add time dependent field to element" begin
    el = Element(Seg2, [1, 2])
    u1 = Vector{Float64}[[0.0, 0.0], [0.0, 0.0]]
    u2 = Vector{Float64}[[1.0, 1.0], [1.0, 1.0]]
    update!(el, "displacement", 0.0 => u1)
    update!(el, "displacement", 1.0 => u2)
    @test length(el["displacement"]) == 2
    @test isapprox(el("displacement", [0.0], 0.0), [0.0, 0.0])
    @test isapprox(el("displacement", [0.0], 0.5), [0.5, 0.5])
    @test isapprox(el("displacement", [0.0], 1.0), [1.0, 1.0])
    el2 = Element(Poi1, [1])
    update!(el2, "force 1", 0.0 => 1.0)
end

@testset "add CVTV field to element" begin
    el = Element(Seg2, [1, 2])
    el["my field"] = (xi,t) -> xi[1]*t
    v = el("my field", [1.0], 2.0)
    @test isapprox(v, 2.0)
end

@testset "interpolate DCTI from element" begin
    el = Element(Seg2, [1, 2])
    update!(el, "foobar", 1.0)
    fb = el("foobar", [0.0], 0.0)
    @test isa(fb, Float64)
    @test isapprox(fb, 1.0)
end

@testset "add elements to elements" begin
    el1 = Element(Seg2, [1, 2])
    el2 = Element(Seg2, [3, 4])
    update!(el1, "master elements", [el2])
    lst = el1("master elements", 0.0)
    @test isa(lst, Vector)
end

@testset "extend basis" begin
    el = Element(Quad4, [1, 2, 3, 4])
    expected = [
        0.25 0.00 0.25 0.00 0.25 0.00 0.25 0.00
        0.00 0.25 0.00 0.25 0.00 0.25 0.00 0.25]
    @test isapprox(el([0.0, 0.0], 0.0, 2), expected)
end

@testset "group elements" begin
    e1 = Element(Seg2, [1, 2])
    e2 = Element(Quad4, [1, 2, 3, 4])
    elements = [e1, e2]
    r = group_by_element_type(elements)
    @test length(r) == 2
    @test first(r[Element{Seg2}]) == e1
    @test first(r[Element{Quad4}]) == e2
    @test get_element_type(e1) == Seg2
    e1.id = 1
    @test get_element_id(e1) == 1
    @test is_element_type(e1, Seg2)
    seg2_elements = filter_by_element_type(elements, Seg2)
end

@testset "inverse isoparametric mapping" begin
    el = Element(Quad4, [1, 2, 3, 4])
    X = Dict(
        1 => [0.0, 0.0],
        2 => [1.0, 0.0],
        3 => [1.0, 1.0],
        4 => [0.0, 1.0])
    update!(el, "geometry", X)
    time = 0.0
    X1 = el("geometry", [0.1, 0.2], time)
    xi = get_local_coordinates(el, X1, time)
    X2 = el("geometry", xi, time)
    info("X1 = $X1, X2 = $X2")
    @test isapprox(X1, X2)
end

@testset "inside of linear element" begin
    el = Element(Quad4, [1, 2, 3, 4])
    X = Dict{Int64, Vector{Float64}}(
        1 => [0.0, 0.0],
        2 => [1.0, 0.0],
        3 => [1.0, 1.0],
        4 => [0.0, 1.0])
    update!(el, "geometry", X)
    time = 0.0
    @test inside(el, [0.5, 0.5], time) == true
    @test inside(el, [1.0, 0.5], time) == true
    @test inside(el, [1.0, 1.0], time) == true
    @test inside(el, [1.01, 1.0], time) == false
    @test inside(el, [1.0, 1.01], time) == false
end

@testset "inside of quadratic element" begin
    el = Element(Tri6, [1, 2, 3, 4, 5, 6])
    X = Dict{Int64, Vector{Float64}}(
        1 => [0.0, 0.0],
        2 => [1.0, 0.0],
        3 => [0.0, 1.0],
        4 => [0.5, 0.2],
        5 => [0.8, 0.6],
        6 => [-0.2, 0.5])
    update!(el, "geometry", X)
    p = [0.94, 0.3] # visually checked to be inside
    @test inside(el, p, 0.0) == true
    p = [-0.2, 0.8] # visually checked to be outside
    @test inside(el, p, 0.0) == false
end

@testset "dict field" begin
    el = Element(Seg2, [1, 2])
    X = Dict(1 => [0.0, 0.0], 2 => [1.0, 0.0], 3 => [0.5, 0.5])
    f = field(X)
    #update!(el, "geometry", X)
    el["geometry"] = f
    @test isapprox(el("geometry")[1], [0.0, 0.0])
    @test isapprox(el("geometry", 0.0)[1], [0.0, 0.0])
    @test isapprox(el("geometry", [0.0], 0.0), [0.5, 0.0])
end

@testset "test nodal element" begin
    el = Element(Poi1, [1])
    @test get_basis(el, [0.0], 1.0) == [1]
    @test get_dbasis(el, [0.0], 1.0) == [0]
    @test isapprox(el([0.0], 0.0, Val{:detJ}), 1.0)
    @test get_integration_order(el.properties) == 1
    ips = get_integration_points(el.properties, 1)
    @test length(ips) == 1
    @test length(ips[1]) == 2
    @test size(el) == (0, 1)
    @test length(el) == 1
    @test get_reference_element_coordinates(Poi1) == Vector{Float64}[[0.0]]
end

@testset "get reference coordinates" begin
    el = Element(Seg2, [1, 2])
    X = get_reference_coordinates(el)
    @test isapprox(X[1][1], -1.0)
    @test isapprox(X[2][1],  1.0)
end

@testset "function field" begin
    function f(element, xi, time)
        return xi[1]*time
    end
    e1 = Element(Seg2, [1, 2])
    e1["f"] = f
    @test isapprox(e1("f", [1.0], 1.0), 1.0)
end

@testset "gradient, jacobian and determinant of jacobian" begin
    X = Dict(
             1 => [0.0, 0.0],
             2 => [2.0, 0.0],
             3 => [2.0, 2.0],
             4 => [0.0, 2.0])
    element = Element(Quad4, [1, 2, 3, 4])
    element.fields["geometry"] = field(X)
    J = element([0.0, 0.0], 0.0, Val{:Jacobian})
    detJ = element([0.0, 0.0], 0.0, Val{:detJ})
    @test isapprox(det(J), detJ)
    @test isapprox(detJ, 1.0)
    element2 = Element(Seg2, [1, 2])
    update!(element2, "geometry", X)
    detJ = element2([0.0], 0.0, Val{:detJ})
    @test isapprox(detJ, 1.0)
    X = Dict(
             1 => [0.0, 0.0, 0.0],
             2 => [2.0, 0.0, 0.0],
             3 => [2.0, 2.0, 0.0],
             4 => [0.0, 2.0, 0.0])
    element3 = Element(Quad4, [1, 2, 3, 4])
    update!(element3, "geometry", X)
    detJ = element3([0.0, 0.0], 0.0, Val{:detJ})
    @test isapprox(detJ, 1.0)
    grad = element([0.0, 0.0], 0.0, Val{:Grad})
    @test isapprox(grad, 1/4*[-1.0 1.0 1.0 -1.0; -1.0 -1.0 1.0 1.0])
    update!(element3, "youngs modulus", 1.0)
    @test isapprox(element3("youngs modulus", [0.0], 0.0), 1.0)
    update!(element3, "test", [1.0, 2.0, 3.0])
    test_val = element3("test", [0.0], 1.0)
end

@testset "update data to element" begin
    element = Element(Seg2, [1, 2])
    elements = Element[element]

#=
    # constant fields (not changing on element area)

    # time intependent

    # scalar + updating scalar
    update!(elements, "scalar 1", 0.0)
    @test isapprox(element("scalar 1", (0.0), 0.0), 0.0)
    update!(elements, "scalar 1", 1.0)
    @test isapprox(element("scalar 1", (0.0), 0.0), 1.0)

    # vector + updating vector
    update!(elements, "vector 1", [1.0, 2.0])
    @test isapprox(element("vector 1", (0.0), 0.0), [1.0, 2.0])
    update!(elements, "vector 1", [2.0, 1.0])
    @test isapprox(element("vector 1", (0.0), 0.0), [2.0, 1.0])

    # tensor + updating tensor
    update!(elements, "tensor 1", [1.0 2.0; 3.0 4.0])
    @test isapprox(element("tensor 1", [1.0 2.0; 3.0 4.0]))
    update!(elements, "tensor 1", [4.0 3.0; 2.0 1.0])
    @test isapprox(element("tensor 1", [4.0 3.0; 2.0 1.0]))

    # time dependent

    # scalar + interpolate
    update!(elements, "scalar 2", 0.0 => 0.0)
    @test isapprox(elements("scalar 2", (0.0), 0.5), 0.0)
    update!(elements, "scalar 2", 0.0 => 1.0)
    @test isapprox(elements("scalar 2", (0.0), 0.5), 1.0)
    update!(elements, "scalar 2", 1.0 => 0.0)
    @test isapprox(elements("scalar 2", (0.0), 0.5), 0.5)

    # vector + interpolate
    update!(elements, "vector 2", 0.0 => [0.0, 0.0])
    @test isapprox(elements("vector 2", (0.0), 0.5), [0.0, 0.0])
    update!(elements, "vector 2", 0.0 => [0.0, 1.0])
    @test isapprox(elements("vector 2", (0.0), 0.5), [0.0, 1.0])
    update!(elements, "vector 2", 1.0 => [1.0, 0.0])
    @test isapprox(elements("vector 2", (0.0), 0.5), [0.5, 0.5])

    # tensor + interpolate
    update!(elements, "tensor 2", 0.0 => [0.0 0.0; 0.0 0.0])
    @test isapprox(elements("tensor 2", (0.0), 0.5), [0.0 0.0; 0.0 0.0])
    update!(elements, "tensor 2", 0.0 => [1.0 1.0; 1.0 1.0])
    @test isapprox(elements("tensor 2", (0.0), 0.5), [1.0 1.0; 1.0 1.0])
    update!(elements, "tensor 2", 1.0 => [2.0 2.0; 2.0 2.0])
    @test isapprox(elements("tensor 2", (0.0), 0.5), [1.5 1.5; 1.5 1.5])

    # variable fields (changing on element area, needs interpolating)

    # time intependent

    # scalar changing on element area
    update!(elements, "scalar 3", (1.0, 2.0))
    @test isapprox(element("scalar 3", (0.0), 0.0), 1.5)
    update!(elements, "scalar 3", (2.0, 3.0))
    @test isapprox(element("scalar 3", (0.0), 0.0), 2.5)

    # vector changing on element area
    update!(elements, "vector 3", ([1.0, 2.0], [2.0, 3.0]))
    @test isapprox(element("vector 3", (0.0), 0.0), [1.5, 2.5])
    update!(elements, "vector 3", ([2.0, 3.0], [3.0, 4.0]))
    @test isapprox(element("vector 3", (0.0), 0.0), [2.5, 3.5])

    # tensor changing on element area
    update!(elements, "tensor 3", ([0.0 0.0; 0.0 0.0], [1.0 1.0; 1.0 1.0]))
    @test isapprox(element("tensor 3", (0.0), 0.0), [0.5 0.5; 0.5 0.5])
    update!(elements, "tensor 3", ([1.0 1.0; 1.0 1.0], [2.0 2.0; 2.0 2.0]))
    @test isapprox(element("tensor 3", (0.0), 0.0), [1.5 1.5; 1.5 1.5])

    # time dependent

    # scalar changing on element area over time
    update!(elements, "scalar 4", 0.0 => (1.0, 3.0))
    update!(elements, "scalar 4", 1.0 => (3.0, 5.0))
    @test isapprox(element("scalar 4", (0.0), 0.5), 3.0)

    # vector changing on element area over time
    update!(elements, "vector 4", 0.0 => ([1.0, 1.0], [3.0, 3.0]))
    update!(elements, "vector 4", 1.0 => ([3.0, 3.0], [5.0, 5.0]))
    @test isapprox(element("scalar 4", (0.0), 0.5), [3.0, 3.0])

    # tensor changing on element area over time
    update!(elements, "tensor 4", 0.0 => ([1.0 1.0; 1.0 1.0], [3.0 3.0; 3.0 3.0]))
    update!(elements, "tensor 4", 1.0 => ([3.0 3.0; 3.0 3.0], [5.0 5.0; 5.0 5.0]))
    @test isapprox(element("tensor 4", (0.0), 0.5), [3.0 3.0; 3.0 3.0])
=#

end

@testset "calculate displacement gradient" begin
    X = Dict(
             1 => [0.0, 0.0],
             2 => [1.0, 0.0],
             3 => [1.0, 1.0],
             4 => [0.0, 1.0])
    u = Dict(
             1 => [0.0, 0.0],
             2 => [1.0, -1.0],
             3 => [2.0, 3.0],
             4 => [0.0, 0.0])
    element = Element(Quad4, [1, 2, 3, 4])
    update!(element, "geometry", X)
    update!(element, "displacement", 0.0 => u)
    gradu = element("displacement", (0.0, 0.0), 0.0, Val{:Grad})
    gradu_expected = [1.5 0.5; 1.0 2.0]
    @test isapprox(gradu, gradu_expected)

    element = Element(Quad4, [1, 2, 3, 4])
    X = ([0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0])
    u = ([0.0,0.0], [1.0,-1.0], [2.0,3.0], [0.0,0.0])
    element["geometry"] = X
    element["displacement"] = u
    gradu = element("displacement", (0.0, 0.0), 0.0, Val{:Grad})
    gradu_expected = [1.5 0.5; 1.0 2.0]
    @test isapprox(gradu, gradu_expected)
    @test isapprox(element("geometry", (0.0,0.0), 0.0), [0.5,0.5])
end

@testset "interpolate field" begin
    element = Element(Quad4, [1,2,3,4])
    update!(element, "test", 0.0 => 0.0)
    update!(element, "test", 1.0 => 1.0)
    @test isapprox(element("test", (0.0, 0.0), 0.5), 0.5)
end

@testset "get integration points for Quad4" begin
    el = Element(Quad4, [1, 2, 3, 4])
    ips = get_integration_points(el)
    @test length(ips) == 4
end

@testset "test nonconverging inverse isoparametric mapping" begin
    X = Dict(1 => [0.0, 0.0], 2 => [1.0, 0.0], 3 => [1.0, 1.0], 4 => [0.0, 1.0])
    el = Element(Quad4, [1, 2, 3, 4])
    update!(el, "geometry", X)
    @test_throws Exception get_local_coordinates(el, [0.1, 0.1], 0.0; max_iterations=0)
end

@testset "analytical functions as fields" begin
    f(xi, time) = xi[1]*time
    g(element, xi, time) = element("geometry", xi, time)*time
    X = Dict(1 => [0.0, 0.0], 2 => [1.0, 0.0], 3 => [1.0, 1.0], 4 => [0.0, 1.0])
    el = Element(Quad4, [1, 2, 3, 4])
    update!(el, "geometry", X)
    update!(el, "f", f)
    update!(el, "g", g)
    @test isapprox(el("f", (0.5, 0.5), 2.0), 1.0)
    @test isapprox(el("g", (0.0, 0.0), 2.0), [1.0, 1.0])
end

@testset "add field to element using update!" begin
    el = Element(Seg2, [1, 2])
    f = field(1)
    update!(el, "f", f)
    @test f === el.fields["f"]
end

@testset "update DVTI/DVTV with Dict" begin
    el = Element(Seg2, [1, 2])
    update!(el, "u", 0.0 => (0.0, 0.0))
    update!(el, "v", (0.0, 0.0))
    u = Dict(1 => 1.0, 2 => 1.0)
    update!(el, "u", 0.0 => u)
    update!(el, "v", u)
    @test isapprox(el("u", (0.0,), 0.0), 1.0)
    @test isapprox(el("v", (0.0,), 0.0), 1.0)
end


text:
@testset "Common" begin
    @test get_code_center_frequency_ratio(GPSL1(use_gpu = Val(true))) ≈ 1/1540
end


text:
using jInv.Mesh
using Test

# setup 3D problem
nc = [5, 7, 8]
x0 = rand(3)
domain = [x0[1], 4, x0[2], 2, x0[3], 6]
h   = (domain[2:2:end]-domain[1:2:end])./nc
h1  = h[1]*ones(nc[1])
h2  = h[2]*ones(nc[2])
h3  = h[3]*ones(nc[3])

Mt = getTensorMesh3D(h1,h2,h3,x0)
Mr = getRegularMesh(domain,nc)
Mt2 = getTensorMesh3D(h1+rand(nc[1]),h2+rand(nc[2]),h3+rand(nc[3]),x0)

Meshes = (Mt,Mr,Mt2)

for k=1:length(Meshes)
	M = Meshes[k]
	print("\ttesting differential operators for $(typeof(M))...")
	@test norm(getDivergenceMatrix(M)*getCurlMatrix(M),Inf)<1e-12

	Lap = getNodalLaplacianMatrix(M)
	G   = getNodalGradientMatrix(M)
	@test norm(Lap*ones(size(Lap,2))/M.nc) < 1e-10
	@test norm(G'*G- Lap,Inf) <1e-12
	print("passed!\n")
end


text:
using Wavelets

wf = wavelet(WT.db4, WT.Filter)
wl = wavelet(WT.cdf97, WT.Lifting)
N = 1024;
x0 = rand(N,N);
L = maxtransformlevels(N)  # int(log2(N)-2)
tn = 20