text:
using ITensors,
      Test,
      Random

Random.seed!(1234)

@testset "BlockSparse ITensor" begin

  @testset "Constructor" begin
    i = Index([QN(0)=>1,QN(1)=>2],"i")
    j = Index([QN(0)=>3,QN(1)=>4,QN(2)=>5],"j")

    A = ITensor(QN(0),i,dag(j))

    @test flux(A) == QN(0)
    @test nnzblocks(A) == 2
  end

  @testset "Constructor (from Tuple)" begin
    i = Index([QN(0)=>1,QN(1)=>2],"i")
    j = Index([QN(0)=>3,QN(1)=>4,QN(2)=>5],"j")

    A = ITensor(QN(0), (i, dag(j)))

    @test flux(A) == QN(0)
    @test nnzblocks(A) == 2
  end

  @testset "Constructor (no flux specified)" begin
    i = Index([QN(0)=>1, QN(1)=>2], "i")
    j = Index([QN(0)=>3, QN(1)=>4, QN(2)=>5], "j")

    A = ITensor(i, dag(j))

    @test flux(A) == QN(0)
    @test nnzblocks(A) == 2
  end

  @testset "Constructor (Tuple, no flux specified)" begin
    i = Index([QN(0)=>1, QN(1)=>2], "i")
    j = Index([QN(0)=>3, QN(1)=>4, QN(2)=>5], "j")

    A = ITensor((i, dag(j)))

    @test flux(A) == QN(0)
    @test nnzblocks(A) == 2
  end

  @testset "No indices getindex" begin
    T = ITensor(QN())
    @test order(T) == 0
    @test flux(T) == nothing
    @test nnzblocks(T) == 1
    @test T[] == 0

    s = Index(QN(-1)=>1,QN(1)=>1)
    A = emptyITensor(s, dag(s'))
    B = emptyITensor(s', dag(s))
    A[1, 1] = 1
    B[2, 2] = 1
    C = A * B
    @test order(C) == 0
    @test flux(C) == nothing
    @test nnzblocks(C) == 0
    @test C[] == 0
  end

  @testset "Empty constructor" begin
    i = Index([QN(0)=>1,QN(1)=>2],"i")

    A = emptyITensor(i,dag(i'))

    @test nnzblocks(A) == 0
    @test nnz(A) == 0
    @test hasinds(A,i,i')
    @test isnothing(flux(A))

    A[i(1),i'(1)] = 1.0

    @test nnzblocks(A) == 1
    @test nnz(A) == 1
    @test flux(A) == QN(0)

    A[i(2),i'(2)] = 1.0

    @test nnzblocks(A) == 2
    @test nnz(A) == 5
    @test flux(A) == QN(0)
  end

  @testset "Check flux when setting elements" begin
    i = Index(QN(0)=>1,QN(1)=>1; tags="i")
    A = randomITensor(QN(0),i,dag(i'))
    @test_throws ErrorException A[i(1),i'(2)] = 1.0
  end

  @testset "Random constructor" begin
    i = Index([QN(0)=>1,QN(1)=>2],"i")
    j = Index([QN(0)=>3,QN(1)=>4,QN(2)=>5],"j")

    A = randomITensor(QN(1),i,dag(j))

    @test flux(A) == QN(1)
    @test nnzblocks(A) == 1

    B = randomITensor(i,dag(j))

    @test flux(B) == QN()
    @test nnzblocks(B) == 2
  end

  @testset "QN setelt" begin
    i = Index(QN(0)=>2,QN(1)=>2,tags="i")

    T = setelt(i(1))
    @test T[i(1)] ≈ 1.0
    @test T[i(2)] ≈ 0.0
    @test T[i(3)] ≈ 0.0
    @test T[i(4)] ≈ 0.0

    T = setelt(i(2))
    @test T[i(1)] ≈ 0.0
    @test T[i(2)] ≈ 1.0
    @test T[i(3)] ≈ 0.0
    @test T[i(4)] ≈ 0.0

    # Test setelt taking Pair{Index,Int}
    T = setelt(i=>3)
    @test T[i(1)] ≈ 0.0
    @test T[i(2)] ≈ 0.0
    @test T[i(3)] ≈ 1.0
    @test T[i(4)] ≈ 0.0
  end

  @testset "setindex!" begin

    @testset "Test 1" 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")
      A = emptyITensor(s1,s2)

      @test nnzblocks(A) == 0
      @test nnz(A) == 0
      @test hasinds(A,s1,s2)
      @test isnothing(flux(A))

      A[2,1] = 1.0/sqrt(2)

      @test nnzblocks(A) == 1
      @test nnz(A) == 1
      @test A[s1(2),s2(1)] ≈ 1.0/sqrt(2)
      @test flux(A) == QN("N",1,-1)

      A[1,2] = 1.0/sqrt(2)

      @test nnzblocks(A) == 2
      @test nnz(A) == 2
      @test A[s1(2),s2(1)] ≈ 1.0/sqrt(2)
      @test A[s1(1),s2(2)] ≈ 1.0/sqrt(2)
      @test flux(A) == QN("N",1,-1)
    end

    @testset "Test 2" 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")
      A = emptyITensor(s1,s2)

      @test nnzblocks(A) == 0
      @test nnz(A) == 0
      @test hasinds(A,s1,s2)
      @test isnothing(flux(A))

      A[1,2] = 1.0/sqrt(2)

      @test nnzblocks(A) == 1
      @test nnz(A) == 1
      @test A[s1(1),s2(2)] ≈ 1.0/sqrt(2)
      @test flux(A) == QN("N",1,-1)

      A[2,1] = 1.0/sqrt(2)

      @test nnzblocks(A) == 2
      @test nnz(A) == 2
      @test A[s1(2),s2(1)] ≈ 1.0/sqrt(2)
      @test A[s1(1),s2(2)] ≈ 1.0/sqrt(2)
      @test flux(A) == QN("N",1,-1)
    end

  end

  @testset "Multiply by scalar" begin
    i = Index([QN(0)=>1,QN(1)=>2],"i")
    j = Index([QN(0)=>3,QN(1)=>4],"j")

    A = randomITensor(QN(0),i,dag(j))

    @test flux(A) == QN(0)
    @test nnzblocks(A) == 2

    B = 2*A

    @test flux(B) == QN(0)
    @test nnzblocks(B) == 2

    for ii in dim(i), jj in dim(j)
      @test 2*A[i=>ii,j=>jj] == B[i=>ii,j=>jj]
    end
  end

  @testset "Copy" begin
    s = Index([QN(0)=>1,QN(1)=>1],"s")
    T = randomITensor(QN(0),s,s')
    cT = copy(T)
    for ss in dim(s), ssp in dim(s')
      @test T[s=>ss,s'=>ssp] == cT[s=>ss,s'=>ssp]
    end
  end

  @testset "Permute" begin
    i = Index([QN(0)=>1,QN(1)=>2],"i")
    j = Index([QN(0)=>3,QN(1)=>4],"j")

    A = randomITensor(QN(1),i,dag(j))

    @test flux(A) == QN(1)
    @test nnzblocks(A) == 1

    B = permute(A,j,i)

    @test flux(B) == QN(1)
    @test nnzblocks(B) == 1

    for ii in dim(i), jj in dim(j)
      @test A[ii,jj] == B[jj,ii]
    end
  end

  @testset "Contraction" begin
    i = Index([QN(0)=>1,QN(1)=>2],"i")
    j = Index([QN(0)=>3,QN(1)=>4],"j")

    A = randomITensor(QN(0),i,dag(j))

    @test flux(A) == QN(0)
    @test nnzblocks(A) == 2

    B = randomITensor(QN(1),j,dag(i)')

    @test flux(B) == QN(1)
    @test nnzblocks(B) == 1

    C = A*B

    @test hasinds(C,i,i')
    @test flux(C) == QN(1)
    @test nnzblocks(C) == 1
  end

  @testset "Combine and uncombine" begin

    @testset "Combine no indices" begin
      i1 = Index([QN(0,2)=>2,QN(1,2)=>2],"i1")
      A = randomITensor(QN(),i1,dag(i1'))

      C = combiner()
      c = combinedind(C)
      @test isnothing(c)
      AC = A*C
      @test nnz(AC) == nnz(A)
      @test nnzblocks(AC) == nnzblocks(A)
      @test hassameinds(AC,A)
      @test norm(AC-A*C) ≈ 0.0
      Ap = AC*dag(C)
      @test nnz(Ap) == nnz(A)
      @test nnzblocks(Ap) == nnzblocks(A)
      @test hassameinds(Ap,A)
      @test norm(A-Ap) ≈ 0.0
    end

    @testset "Combine set direction" begin
      i1 = Index([QN(0)=>2,QN(1)=>3],"i1")
      A = randomITensor(i1',dag(i1))
      @test isnothing(ITensors.checkflux(A))
      C = combiner(dag(i1); dir = ITensors.Out)
      c = combinedind(C)
      @test dir(c) == ITensors.Out
      AC = A * C
      @test nnz(AC) == nnz(A)
      @test nnzblocks(AC) == nnzblocks(A)
      @test isnothing(ITensors.checkflux(AC))
      Ap = AC*dag(C)
      @test nnz(Ap) == nnz(A)
      @test nnzblocks(Ap) == nnzblocks(A)
      @test hassameinds(Ap, A)
      @test isnothing(ITensors.checkflux(AC))
      @test A ≈ Ap
    end

    @testset "Order 2 (IndexSet constructor)" begin
      i1 = Index([QN(0,2)=>2,QN(1,2)=>2],"i1")
      A = randomITensor(QN(),i1,dag(i1'))

      iss = [i1,
             dag(i1'),
             (i1,dag(i1')),
             (dag(i1'),i1)]

      for is in iss
        C = combiner(IndexSet(is); tags="c")
        AC = A*C
        @test nnz(AC) == nnz(A)
        Ap = AC*dag(C)
        @test nnz(Ap) == nnz(A)
        @test nnzblocks(Ap) == nnzblocks(A)
        @test norm(A-Ap) ≈ 0.0
      end
    end

    @testset "Order 2" begin
      i1 = Index([QN(0,2)=>2,QN(1,2)=>2],"i1")
      A = randomITensor(QN(), (i1, dag(i1')))

      iss = [i1,
             dag(i1'),
             (i1,dag(i1')),
             (dag(i1'),i1)]

      for is in iss
        C = combiner(is; tags="c")
        AC = A*C
        @test nnz(AC) == nnz(A)
        Ap = AC*dag(C)
        @test nnz(Ap) == nnz(A)
        @test nnzblocks(Ap) == nnzblocks(A)
        @test norm(A-Ap) ≈ 0.0
      end
    end

    @testset "Order 3, Combine 2" begin
      i = Index([QN(0)=>2,QN(1)=>2],"i")

      A = randomITensor(QN(0),i,dag(i)',dag(i)'')

      C = combiner(i,dag(i)'')
      c = combinedind(C)

      AC = A*C

      @test hasinds(AC,c,i')
      @test nnz(AC) == nnz(A)

      for b in nzblocks(AC)
        @test flux(AC,b) == QN(0)
      end

      B = ITensor(QN(0),i',c)
      @test nnz(B) == nnz(AC)
      @test nnzblocks(B) == nnzblocks(AC)

      Ap = AC*dag(C)

      @test norm(A-Ap) == 0
      @test nnz(A) == nnz(Ap)
      @test nnzblocks(A) == nnzblocks(Ap)
      @test hassameinds(A,Ap)
    end

    @testset "Order 3" begin
      i1 = Index([QN(0,2)=>2,QN(1,2)=>2],"i1")
      i2 = settags(i1,"i2")
      A = randomITensor(QN(),i1,i2,dag(i1'))

      iss = [i1,
             i2,
             dag(i1'),
             (i1,i2),
             (i2,i1),
             (i1,dag(i1')),
             (dag(i1'),i1),
             (i2,dag(i1')),
             (dag(i1'),i2),
             (i1,i2,dag(i1')),
             (i1,dag(i1'),i2),
             (i2,i1,dag(i1')),
             (i2,dag(i1'),i1),
             (dag(i1'),i1,i2),
             (dag(i1'),i2,i1)]

      for is in iss
        C = combiner(is; tags="c")
        AC = A*C
        @test nnz(AC) == nnz(A)
        Ap = AC*dag(C)
        @test nnz(Ap) == nnz(A)
        @test nnzblocks(Ap) == nnzblocks(A)
        @test norm(A-AC*dag(C)) ≈ 0.0
      end
    end

    @testset "Order 4" begin
      i1 = Index([QN(0,2)=>2,QN(1,2)=>2],"i1")
      i2 = settags(i1,"i2")
      A = randomITensor(QN(),i1,i2,dag(i1'),dag(i2'))

      iss = [i1,
             i2,
             dag(i1'),
             dag(i2'),
             (i1,i2),
             (i2,i1),
             (i1,dag(i1')),
             (dag(i1'),i1),
             (i1,dag(i2')),
             (dag(i2'),i1),
             (i2,dag(i1')),
             (dag(i1'),i2),
             (i2,dag(i2')),
             (dag(i2'),i2),
             (dag(i1'),dag(i2')),
             (dag(i2'),dag(i1')),
             (i1,i2,dag(i1')),
             (i1,dag(i1'),i2),
             (i2,i1,dag(i1')),
             (i2,dag(i1'),i1),
             (dag(i1'),i1,i2),
             (dag(i1'),i2,i1),
             (i1,dag(i1'),dag(i2')),
             (i1,dag(i2'),dag(i1')),
             (dag(i1'),i1,dag(i2')),
             (dag(i1'),dag(i2'),i1),
             (dag(i2'),i1,dag(i1')),
             (dag(i2'),dag(i1'),i1),
             (i1,i2,dag(i1'),dag(i2')),
             (i1,i2,dag(i2'),dag(i1')),
             (i1,dag(i1'),i2,dag(i2')),
             (i1,dag(i1'),dag(i2'),i2),
             (i1,dag(i2'),i2,dag(i1')),
             (i1,dag(i2'),dag(i1'),i2),
             (i2,i1,dag(i1'),dag(i2')),
             (i2,i1,dag(i2'),dag(i1')),
             (i2,dag(i1'),i1,dag(i2')),
             (i2,dag(i1'),dag(i2'),i1),
             (i2,dag(i2'),i1,dag(i1')),
             (i2,dag(i2'),dag(i1'),i1),
             (dag(i1'),i2,i1,dag(i2')),
             (dag(i1'),i2,dag(i2'),i1),
             (dag(i1'),i1,i2,dag(i2')),
             (dag(i1'),i1,dag(i2'),i2),
             (dag(i1'),dag(i2'),i2,i1),
             (dag(i1'),dag(i2'),i1,i2),
             (dag(i2'),i1,dag(i1'),i2),
             (dag(i2'),i1,i2,dag(i1')),
             (dag(i2'),dag(i1'),i1,i2),
             (dag(i2'),dag(i1'),i2,i1),
             (dag(i2'),i2,i1,dag(i1')),
             (dag(i2'),i2,dag(i1'),i1)]

      for is in iss
        C = combiner(is; tags="c")
        AC = A*C
        @test nnz(AC) == nnz(A)
        Ap = AC*dag(C)
        @test nnz(Ap) == nnz(A)
        @test nnzblocks(Ap) == nnzblocks(A)
        @test norm(A-Ap) ≈ 0.0
      end
    end

    @testset "Order 4, Combine 2, Example 1" begin
      s1 = Index([QN(("Sz", 0),("Nf",0))=>1,
                  QN(("Sz",+1),("Nf",1))=>1,
                  QN(("Sz",-1),("Nf",1))=>1,
                  QN(("Sz", 0),("Nf",2))=>1],"site,n=1");
      s2 = replacetags(s1,"n=1","n=2")

      A = randomITensor(QN(),s1,s2,dag(s1)',dag(s2)')

      C = combiner(dag(s1)',dag(s2)')
      c = combinedind(C)

      AC = A*C

      @test norm(AC) ≈ norm(A)
      @test hasinds(AC,s1,s2,c)
      @test nnz(AC) == nnz(A)
      for b in nzblocks(AC)
        @test flux(AC,b) == QN()
      end

      @test nnzblocks(AC) < nnz(A)

      B = ITensor(QN(),s1,s2,c)
      @test nnz(B) == nnz(AC)
      @test nnzblocks(B) == nnzblocks(AC)

      Ap = AC*dag(C)

      @test hassameinds(A,Ap)
      @test norm(A-Ap) == 0
      @test nnz(A) == nnz(Ap)
      @test nnzblocks(A) == nnzblocks(Ap)
    end

    @testset "Order 4, Combine 2, Example 2" begin
      s1 = Index([QN(("Nf",0))=>1,
                  QN(("Nf",1))=>1],"site,n=1")
      s2 = replacetags(s1,"n=1","n=2")

      A = randomITensor(QN(),dag(s2)',s2,dag(s1)',s1)

      C = combiner(dag(s2)',dag(s1)')
      c = combinedind(C)

      AC = A*C

      @test norm(AC) ≈ norm(A)
      @test hasinds(AC,s1,s2,c)
      @test nnz(AC) == nnz(A)
      for b in nzblocks(AC)
        @test flux(AC,b) == QN()
      end

      B = ITensor(QN(),s1,s2,c)
      @test nnzblocks(B) == nnzblocks(AC)

      Ap = AC*dag(C)

      @test hassameinds(A,Ap)
      @test norm(A-Ap) == 0
      @test nnz(A) == nnz(Ap)
      @test nnzblocks(A) == nnzblocks(Ap)
    end

    @testset "Order 4, Combine 2, Example 3" begin
      s1 = Index([QN(("Nf",0))=>1,
                  QN(("Nf",1))=>1],"site,n=1")
      s2 = replacetags(s1,"n=1","n=2")

      A = randomITensor(QN(),dag(s1)',s2,dag(s2)',s1)

      C = combiner(dag(s2)',dag(s1)')
      c = combinedind(C)

      AC = A*C

      @test norm(AC) ≈ norm(A)
      @test hasinds(AC,s1,s2,c)
      @test nnz(AC) == nnz(A)
      for b in nzblocks(AC)
        @test flux(AC,b) == QN()
      end

      B = ITensor(QN(),s1,s2,c)
      @test nnzblocks(B) == nnzblocks(AC)

      Ap = AC*dag(C)

      @test hassameinds(A,Ap)
      @test norm(A-Ap) == 0
      @test nnz(A) == nnz(Ap)
      @test nnzblocks(A) == nnzblocks(Ap)
    end
  end

  @testset "Check that combiner commutes" begin
    i = Index(QN(0,2)=>2,QN(1,2)=>2; tags="i")
    j = settags(i,"j")
    A = randomITensor(QN(0,2),i,j,dag(i'),dag(j'))
    C = combiner(i,j)
    @test norm(A*dag(C')*C-A*C*dag(C')) ≈ 0.0
  end

  @testset "Combiner for block deficient ITensor" begin
    i = Index(QN(0,2)=>2,QN(1,2)=>2; tags="i")
    j = settags(i,"j")
    A = emptyITensor(i,j,dag(i'))
    A[1,1,1] = 1.0
    C = combiner(i,j; tags="c")
    AC = A*C
    Ap = AC*dag(C)
    @test norm(A-Ap) ≈ 0.0
    @test norm(Ap-A) ≈ 0.0
  end

  @testset "Combine Complex ITensor" begin
    s1 = Index(QN(("Sz",1)) => 1,QN(("Sz",-1)) => 1;tags="S=1/2,Site,n=1")
    s2 = Index(QN(("Sz",1)) => 1,QN(("Sz",-1)) => 1;tags="S=1/2,Site,n=2")

    T = randomITensor(ComplexF64,QN("Sz",0),s1,s2)

    C = combiner(s1,s2)
    CT = C*T
    @test norm(CT) ≈ norm(T)
    TT = dag(C)*CT
    @test TT ≈ T
  end

  @testset "Combiner bug #395" begin
    i1 = Index([QN(0)=>1,QN(1)=>2],"i1")
    i2 = Index([QN(0)=>1,QN(1)=>2],"i2")
    A = randomITensor(QN(),i1,i2,dag(i1)', dag(i2)')
    CL = combiner(i1,i2)
    CR = combiner(dag(i1)',dag(i2)')
    AC = A * CR * CL
    @test AC * dag(CR) * dag(CL) ≈ A
  end

  @testset "Contract to scalar" begin
    i = Index([QN(0)=>1,QN(1)=>1],"i")
    A = randomITensor(QN(0),i,dag(i'))

    c = A*dag(A)

    @test nnz(c) == 1
    @test nnzblocks(c) == 1
    @test c[] isa Float64
    @test c[] ≈ norm(A)^2
  end

  @testset "eigen" begin

    @testset "eigen hermitian" begin
      i = Index(QN(0)=>2,QN(1)=>3,QN(2)=>4; tags="i")
      j = settags(i,"j")
      k = settags(i,"k")
      l = settags(i,"l")

      A = randomITensor(QN(),i,j,dag(k),dag(l))
      A = A*prime(dag(A),(i,j))

      F = eigen(A; ishermitian=true, tags="x")

      D, U = F
      Ut = F.Vt

      @test store(U) isa NDTensors.BlockSparse
      @test store(D) isa NDTensors.DiagBlockSparse

      u = commonind(D,U)
      up = uniqueind(D,U)

      @test hastags(u,"x")
      @test plev(u) == 0
      @test hastags(up,"x")
      @test plev(up) == 1

      @test hassameinds(U, (i,j,u))
      @test hassameinds(D, (u,up))

      @test norm(A - dag(U) * D * U') ≈ 0.0 atol=1e-11
      @test norm(A - dag(U) * D * Ut) ≈ 0.0 atol=1e-11
      @test norm(A * U - U' * D) ≈ 0.0 atol=1e-11
      @test norm(A * U - Ut * D) ≈ 0.0 atol=1e-11
    end

    @testset "eigen hermitian (truncate)" begin
      i = Index(QN(0)=>2,QN(1)=>3,QN(2)=>4; tags="i")
      j = settags(i,"j")
      k = settags(i,"k")
      l = settags(i,"l")

      A = randomITensor(QN(),i,j,dag(k),dag(l))
      A = A*prime(dag(A),(i,j))
      for i = 1:4
        A = mapprime(A*A',2,1)
      end
      A = A/norm(A)

      cutoff = 1e-5
      F = eigen(A; ishermitian=true,
                   tags="x",
                   cutoff=cutoff)

      D, U, spec = F
      Ut = F.Vt

      @test store(U) isa NDTensors.BlockSparse
      @test store(D) isa NDTensors.DiagBlockSparse

      u = commonind(D, U)
      up = uniqueind(D, U)

      @test hastags(u,"x")
      @test plev(u) == 0
      @test hastags(up,"x")
      @test plev(up) == 1

      @test hassameinds(U,(i,j,u))
      @test hassameinds(D,(u,up))

      for b in nzblocks(A)
        @test flux(A,b)==QN(0)
      end
      for b in nzblocks(U)
        @test flux(U,b)==QN(0)
      end
      for b in nzblocks(D)
        @test flux(D,b)==QN(0)
      end

      Ap = dag(U) * D * U'

      @test norm(Ap-A) ≤ 1e-2
      @test norm(dag(U) * D * Ut - A) ≤ 1e-2
      @test minimum(dims(D)) == length(spec.eigs)
      @test minimum(dims(D)) < dim(i)*dim(j)

      @test spec.truncerr ≤ cutoff
      err = sqrt(1-(Ap*dag(Ap))[]/(A*dag(A))[])
      @test err ≤ cutoff
      @test err ≈ spec.truncerr rtol=1e-1
		end

    @testset "eigen non-hermitian" begin
      i = Index(QN(0)=>2,QN(1)=>3,QN(2)=>4; tags="i")
      j = settags(i,"j")

      A = randomITensor(QN(),i,j,dag(i'),dag(j'))

      F = eigen(A; tags="x")

      D, U = F
      Ut = F.Vt

      @test store(U) isa NDTensors.BlockSparse
      @test store(D) isa NDTensors.DiagBlockSparse

      u = commonind(D,U)
      up = uniqueind(D,U)

      @test hastags(u,"x")
      @test plev(u) == 0
      @test hastags(up,"x")
      @test plev(up) == 1

      @test norm(A - U' * D * dag(U)) ≉ 0.0 atol=1e-12
      @test norm(A - Ut * D * dag(U)) ≉ 0.0 atol=1e-12
      @test norm(A * U - U' * D) ≈ 0.0 atol=1e-12
      @test norm(A * U - Ut * D) ≈ 0.0 atol=1e-12
    end

    @testset "eigen non-hermitian (general inds)" begin
      i = Index(QN(0) => 2,
                QN(1) => 3,
                QN(2) => 4; tags="i")
      j = settags(i, "j")
      ĩ, j̃ = sim(i), sim(j)

      A = randomITensor(QN(), i, j, dag(ĩ),dag(j̃))

      F = eigen(A, (i, j), (ĩ, j̃); lefttags = "x", righttags = "y")

      D, U = F
      Ut = F.Vt

      @test store(U) isa NDTensors.BlockSparse
      @test store(D) isa NDTensors.DiagBlockSparse

      l = uniqueind(D, U)
      r = commonind(D, U)

      @test F.l == l
      @test F.r == r

      @test hastags(l, "x")
      @test plev(l) == 0
      @test hastags(r, "y")
      @test plev(r) == 0

      @test hassameinds(U, (ĩ, j̃, r))
      @test hassameinds(Ut, (i, j, l))

      @test norm(A * U - Ut * D) ≈ 0.0 atol=1e-12
      @test norm(A - Ut * D * dag(U)) ≉ 0.0 atol=1e-12
    end

    @testset "eigen mixed arrows" begin
      i1 = Index([QN(0)=>1,QN(1)=>2],"i1")
      i2 = Index([QN(0)=>1,QN(1)=>2],"i2")
      A = randomITensor(i1, i2, dag(i1)', dag(i2)')
      F = eigen(A, (i1, i1'), (i2', i2))
      D, U = F
      Ut = F.Vt
      @test norm(A * U - Ut * D) ≈ 0.0 atol=1e-12
    end

  end

  @testset "svd" for ElT ∈ (Float64, ComplexF64)

    @testset "svd example 1" begin
      i = Index(QN(0)=>2,QN(1)=>2; tags="i")
      j = Index(QN(0)=>2,QN(1)=>2; tags="j")
      A = randomITensor(ElT, QN(0),i,dag(j))
      for b in nzblocks(A)
        @test flux(A,b)==QN(0)
      end
      U,S,V = svd(A,i)

      @test store(U) isa NDTensors.BlockSparse
      @test store(S) isa NDTensors.DiagBlockSparse
      @test store(V) isa NDTensors.BlockSparse

      for b in nzblocks(U)
        @test flux(U,b)==QN(0)
      end
      for b in nzblocks(S)
        @test flux(S,b)==QN(0)
      end
      for b in nzblocks(V)
        @test flux(V,b)==QN(0)
      end
      @test isapprox(norm(U*S*V-A),0.0; atol=1e-14)
    end

    @testset "svd example 2" begin
      i = Index(QN(0)=>5,QN(1)=>6; tags="i")
      j = Index(QN(-1)=>2,QN(0)=>3,QN(1)=>4; tags="j")
      A = randomITensor(ElT, QN(0),i,j)
      for b in nzblocks(A)
        @test flux(A,b)==QN(0)
      end
      U,S,V = svd(A,i)

      @test store(U) isa NDTensors.BlockSparse
      @test store(S) isa NDTensors.DiagBlockSparse
      @test store(V) isa NDTensors.BlockSparse

      for b in nzblocks(U)
        @test flux(U,b)==QN(0)
      end
      for b in nzblocks(S)
        @test flux(S,b)==QN(0)
      end
      for b in nzblocks(V)
        @test flux(V,b)==QN(0)
      end
      @test isapprox(norm(U*S*V-A),0.0; atol=1e-14)
    end

    @testset "svd example 3" begin
      i = Index(QN(0)=>5,QN(1)=>6; tags="i")
      j = Index(QN(-1)=>2,QN(0)=>3,QN(1)=>4; tags="j")
      A = randomITensor(ElT, QN(0),i,dag(j))
      for b in nzblocks(A)
        @test flux(A,b)==QN(0)
      end
      U,S,V = svd(A,i)

      @test store(U) isa NDTensors.BlockSparse
      @test store(S) isa NDTensors.DiagBlockSparse
      @test store(V) isa NDTensors.BlockSparse

      for b in nzblocks(U)
        @test flux(U,b)==QN(0)
      end
      for b in nzblocks(S)
        @test flux(S,b)==QN(0)
      end
      for b in nzblocks(V)
        @test flux(V,b)==QN(0)
      end
      @test isapprox(norm(U*S*V-A),0.0; atol=1e-14)
    end

    @testset "svd example 4" begin
			i = Index(QN(0,2)=>2,QN(1,2)=>2; tags="i")
			j = settags(i,"j")

			A = randomITensor(ElT, QN(0,2),i,j,dag(i'),dag(j'))

			U,S,V = svd(A,i,j)

      @test store(U) isa NDTensors.BlockSparse
      @test store(S) isa NDTensors.DiagBlockSparse
      @test store(V) isa NDTensors.BlockSparse

      for b in nzblocks(A)
        @test flux(A,b)==QN(0,2)
      end
      U,S,V = svd(A,i)
      for b in nzblocks(U)
        @test flux(U,b)==QN(0,2)
      end
      for b in nzblocks(S)
        @test flux(S,b)==QN(0,2)
      end
      for b in nzblocks(V)
        @test flux(V,b)==QN(0,2)
      end
      @test isapprox(norm(U*S*V-A),0.0; atol=1e-14)
    end

    @testset "svd example 5" begin
			i = Index(QN(0,2)=>2,QN(1,2)=>2; tags="i")
			j = settags(i,"j")

			A = randomITensor(ElT, QN(1,2),i,j,dag(i'),dag(j'))

			U,S,V = svd(A,i,j)

      @test store(U) isa NDTensors.BlockSparse
      @test store(S) isa NDTensors.DiagBlockSparse
      @test store(V) isa NDTensors.BlockSparse

      for b in nzblocks(A)
        @test flux(A,b)==QN(1,2)
      end
      U,S,V = svd(A,i)
      for b in nzblocks(U)
        @test flux(U,b)==QN(0,2)
      end
      for b in nzblocks(S)
        @test flux(S,b)==QN(1,2)
      end
      for b in nzblocks(V)
        @test flux(V,b)==QN(0,2)
      end
      @test isapprox(norm(U*S*V-A),0.0; atol=1e-14)
    end

    @testset "svd example 6" begin
			i = Index(QN(0,2)=>2,QN(1,2)=>2; tags="i")
			j = settags(i,"j")

			A = randomITensor(ElT, QN(1,2),i,j,dag(i'),dag(j'))

			U,S,V = svd(A,i,i')

      @test store(U) isa NDTensors.BlockSparse
      @test store(S) isa NDTensors.DiagBlockSparse
      @test store(V) isa NDTensors.BlockSparse

      for b in nzblocks(A)
        @test flux(A,b)==QN(1,2)
      end
      U,S,V = svd(A,i)
      for b in nzblocks(U)
        @test flux(U,b)==QN(0,2)
      end
      for b in nzblocks(S)
        @test flux(S,b)==QN(1,2)
      end
      for b in nzblocks(V)
        @test flux(V,b)==QN(0,2)
      end
      @test isapprox(norm(U*S*V-A),0.0; atol=1e-14)
    end

    @testset "svd truncation example 1" begin
      i = Index(QN(0)=>2,QN(1)=>3; tags="i")
      j = settags(i,"j")
      A = randomITensor(ElT, QN(0),i,j,dag(i'),dag(j'))
      for i = 1:4
        A = mapprime(A*A',2,1)
      end
      A = A/norm(A)

      cutoff = 1e-5
      U,S,V,spec = svd(A,i,j; utags="x", vtags="y", cutoff=cutoff)

      @test store(U) isa NDTensors.BlockSparse
      @test store(S) isa NDTensors.DiagBlockSparse
      @test store(V) isa NDTensors.BlockSparse

      u = commonind(S,U)
      v = commonind(S,V)

      @test hastags(u,"x")
      @test hastags(v,"y")

      @test hassameinds(U,(i,j,u))
      @test hassameinds(V,(i',j',v))

      for b in nzblocks(A)
        @test flux(A,b)==QN(0)
      end
      for b in nzblocks(U)
        @test flux(U,b)==QN(0)
      end
      for b in nzblocks(S)
        @test flux(S,b)==QN(0)
      end
      for b in nzblocks(V)
        @test flux(V,b)==QN(0)
      end

      Ap = U*S*V

      @test norm(Ap-A) ≤ 1e-2
      @test minimum(dims(S)) == length(spec.eigs)
      @test minimum(dims(S)) < dim(i)*dim(j)

      @test spec.truncerr ≤ cutoff
      err = real(1-(Ap*dag(Ap))[]/(A*dag(A))[])
      @test err ≤ cutoff
      @test isapprox(err,spec.truncerr; rtol=1e-6)
    end

    @testset "svd truncation example 2" begin
      i = Index(QN(0)=>3,QN(1)=>2; tags="i")
      j = settags(i,"j")
      A = randomITensor(ElT, QN(0),i,j,dag(i'),dag(j'))

      maxdim = 4
      U,S,V,spec = svd(A,i,j; utags="x", vtags="y", maxdim=maxdim)

      @test store(U) isa NDTensors.BlockSparse
      @test store(S) isa NDTensors.DiagBlockSparse
      @test store(V) isa NDTensors.BlockSparse

      u = commonind(S,U)
      v = commonind(S,V)

      @test hastags(u,"x")
      @test hastags(v,"y")

      @test hassameinds(U,(i,j,u))
      @test hassameinds(V,(i',j',v))

      for b in nzblocks(A)
        @test flux(A,b)==QN(0)
      end
      for b in nzblocks(U)
        @test flux(U,b)==QN(0)
      end
      for b in nzblocks(S)
        @test flux(S,b)==QN(0)
      end
      for b in nzblocks(V)
        @test flux(V,b)==QN(0)
      end

      @test minimum(dims(S)) == maxdim
      @test minimum(dims(S)) == length(spec.eigs)
      @test minimum(dims(S)) < dim(i)*dim(j)

      Ap = U*S*V
      err = 1-(Ap*dag(Ap))[]/(A*dag(A))[]
      @test isapprox(err,spec.truncerr; rtol=1e-6)
    end

    @testset "svd truncation example 3" begin
      i = Index(QN(0)=>2,QN(1)=>3,QN(2)=>4; tags="i")
      j = settags(i,"j")
      A = randomITensor(ElT, QN(1),i,j,dag(i'),dag(j'))

      maxdim = 4
      U,S,V,spec = svd(A,i,j; utags="x", vtags="y", maxdim=maxdim)

      @test store(U) isa NDTensors.BlockSparse
      @test store(S) isa NDTensors.DiagBlockSparse
      @test store(V) isa NDTensors.BlockSparse

      u = commonind(S,U)
      v = commonind(S,V)

      @test hastags(u,"x")
      @test hastags(v,"y")

      @test hassameinds(U,(i,j,u))
      @test hassameinds(V,(i',j',v))

      for b in nzblocks(A)
        @test flux(A,b)==QN(1)
      end
      for b in nzblocks(U)
        @test flux(U,b)==QN(0)
      end
      for b in nzblocks(S)
        @test flux(S,b)==QN(1)
      end
      for b in nzblocks(V)
        @test flux(V,b)==QN(0)
      end

      @test minimum(dims(S)) == maxdim
      @test minimum(dims(S)) == length(spec.eigs)
      @test minimum(dims(S)) < dim(i)*dim(j)

      Ap = U*S*V
      err = 1-(Ap*dag(Ap))[]/(A*dag(A))[]
      @test isapprox(err,spec.truncerr; rtol=1e-6)
    end

    @testset "svd truncation example 4" begin
      i = Index(QN(0,2)=>3,QN(1,2)=>4; tags="i")
      j = settags(i,"j")
      A = randomITensor(ElT, QN(1,2),i,j,dag(i'),dag(j'))

      maxdim = 4
      U,S,V,spec = svd(A,i,j; utags="x", vtags="y", maxdim=maxdim)

      @test store(U) isa NDTensors.BlockSparse
      @test store(S) isa NDTensors.DiagBlockSparse
      @test store(V) isa NDTensors.BlockSparse

      u = commonind(S,U)
      v = commonind(S,V)

      @test hastags(u,"x")
      @test hastags(v,"y")

      @test hassameinds(U,(i,j,u))
      @test hassameinds(V,(i',j',v))

      for b in nzblocks(A)
        @test flux(A,b)==QN(1,2)
      end
      for b in nzblocks(U)
        @test flux(U,b)==QN(0,2)
      end
      for b in nzblocks(S)
        @test flux(S,b)==QN(1,2)
      end
      for b in nzblocks(V)
        @test flux(V,b)==QN(0,2)
      end

      @test minimum(dims(S)) == maxdim
      @test minimum(dims(S)) == length(spec.eigs)
      @test minimum(dims(S)) < dim(i)*dim(j)

      Ap = U*S*V
      err = 1-(Ap*dag(Ap))[]/(A*dag(A))[]
      @test isapprox(err,spec.truncerr; rtol=1e-6)
    end

    @testset "svd truncation example 5" begin
      i = Index(QN(0,2)=>2,QN(1,2)=>3; tags="i")
      j = settags(i,"j")
      A = randomITensor(ElT, QN(1,2),i,j,dag(i'),dag(j'))

      maxdim = 4
      U,S,V,spec = svd(A,i,j'; utags="x", vtags="y", maxdim=maxdim)

      @test store(U) isa NDTensors.BlockSparse
      @test store(S) isa NDTensors.DiagBlockSparse
      @test store(V) isa NDTensors.BlockSparse

      u = commonind(S,U)
      v = commonind(S,V)

      @test hastags(u,"x")
      @test hastags(v,"y")

      @test hassameinds(U,(i,j',u))
      @test hassameinds(V,(i',j,v))

      for b in nzblocks(A)
        @test flux(A,b)==QN(1,2)
      end
      for b in nzblocks(U)
        @test flux(U,b)==QN(0,2)
      end
      for b in nzblocks(S)
        @test flux(S,b)==QN(1,2)
      end
      for b in nzblocks(V)
        @test flux(V,b)==QN(0,2)
      end

      @test minimum(dims(S)) == maxdim
      @test minimum(dims(S)) == length(spec.eigs)
      @test minimum(dims(S)) < dim(i)*dim(j)

      Ap = U*S*V
      err = 1-(Ap*dag(Ap))[]/(A*dag(A))[]
      @test isapprox(err,spec.truncerr; rtol=1e-6)
    end

    @testset "issue #231" begin
      l = Index(QN("Nf",-1,-1)=>2,
            QN("Nf", 0,-1)=>4,
            QN("Nf",+1,-1)=>2,
            tags="CMB,Link")
      s = Index(QN("Nf",0,-1)=>1,
                QN("Nf",1,-1)=>1,
                tags="Fermion,Site,n=4")
      r = Index(QN("Nf",1,-1)=>2,
                QN("Nf",0,-1)=>1,
                QN("Nf",1,-1)=>2,
                tags="Link,u")

      A = emptyITensor(ElT, l,s,dag(r))

      addblock!(A,(2,1,2))
      addblock!(A,(1,2,2))
      addblock!(A,(2,2,3))

      for b in nzblocks(A)
        @test flux(A,b)==QN()
      end

      U,S,V = svd(A,l,s)

      for b in nzblocks(U)
        @test flux(U,b)==QN()
      end
      for b in nzblocks(S)
        @test flux(S,b)==QN()
      end
      for b in nzblocks(V)
        @test flux(V,b)==QN()
      end
      @test norm(U*S*V-A) ≈ 0 atol=1e-15
    end

    @testset "SVD no truncate bug" begin
      s = Index(QN("Sz",-4) => 1,
                QN("Sz",-2) => 4,
                QN("Sz", 0) => 6,
                QN("Sz", 2) => 4,
                QN("Sz", 4) => 1)
      A = emptyITensor(ElT, s, s')
      addblock!(A, (5,2))
      addblock!(A, (4,3))
      addblock!(A, (3,4))
      addblock!(A, (2,5))
      randn!(A)
      U,S,V = svd(A,s)
      @test U*S*V ≈ A
    end

    @testset "SVD no truncate" begin
      s = Index(QN("Sz",-4) => 1,
                QN("Sz",-2) => 4,
                QN("Sz", 0) => 6,
                QN("Sz", 2) => 4,
                QN("Sz", 4) => 1)
      A = emptyITensor(ElT, s, s')
      addblock!(A, (5,1))
      addblock!(A, (4,2))
      addblock!(A, (3,3))
      addblock!(A, (2,4))
      addblock!(A, (1,5))
      U,S,V = svd(A, s)
      @test dims(S) == dims(A)
      @test U*S*V ≈ A
    end

    @testset "SVD truncate zeros" begin
      s = Index(QN("Sz",-4) => 1,
                QN("Sz",-2) => 4,
                QN("Sz", 0) => 6,
                QN("Sz", 2) => 4,
                QN("Sz", 4) => 1)
      A = emptyITensor(ElT, s, s')
      addblock!(A, (5,1))
      addblock!(A, (4,2))
      addblock!(A, (3,3))
      addblock!(A, (2,4))
      addblock!(A, (1,5))
      U,S,V = svd(A, s; cutoff=0)
      @test dims(S) == (0,0)
      @test U*S*V ≈ A
    end

  end

  @testset "Replace Index" begin
    i = Index([QN(0)=>1,QN(1)=>2],"i")
    j = Index([QN(0)=>3,QN(1)=>4,QN(2)=>5],"j")

    T1 = randomITensor(QN(1),i,j)
    T2 = copy(T1)

    k = Index([QN(0)=>1,QN(1)=>2],"k")

    replaceind!(T1,i,k)
    @test hasind(T1,k)
    @test dir(inds(T1)[1]) == dir(i)

    # Check that replaceind! keeps
    # original Arrow direction
    replaceind!(T2,i,dag(k))
    @test hasind(T2,k)
    @test dir(inds(T2)[1]) == dir(i)
    @test dir(inds(T2)[1]) != dir(dag(k))
  end

@testset "BlockSparse dag copy behavior" begin
  i = Index(QN(0)=>2,QN(1)=>2,tags="i")
  j = Index(QN(0)=>2,QN(1)=>2,tags="j")

  v1 = randomITensor(QN(1),i,j)
  orig_elt = v1[1,3]
  cv1 = dag(v1;always_copy=false)
  cv1[1,3] = 123.45
  @test v1[1,3] ≈ cv1[1,3]

  v2 = randomITensor(QN(1),i,j)
  orig_elt = v2[1,3]
  cv2 = dag(v2;always_copy=true)
  cv2[1,3] = 123.45
  @test v2[1,3] ≈ orig_elt

  v3 = randomITensor(ComplexF64,QN(1),i,j)
  orig_elt = v3[1,3]
  cv3 = dag(v3;always_copy=false)
  cv3[1,3] = 123.45
  @test v3[1,3] ≈ orig_elt

  v4 = randomITensor(ComplexF64,QN(1),i,j)
  orig_elt = v4[1,3]
  cv4 = dag(v4;always_copy=true)
  cv4[1,3] = 123.45
  @test v4[1,3] ≈ orig_elt

end

@testset "exponentiate" begin

  @testset "Simple arrows" begin
    i1 = Index([QN(0)=>1,QN(1)=>2],"i1")
    i2 = Index([QN(0)=>1,QN(1)=>2],"i2")
    A = randomITensor(QN(),i1,i2,dag(i1)', dag(i2)')
    Aexp = exp(A)
    Amat = Array(A, i1,i2, i1', i2')
    Amatexp = reshape(exp(reshape(Amat,9,9)),3,3,3,3)
    @test isapprox(norm(Array(Aexp,i1,i2,i1',i2') - Amatexp),0.0, atol=1e-14)
    @test flux(Aexp) == QN()
    @test length(setdiff(inds(Aexp),inds(A)))==0

    @test isapprox(norm(exp(A, (i1,i2), (i1',i2')) - Aexp),0.0, atol=5e-14)

    # test the case where indices are permuted
    A = randomITensor(QN(),i1,dag(i1)', dag(i2)',i2)
    Aexp = exp(A, (i1,i2), (i1',i2'))
    Amat = Array(A, i1,i2,i1', i2')
    Amatexp = reshape(exp(reshape(Amat,9,9)),3,3,3,3)
    @test isapprox(norm(Array(Aexp,i1,i2,i1',i2') - Amatexp),0.0,atol=1e-14)

    # test exponentiation in the Hermitian case
    i1 = Index([QN(0)=>2,QN(1)=>2,QN(2)=>3],"i1")
    A = randomITensor(QN(),i1,dag(i1)')
    Ad = dag(swapinds(A,IndexSet(i1),IndexSet(dag(i1)')))
    Ah = A+Ad + 1e-10*randomITensor(QN(),i1,dag(i1)')
    Amat = Array(Ah ,i1', i1)
    Aexp = exp(Ah; ishermitian=true)
    Amatexp= exp(LinearAlgebra.Hermitian(Amat))
    @test isapprox(norm(Array(Aexp,i1,i1')-Amatexp),0.0,atol=5e-14)
  end

  @testset "Mixed arrows" begin
    i1 = Index([QN(0)=>1,QN(1)=>2],"i1")
    i2 = Index([QN(0)=>1,QN(1)=>2],"i2")
    A = randomITensor(i1, i2, dag(i1)', dag(i2)')
    expA = exp(A, (i1, i1'), (i2', i2))
    @test exp(dense(A), (i1, i1'), (i2', i2)) ≈ dense(expA)
  end

end

end

nothing


text:
# ------------ Call sites ------------ #

@inline function (ctx::ParameterBackpropagateContext)(c::typeof(plate),
                                                      addr::T,
                                                      call::Function,
                                                      args::Vector) where T <: Address
    #visit!(ctx.visited, addr)
    vcl = get_sub(ctx.tr, addr)
    param_grads = Gradients()
    params = get_sub(ctx.initial_params, addr)
    ret = simulate_call_pullback(params, param_grads, vcl, args)
    ctx.param_grads.tree[addr] = param_grads
    return ret
end

@inline function (ctx::ChoiceBackpropagateContext)(c::typeof(plate),
                                                   addr::T,
                                                   call::Function,
                                                   args::Vector) where T <: Address
    #visit!(ctx.visited, addr)
    vcl = get_sub(ctx.tr, addr)
    choice_grads = Gradients()
    params = get_sub(ctx.initial_params, addr)
    ret = simulate_choice_pullback(params, choice_grads, get_sub(ctx.select, addr), vcl, args)
    ctx.choice_grads.tree[addr] = choice_grads
    return ret
end

ext: jl


text:
module SPH
    import WAV

    export sphread, sphreadheader, sphwrite, sph2wav

    # isip.piconepress.com/projects/speech/software/tutorials/production/fundamentals/v1.0/section_02/text/nist_sphere.text
    function sphreadheader(io::IO)
        sph = Array{UInt8}(undef, 16)
        read!(io, sph)
        version = sph[1:8]
        header_len = parse(Int, String(sph[9:16]))
        h = Array{UInt8}(undef, header_len - 16)
        read!(io, h)
        header_lines = strip.(split(String(h), "\n"))
        header = Dict{String, Any}()

        for l in header_lines
            l = rstrip(l, '\0')
            if isempty(l)
                continue
            end
            toks = split(l)
            i = 1
            if startswith(toks[i], ';')
                i += 1
            end
            if length(toks[i:end]) < 3
                continue
            end
            (obj, t, val) = toks[i:i+2]
            if t == "-i"
                header[obj] = parse(Int, val)
            elseif t == "-r"
                header[obj] = parse(Float32, val)
            elseif val == "TRUE"
                header[obj] = true
            elseif val == "FALSE"
                header[obj] = false
            else
                header[obj] = val
            end
        end
        return header, header_len
    end

    function sphreadheader(filename::AbstractString)
        open(filename, "r") do io
            sphreadheader(io)
        end
    end

    function get_wav_format(header)
        fmt = get(header, "sample_coding", "pcm")
        cc = get(
            Dict(
                "ulaw" => WAV.WAVE_FORMAT_MULAW,
                "alaw" => WAV.WAVE_FORMAT_ALAW
            ),
            fmt,
            WAV.WAVE_FORMAT_PCM
        )
        fs = header["sample_rate"]
        nchan = header["channel_count"]
        nb = header["sample_n_bytes"] * 8
        ba = header["sample_n_bytes"] * nchan
        bps = fs * ba
        ext = WAV.WAVFormatExtension()
        return WAV.WAVFormat(cc, nchan, fs, bps, ba, nb, ext)
    end

    function sphread(io::IO; subrange=(:), format="double")
        header, header_len = sphreadheader(io)
        seekstart(io)
        data = read(io)[(header_len + 1):end]
        buf = IOBuffer(data)
        fmt = get_wav_format(header)
        # use the internal WAV functions for doing ulaw and alaw
        samples = WAV.read_data(buf, length(data), fmt, format, subrange)
        return header, samples
    end

    function sphread(filename::AbstractString; subrange=(:), format="double")
        open(filename, "r") do io
            sphread(io; subrange=subrange, format=format)
        end
    end

    function sphwrite(header, samples, io::IO)
        header_len = 1024
        bytes_written = 0
        line = "NIST_1A\n"
        write(io, line)
        bytes_written += length(line)
        line = "   $header_len\n"
        write(io, line)
        bytes_written += length(line)
        for obj in keys(header)
            val = header[obj]
            s = ""
            T = typeof(val)
            if T <: Integer
                s = "$val"
                t = "-i"
            elseif T <: Real
                s = "$val"
                t = "-r"
            elseif T <: Bool
                s = uppercase("$val")
                t = "-s$(length(s))"
            else
                @assert T <: AbstractString
                s = val
                t = "-s$(length(s))"
            end
            line = string(join([obj, t, s], " "), "\n")
            write(io, line)
            bytes_written += length(line)
        end
        line = "end_head\n"
        write(io, line)
        bytes_written += length(line)
        write(io, repeat(" ", header_len - bytes_written))
        fmt = get_wav_format(header)
        WAV.write_data(io, fmt, samples)
    end

    function sphwrite(header, samples, filename::AbstractString)
        open(filename, "w") do io
            sphwrite(header, samples, io)
        end
    end

    function sph2wav(sph::IO, wav::IO)
        header, samples = sphread(sph)
        fmt = get_wav_format(header)
        WAV.wavwrite(samples, wav; Fs=fmt.sample_rate, nbits=fmt.nbits, compression=fmt.compression_code)
    end

    function sph2wav(sph_filename::AbstractString, wav_filename::AbstractString)
        open(sph_filename, "r") do sph
            open(wav_filename, "w") do wav
                sph2wav(sph, wav)
            end
        end
    end
end


text:
#Demo 11: Digital buffer "Simulation"
#-------------------------------------------------------------------------------
using InspectDR
using Colors
using Random

#==Input
===============================================================================#

#Constants
#-------------------------------------------------------------------------------
black = RGB24(0, 0, 0)
white = RGB24(1, 1, 1)
red = RGB24(1, 0, 0)
green = RGB24(0, 1, 0)
blue = RGB24(0, 0, 1)

line_default = line(color=blue, width=3)
line_supply = line(color=red, width=5)
line_thresh = line(style=:dash, width=2.5, color=RGB24(.4,.4,.4))

#Input data
#-------------------------------------------------------------------------------
VSUPPLY = 1.0
tbit = 1e-9 #Bit period
osr = 40 #samples per bit

Random.seed!(1234) #Generate consistent "random" patterns
#pat = rand(Bool, 8)*1 #Random bit pattern
pat = [0,0,1,1,0,0,0,1,0,1,0]

nsamples = length(pat)
t = collect(range(0, stop=nsamples*tbit, length=nsamples*osr+1)[1:end-1]) #Time
t_1bit = t[1:osr] #Time for 1 bit

#=="Simulation" (really rudimentary to keep things simple)
===============================================================================#
expdecay(Δt::Vector, V0, V∞, τ) = V∞ .+ (V0-V∞)*exp.(-Δt/τ)

function EmulDigtialWfrm(pat::Vector, τrise, τfall)
	vrise = expdecay(t_1bit, 0, VSUPPLY, τrise)
	vfall = expdecay(t_1bit, VSUPPLY, 0, τfall)
	vsup_1bit = VSUPPLY*fill(1, size(t_1bit))
	nsamples = length(pat)

	#Generate waveform:
	result = Vector{Float64}(undef, nsamples*osr)
	padpat = prepend!(copy(pat), [0])
	delta = padpat[2:end] - padpat[1:end-1] #Difference between consecutive bits (i.e. transitions)
	for ibit in 1:length(pat)
		istart = (ibit-1)*osr
		if +1 == delta[ibit]
			data = vrise
		elseif -1 == delta[ibit]
			data = vfall
		else
			data = pat[ibit]*vsup_1bit
		end

		for ioffset in 1:osr
			result[istart+ioffset] = data[ioffset]
		end
	end

	return result
end

vpath1 = EmulDigtialWfrm(pat, tbit/8, tbit/20)
#Path 2 is properly balanced:
vpath2 = EmulDigtialWfrm(pat, tbit/20, tbit/20)

vgnd = 0*t
vsup = vgnd .+ VSUPPLY

#==Generate plot
===============================================================================#
mplot = InspectDR.Multiplot(title="Digital Buffer \"Simulation\"")
mplot.layout[:ncolumns] = 1
ext_digital = InspectDR.PExtents1D(-0.2, 1.2) #"Digital" signals
v10 = 0.1*VSUPPLY; v90 = 0.9*VSUPPLY #10/90 digital thresholds

#To be plotted:
siginfolist = [
	("path1", vpath1, ext_digital),
	("path2", vpath2, ext_digital),
]

plot = add(mplot, InspectDR.transientplot([:lin for i in siginfolist],
	title="Measure trise/tfall with `r` & `d` keys",
	ylabels=["(V)" for i in siginfolist])
)
plot.layout[:enable_legend] = true
#plot.layout[:halloc_legend] = 150

let wfrm #HIDEWARN_0.7
	for (i, siginfo) in enumerate(siginfolist)
		id, sig, ext = siginfo
		plot.strips[i].yext_full = ext
		plot.strips[i].grid = InspectDR.GridRect() #No grid
#		plot.strips[i].grid = InspectDR.GridRect(vmajor=true, vminor=true, hmajor=false, hminor=false)

		#Draw supply lines under data:
		wfrm = add(plot, t, vgnd, strip=i)
			wfrm.line = line_supply
		wfrm = add(plot, t, vsup, strip=i)
			wfrm.line = line_supply

		#Draw current signal:
		wfrm = add(plot, t, sig, id=id, strip=i)
			wfrm.line = line_default

		#Add threshold markers:
		add(plot, hmarker(v10, line_thresh, strip=i))
		add(plot, hmarker(v90, line_thresh, strip=i))
	end
end

gplot = display(InspectDR.GtkDisplay(), mplot)

mplot.layout[:halloc_plot] = 700; mplot.layout[:valloc_plot] = 450
InspectDR.write_png("export_simdig.png", mplot)

:DONE


text:
function construct_network!(
    psi_container::PSIContainer,
    sys::PSY.System,
    ::Type{CopperPlatePowerModel},
)
    buses = PSY.get_components(PSY.Bus, sys)
    bus_count = length(buses)

    get_slack_variables(psi_container.settings) &&
        add_slacks!(psi_container, CopperPlatePowerModel)
    copper_plate(psi_container, :nodal_balance_active, bus_count)

    return
end

function construct_network!(
    psi_container::PSIContainer,
    sys::PSY.System,
    ::Type{StandardPTDFModel},
)
    buses = PSY.get_components(PSY.Bus, sys)
    ac_branches = get_available_components(PSY.ACBranch, sys)
    ptdf = get_PTDF(psi_container)

    if isnothing(ptdf)
        throw(ArgumentError("no PTDF matrix supplied"))
    end

    get_slack_variables(psi_container.settings) &&
        add_slacks!(psi_container, StandardPTDFModel)

    ptdf_networkflow(psi_container, ac_branches, buses, :nodal_balance_active, ptdf)

    dc_branches = get_available_components(PSY.DCBranch, sys)
    dc_branch_types = typeof.(dc_branches)
    for btype in Set(dc_branch_types)
        typed_dc_branches = IS.FlattenIteratorWrapper(
            btype,
            Vector([[b for b in dc_branches if typeof(b) == btype]]),
        )
        flow_variables!(psi_container, StandardPTDFModel, typed_dc_branches)
    end
    return
end

function construct_network!(
    psi_container::PSIContainer,
    sys::PSY.System,
    ::Type{T},
) where {T <: PM.AbstractPowerModel}
    incompat_list = [
        PM.SDPWRMPowerModel,
        PM.SparseSDPWRMPowerModel,
        PM.SOCBFPowerModel,
        PM.SOCBFConicPowerModel,
    ]
    if T in incompat_list
        throw(ArgumentError("$(T) formulation is not currently supported in PowerSimulations"))
    end

    get_slack_variables(psi_container.settings) && add_slacks!(psi_container, T)

    powermodels_network!(psi_container, T, sys)
    add_pm_var_refs!(psi_container, T, sys)
    return
end


text:
using Mux

@app test = (                                       # 1
  Mux.defaults,
  page(respond("
Hello World from Mux!
")),
  page("/about", respond("
About Mux
")),    # 2
  page("/user/:user", req -> "
Hello, $(req[:params][:user])!
"),
  Mux.notfound())

serve(test)                                         # 3


text:
# Add support for analytical derivatives of min and max using Heaviside function (taken from SymPy)
function differentiate(::SymbolParameter{:min}, args, wrt)
   k = 1
   derivative = "("*string(differentiate(args[k],wrt))*")"
   derivative *= "*heaviside(min("*replace(replace(string(args[[1:(k-1), (k+1):end]]), r"\{|\}", ""), ":","")*") - ($(args[k]))) + "
   for k in 2:length(args)
       derivative *= string(differentiate(args[k],wrt))
       derivative *= "*heaviside(min("*replace(replace(string(args[[1:(k-1), (k+1):end]]), r"\{|\}", ""), ":","")*") - ($(args[k]))) + "
   end
   parse(derivative[1:(end - 2)])
end

function differentiate(::SymbolParameter{:max}, args, wrt)
   k = 1
   derivative = "("*string(differentiate(args[k],wrt))*")"
   derivative *= "*heaviside(-max("*replace(replace(string(args[[1:(k-1), (k+1):end]]), r"\{|\}", ""), ":","")*") + ($(args[k]))) + "
   for k in 2:length(args)
       derivative *= string(differentiate(args[k],wrt))
       derivative *= "*heaviside(-max("*replace(replace(string(args[[1:(k-1), (k+1):end]]), r"\{|\}", ""), ":","")*") + ($(args[k]))) + "
   end
   parse(derivative[1:(end - 2)])
end

function differentiate(::SymbolParameter{:heaviside}, args, wrt)
   derivative = "("*string(differentiate(args[1],wrt))*")"
   derivative *= "*dirac("*replace(replace(string(args), r"\{|\}", ""), ":","")*")"
   parse(derivative)
end

differentiate(::SymbolParameter{:dirac}, args, wrt) = 0


text:
import Base.getindex

type DelayedBuffer
    data::Array{Float64,1}
    w_len::Int
    w_inc::Int
    delay::Int
    z::Array{Float64,1}

    function DelayedBuffer(data::Array{Float64,1},w_len::Int,w_inc::Int)
        delay = w_len-w_inc
        z = zeros(w_len)
        new(data,w_len,w_inc,delay,z)
    end
end

function getindex(x::DelayedBuffer, i::Integer)
    ll = 1+(i-1)*x.w_inc-x.delay
    rr = (i-1)*x.w_inc+x.w_len-x.delay
    if ll < 1 # accessing delay indexes
        x.z[1:abs(ll)+1] = 0
        x.z[abs(ll)+2:end] = sub(x.data,1:rr)
    elseif ll > length(x.data)-x.delay
        error("BoundsError")
    elseif rr > length(x.data)
        x.z[1:length(x.data)-ll+1] = sub(x.data,ll:end)
        x.z[length(x.data)-ll+2:end] = 0
    else
        return sub(x.data,ll:rr)
    end
    return x.z
end


text:
module FluxRNNsBench

using ReTest
using FluxRNNs
using Random
using CUDA
using Statistics
using ProgressMeter
using DataFrames
using CSV
using Git

include("bench_utils.jl")

function git_head()
    try
        s = read(`$(Git.git()) rev-parse HEAD`, String)
        s[1:end-1]
    catch
        "0"
    end
end

# rng = Random.MersenneTwister(10)

function get_statistics(t)
    tms = t.times
    mean(tms), var(tms), median(tms), minimum(tms), maximum(tms), t.allocs, t.memory
end

function get_runtimes(create_features::Function, cell; Ns = [2, 20, 200], Ts = [1, 8, 16, 64], cuda=false)

    RUNTIMES = Dict(k=>DataFrames.DataFrame(
        cell=String[],
        gpu=Bool[],
        n=Int[],
        t=Int[],
        mean=Float64[],
        variance=Float64[],
        median=Float64[],
        minimum=Float64[],
        maximum=Float64[],
        allocs=Int[],
        memory=Int[]) for k in [:fw, :bw, :fwbw])

    Random.seed!(10)

    @showprogress "MatrixSize: " for n in Ns
        @showprogress "TemporalLength: " 1 for t in Ts
            x = create_features(n, t)
            rnn_model = cell(n, n)

            timings = run_benchmark_recur(rnn_model, x, cuda=cuda)
            for k in keys(timings)
                push!(RUNTIMES[k], ("$(cell)", cuda, n, t, get_statistics(timings[k])...))
            end
        end
    end
    RUNTIMES
end

function run_benchmark(cells=[RNN, GRU, LSTM]; name="$(git_head())", cuda=false, kwargs...)

    save_dir = joinpath("results", name)
    if !isdir(save_dir)
        mkpath(save_dir)
    end

    for cell in cells
        println("""$(cell)-$(cuda ? "GPU" : "CPU")""")
        cell_svpath = joinpath(save_dir, """$(cell)-$(cuda ? "GPU" : "CPU")""")
        if !isdir(cell_svpath)
            mkdir(cell_svpath)
        end
        rt = get_runtimes(cell; cuda=cuda, kwargs...) do n, t
            [rand(Float32, n, n) for _ in 1:t]
        end
        for k in keys(rt)
            filename = joinpath(cell_svpath, "$(k).csv")
            CSV.write(filename, rt[k])
        end
        println()
        println("$(cell)3d-CPU")
        rt = get_runtimes(cell; cuda=cuda, kwargs...) do n, t
            rand(Float32, n, n, t)
        end
        for k in keys(rt)
            filename = joinpath(cell_svpath, "$(k)_3d.csv")
            CSV.write(filename, rt[k])
        end
    end
end

end


text:
module Solvers

export
    BatchedTridiagonalSolver, solve!,
    FFTBasedPoissonSolver, FourierTridiagonalPoissonSolver, PreconditionedConjugateGradientSolver, PressureSolver,
    solve_for_pressure!, solve_poisson_equation!

using Statistics
using FFTW
using CUDA
using KernelAbstractions
using KernelAbstractions.Extras.LoopInfo: @unroll

using Oceananigans.Architectures: device, CPU, GPU, array_type, arch_array
using Oceananigans.Utils
using Oceananigans.Grids
using Oceananigans.BoundaryConditions
using Oceananigans.Fields

using Oceananigans.Grids: unpack_grid

"""
    ω(M, k)

Return the `M`th root of unity raised to the `k`th power.
"""
@inline ω(M, k) = exp(-2im*π*k/M)

reshaped_size(N, dim) = dim == 1 ? (N, 1, 1) :
                        dim == 2 ? (1, N, 1) :
                        dim == 3 ? (1, 1, N) : nothing

include("batched_tridiagonal_solver.jl")
include("poisson_eigenvalues.jl")
include("index_permutations.jl")
include("discrete_transforms.jl")
include("plan_transforms.jl")
include("fft_based_poisson_solver.jl")
include("fourier_tridiagonal_poisson_solver.jl")
include("preconditioned_conjugate_gradient_solver.jl")
include("solve_for_pressure.jl")

PressureSolver(arch, grid::RegularRectilinearGrid) = FFTBasedPoissonSolver(arch, grid)
PressureSolver(arch, grid::VerticallyStretchedRectilinearGrid) = FourierTridiagonalPoissonSolver(arch, grid)

end


text:
using JPPF
using Test

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


text:
function convert(::Type{Project{Existing}}, gen::ThermalGenEMIS{BuildPhase})
    return ThermalGenEMIS{Existing}(gen.name,
                               gen.tech,
                               gen.decision_year,
                               gen.construction_year,
                               gen.retirement_year,
                               gen.end_life_year,
                               gen.products,
                               gen.finance_data)
end

function convert(::Type{Project{Planned}}, gen::ThermalGenEMIS{BuildPhase})
    return ThermalGenEMIS{Planned}(gen.name,
                               gen.tech,
                               gen.decision_year,
                               gen.construction_year,
                               gen.retirement_year,
                               gen.end_life_year,
                               gen.products,
                               gen.finance_data)
end

function convert(::Type{Project{Queue}}, gen::ThermalGenEMIS{BuildPhase})
    return ThermalGenEMIS{Queue}(gen.name,
                               gen.tech,
                               gen.decision_year,
                               gen.construction_year,
                               gen.retirement_year,
                               gen.end_life_year,
                               gen.products,
                               gen.finance_data)
end

function convert(::Type{Project{Option}}, gen::ThermalGenEMIS{BuildPhase})
    return ThermalGenEMIS{Option}(gen.name,
                               gen.tech,
                               gen.decision_year,
                               gen.construction_year,
                               gen.retirement_year,
                               gen.end_life_year,
                               gen.products,
                               gen.finance_data)
end

function convert(::Type{Project{Retired}}, gen::ThermalGenEMIS{BuildPhase})
    return ThermalGenEMIS{Retired}(gen.name,
                               gen.tech,
                               gen.decision_year,
                               gen.construction_year,
                               gen.retirement_year,
                               gen.end_life_year,
                               gen.products,
                               gen.finance_data)
end

function convert(::Type{Project{Existing}}, gen::RenewableGenEMIS{BuildPhase})
    return RenewableGenEMIS{Existing}(gen.name,
                               gen.tech,
                               gen.decision_year,
                               gen.construction_year,
                               gen.retirement_year,
                               gen.end_life_year,
                               gen.products,
                               gen.finance_data)
end

function convert(::Type{Project{Planned}}, gen::RenewableGenEMIS{BuildPhase})
    return RenewableGenEMIS{Planned}(gen.name,
                               gen.tech,
                               gen.decision_year,
                               gen.construction_year,
                               gen.retirement_year,
                               gen.end_life_year,
                               gen.products,
                               gen.finance_data)
end

function convert(::Type{Project{Queue}}, gen::RenewableGenEMIS{BuildPhase})
    return RenewableGenEMIS{Queue}(gen.name,
                               gen.tech,
                               gen.decision_year,
                               gen.construction_year,
                               gen.retirement_year,
                               gen.end_life_year,
                               gen.products,
                               gen.finance_data)
end

function convert(::Type{Project{Option}}, gen::RenewableGenEMIS{BuildPhase})
    return RenewableGenEMIS{Option}(gen.name,
                               gen.tech,
                               gen.decision_year,
                               gen.construction_year,
                               gen.retirement_year,
                               gen.end_life_year,
                               gen.products,
                               gen.finance_data)
end

function convert(::Type{Project{Retired}}, gen::RenewableGenEMIS{BuildPhase})
    return RenewableGenEMIS{Retired}(gen.name,
                               gen.tech,
                               gen.decision_year,
                               gen.construction_year,
                               gen.retirement_year,
                               gen.end_life_year,
                               gen.products,
                               gen.finance_data)
end

function convert(::Type{Project{Existing}}, gen::HydroGenEMIS{BuildPhase})
    return HydroGenEMIS{Existing}(gen.name,
                               gen.tech,
                               gen.decision_year,
                               gen.construction_year,
                               gen.retirement_year,
                               gen.end_life_year,
                               gen.products,
                               gen.finance_data)
end

function convert(::Type{Project{Planned}}, gen::HydroGenEMIS{BuildPhase})
    return HydroGenEMIS{Planned}(gen.name,
                               gen.tech,
                               gen.decision_year,
                               gen.construction_year,
                               gen.retirement_year,
                               gen.end_life_year,
                               gen.products,
                               gen.finance_data)
end

function convert(::Type{Project{Queue}}, gen::HydroGenEMIS{BuildPhase})
    return HydroGenEMIS{Queue}(gen.name,
                               gen.tech,
                               gen.decision_year,
                               gen.construction_year,
                               gen.retirement_year,
                               gen.end_life_year,
                               gen.products,
                               gen.finance_data)
end

function convert(::Type{Project{Option}}, gen::HydroGenEMIS{BuildPhase})
    return HydroGenEMIS{Option}(gen.name,
                               gen.tech,
                               gen.decision_year,
                               gen.construction_year,
                               gen.retirement_year,
                               gen.end_life_year,
                               gen.products,
                               gen.finance_data)
end

function convert(::Type{Project{Retired}}, gen::HydroGenEMIS{BuildPhase})
    return HydroGenEMIS{Retired}(gen.name,
                               gen.tech,
                               gen.decision_year,
                               gen.construction_year,
                               gen.retirement_year,
                               gen.end_life_year,
                               gen.products,
                               gen.finance_data)
end

function convert(::Type{Project{Existing}}, battery::BatteryEMIS{BuildPhase})
    return BatteryEMIS{Existing}(battery.name,
                               battery.tech,
                               battery.decision_year,
                               battery.construction_year,
                               battery.retirement_year,
                               battery.end_life_year,
                               battery.products,
                               battery.finance_data)
end

function convert(::Type{Project{Planned}}, battery::BatteryEMIS{BuildPhase})
    return BatteryEMIS{Planned}(battery.name,
                               battery.tech,
                               battery.decision_year,
                               battery.construction_year,
                               battery.retirement_year,
                               battery.end_life_year,
                               battery.products,
                               battery.finance_data)
end

function convert(::Type{Project{Queue}}, battery::BatteryEMIS{BuildPhase})
    return BatteryEMIS{Queue}(battery.name,
                               battery.tech,
                               battery.decision_year,
                               battery.construction_year,
                               battery.retirement_year,
                               battery.end_life_year,
                               battery.products,
                               battery.finance_data)
end

function convert(::Type{Project{Option}}, battery::BatteryEMIS{BuildPhase})
    return BatteryEMIS{Option}(battery.name,
                               battery.tech,
                               battery.decision_year,
                               battery.construction_year,
                               battery.retirement_year,
                               battery.end_life_year,
                               battery.products,
                               battery.finance_data)
end

function convert(::Type{Project{Retired}}, battery::BatteryEMIS{BuildPhase})
    return BatteryEMIS{Retired}(battery.name,
                               battery.tech,
                               battery.decision_year,
                               battery.construction_year,
                               battery.retirement_year,
                               battery.end_life_year,
                               battery.products,
                               battery.finance_data)
end

function convert(::Type{Project{Queue}}, gen::ThermalGenEMIS{Option})
    return ThermalGenEMIS{Queue}(gen.name,
                            gen.tech,
                            gen.decision_year,
                            gen.construction_year,
                            gen.retirement_year,
                            gen.end_life_year,
                            gen.products,
                            gen.finance_data)
end

function convert(::Type{Project{Planned}}, gen::ThermalGenEMIS{Queue})
    return ThermalGenEMIS{Planned}(gen.name,
                              gen.tech,
                              gen.decision_year,
                              gen.construction_year,
                              gen.retirement_year,
                              gen.end_life_year,
                              gen.products,
                              gen.finance_data)
end

function convert(::Type{Project{Existing}}, gen::ThermalGenEMIS{Planned})
    return ThermalGenEMIS{Existing}(gen.name,
                              gen.tech,
                              gen.decision_year,
                              gen.construction_year,
                              gen.retirement_year,
                              gen.end_life_year,
                              gen.products,
                              gen.finance_data)
end

function convert(::Type{Project{Retired}}, gen::ThermalGenEMIS{Existing})
    return ThermalGenEMIS{Retired}(gen.name,
                              gen.tech,
                              gen.decision_year,
                              gen.construction_year,
                              gen.retirement_year,
                              gen.end_life_year,
                              gen.products,
                              gen.finance_data)
end

function convert(::Type{Project{Retired}}, gen::ThermalGenEMIS{Planned})
    return ThermalGenEMIS{Retired}(gen.name,
                              gen.tech,
                              gen.decision_year,
                              gen.construction_year,
                              gen.retirement_year,
                              gen.end_life_year,
                              gen.products,
                              gen.finance_data)
end

function convert(::Type{Project{Retired}}, gen::ThermalGenEMIS{Queue})
    return ThermalGenEMIS{Retired}(gen.name,
                              gen.tech,
                              gen.decision_year,
                              gen.construction_year,
                              gen.retirement_year,
                              gen.end_life_year,
                              gen.products,
                              gen.finance_data)
end

function convert(::Type{Project{Queue}}, gen::RenewableGenEMIS{Option})
    return RenewableGenEMIS{Queue}(gen.name,
                            gen.tech,
                            gen.decision_year,
                            gen.construction_year,
                            gen.retirement_year,
                            gen.end_life_year,
                            gen.products,
                            gen.finance_data)
end

function convert(::Type{Project{Planned}}, gen::RenewableGenEMIS{Queue})
    return RenewableGenEMIS{Planned}(gen.name,
                              gen.tech,
                              gen.decision_year,
                              gen.construction_year,
                              gen.retirement_year,
                              gen.end_life_year,
                              gen.products,
                              gen.finance_data)
end

function convert(::Type{Project{Existing}}, gen::RenewableGenEMIS{Planned})
    return RenewableGenEMIS{Existing}(gen.name,
                              gen.tech,
                              gen.decision_year,
                              gen.construction_year,
                              gen.retirement_year,
                              gen.end_life_year,
                              gen.products,
                              gen.finance_data)
end

function convert(::Type{Project{Retired}}, gen::RenewableGenEMIS{Existing})
    return RenewableGenEMIS{Retired}(gen.name,
                              gen.tech,
                              gen.decision_year,
                              gen.construction_year,
                              gen.retirement_year,
                              gen.end_life_year,
                              gen.products,
                              gen.finance_data)
end

function convert(::Type{Project{Retired}}, gen::RenewableGenEMIS{Planned})
    return RenewableGenEMIS{Retired}(gen.name,
                              gen.tech,
                              gen.decision_year,
                              gen.construction_year,
                              gen.retirement_year,
                              gen.end_life_year,
                              gen.products,
                              gen.finance_data)
end

function convert(::Type{Project{Retired}}, gen::RenewableGenEMIS{Queue})
    return RenewableGenEMIS{Retired}(gen.name,
                              gen.tech,
                              gen.decision_year,
                              gen.construction_year,
                              gen.retirement_year,
                              gen.end_life_year,
                              gen.products,
                              gen.finance_data)
end

function convert(::Type{Project{Queue}}, gen::HydroGenEMIS{Option})
    return HydroGenEMIS{Queue}(gen.name,
                            gen.tech,
                            gen.decision_year,
                            gen.construction_year,
                            gen.retirement_year,
                            gen.end_life_year,
                            gen.products,
                            gen.finance_data)
end

function convert(::Type{Project{Planned}}, gen::HydroGenEMIS{Queue})
    return HydroGenEMIS{Planned}(gen.name,
                              gen.tech,
                              gen.decision_year,
                              gen.construction_year,
                              gen.retirement_year,
                              gen.end_life_year,
                              gen.products,
                              gen.finance_data)
end

function convert(::Type{Project{Existing}}, gen::HydroGenEMIS{Planned})
    return HydroGenEMIS{Existing}(gen.name,
                              gen.tech,
                              gen.decision_year,
                              gen.construction_year,
                              gen.retirement_year,
                              gen.end_life_year,
                              gen.products,
                              gen.finance_data)
end

function convert(::Type{Project{Retired}}, gen::HydroGenEMIS{Existing})
    return HydroGenEMIS{Retired}(gen.name,
                              gen.tech,
                              gen.decision_year,
                              gen.construction_year,
                              gen.retirement_year,
                              gen.end_life_year,
                              gen.products,
                              gen.finance_data)
end

function convert(::Type{Project{Retired}}, gen::HydroGenEMIS{Planned})
    return HydroGenEMIS{Retired}(gen.name,
                              gen.tech,
                              gen.decision_year,
                              gen.construction_year,
                              gen.retirement_year,
                              gen.end_life_year,
                              gen.products,
                              gen.finance_data)
end

function convert(::Type{Project{Retired}}, gen::HydroGenEMIS{Queue})
    return HydroGenEMIS{Retired}(gen.name,
                              gen.tech,
                              gen.decision_year,
                              gen.construction_year,
                              gen.retirement_year,
                              gen.end_life_year,
                              gen.products,
                              gen.finance_data)
end

function convert(::Type{Project{Queue}}, battery::BatteryEMIS{Option})
    return BatteryEMIS{Queue}(battery.name,
                            battery.tech,
                            battery.decision_year,
                            battery.construction_year,
                            battery.retirement_year,
                            battery.end_life_year,
                            battery.products,
                            battery.finance_data)
end

function convert(::Type{Project{Planned}}, battery::BatteryEMIS{Queue})
    return BatteryEMIS{Planned}(battery.name,
                            battery.tech,
                            battery.decision_year,
                            battery.construction_year,
                            battery.retirement_year,
                            battery.end_life_year,
                            battery.products,
                            battery.finance_data)
end

function convert(::Type{Project{Existing}}, battery::BatteryEMIS{Planned})
    return BatteryEMIS{Existing}(battery.name,
                            battery.tech,
                            battery.decision_year,
                            battery.construction_year,
                            battery.retirement_year,
                            battery.end_life_year,
                            battery.products,
                            battery.finance_data)
end

function convert(::Type{Project{Retired}}, battery::BatteryEMIS{Existing})
    return BatteryEMIS{Retired}(battery.name,
                            battery.tech,
                            battery.decision_year,
                            battery.construction_year,
                            battery.retirement_year,
                            battery.end_life_year,
                            battery.products,
                            battery.finance_data)
end

function convert(::Type{Project{Retired}}, battery::BatteryEMIS{Planned})
    return BatteryEMIS{Retired}(battery.name,
                            battery.tech,
                            battery.decision_year,
                            battery.construction_year,
                            battery.retirement_year,
                            battery.end_life_year,
                            battery.products,
                            battery.finance_data)
end

function convert(::Type{Project{Retired}}, battery::BatteryEMIS{Queue})
    return BatteryEMIS{Retired}(battery.name,
                            battery.tech,
                            battery.decision_year,
                            battery.construction_year,
                            battery.retirement_year,
                            battery.end_life_year,
                            battery.products,
                            battery.finance_data)
end


text:
# ---
# title: 238. Product of Array Except Self
# id: problem238
# author: Indigo
# date: 2021-06-03
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/product-of-array-except-self/description/>
# hidden: true
# ---
#
# Given an array `nums` of _n_ integers where _n_ > 1,  return an array `output`
# such that `output[i]` is equal to the product of all the elements of `nums`
# except `nums[i]`.
#
# **Example:**
#
#
#
#     Input:  [1,2,3,4]
#     Output: [24,12,8,6]
#
#
# **Constraint:**  It's guaranteed that the product of the elements of any
# prefix or suffix of the array (including the whole array) fits in a 32 bit
# integer.
#
# **Note:** Please solve it **without division** and in O( _n_ ).
#
# **Follow up:**
# Could you solve it with constant space complexity? (The output array **does
# not** count as extra space for the purpose of space complexity analysis.)
#
#
## @lc code=start
using LeetCode

function product_except_self(nums::Vector{Int})
    n = length(nums)
    left = right = 1
    res = fill(1, n)
    @inbounds for i in 1:n
        res[i] *= left
        left *= nums[i]
        res[n - i + 1] *= right
        right *= nums[n - i + 1]
    end
    return res
end
## @lc code=end


text:
# Information structure for linear modeling opeartors
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January, 2017
#

const IntNum = Union{Integer, Tuple{Integer,Integer}, Tuple{Integer,Integer,Integer}}

export Info, compareInfo

# Object for velocity/slowness models
mutable struct Info
    n::IntNum
    nsrc::Integer
    nt::Array{Any,1}
end

# Constructor if nt is not passed as a cell
"""
    Info
        n::Integer
        nsrc::Integer
        nt::Array{Any,1}

Info structure that contains information from which the dimensions of Modeling\\
and Projection operators can be inferred.

Constructor
===========

Input arguments are total number of grid points `n`, number of source positions `nsrc` and\\
number of computational time steps `nt` (either as single integer or cell array):

    Info(n, nsrc, nt)

"""
function Info(n::IntNum, nsrc::Integer, nt::Integer)
    ntCell = Array{Any}(undef, nsrc)
    for j=1:nsrc
        ntCell[j] = nt
    end
    return Info(n,nsrc,ntCell)
end

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

function compareInfo(info_A, info_B)
# Compare two info structures. Return true if structures are the same and false otherwise
    if isequal(info_A.n, info_B.n) && isequal(info_A.nsrc, info_B.nsrc) && isequal(info_A.nt, info_B.nt)
        return true
    else
        return false
    end
end


text:
"""
    function supermilker()

A cow being milked, probably from Lars Smith (lars@csua.berkeley.edu)

# Example
```jldoctest
julia> cowsay("Paying the bills", cow=Cowsay.supermilker)
 __________________
< Paying the bills >
 ------------------
  \\   ^__^
   \\  (oo)\\_______        ________
      (__)\\       )\\/\\    |Super |
          ||----W |       |Milker|
          ||    UDDDDDDDDD|______|
```
"""
function supermilker(;eyes="oo", tongue="  ", thoughts="\\")
the_cow =
"""
  $thoughts   ^__^
   $thoughts  ($eyes)\\_______        ________
      (__)\\       )\\/\\    |Super |
       $tongue ||----W |       |Milker|
          ||    UDDDDDDDDD|______|
"""
return the_cow
end


text:
# Sioux Falls network data
# http://www.bgu.ac.il/~bargera/tntp/

#Link travel time = free flow time * ( 1 + B * (flow/capacity)^Power ).
#Link generalized cost = Link travel time + toll_factor * toll + distance_factor * distance

# Traffic Assignment Data structure
type TA_Data
    network_name::String

    number_of_zones::Int64
    number_of_nodes::Int64
    first_thru_node::Int64
    number_of_links::Int64

    start_node::Array
    end_node::Array
    capacity::Array
    link_length::Array
    free_flow_time::Array
    B::Array
    power::Array
    speed_limit::Array
    toll::Array
    link_type::Array

    total_od_flow::Float64

    travel_demand::Array
    od_pairs::Array

    toll_factor::Float64
    distance_factor::Float64

    best_objective::Float64
end

function load_ta_network(network_name="Sioux Falls")

    toll_factor = 0
    distance_factor = 0

    if network_name == "Sioux Falls"
        network_data_file = "SiouxFalls_net.txt"
        trip_table_file = "SiouxFalls_trips.txt"
        best_objective = 4.231335287107440e6 #42.31335287107440
    elseif network_name == "Sioux Falls simplified"
        network_data_file = "SiouxFalls_net_simplified.txt"
        trip_table_file = "SiouxFalls_trips_simplified.txt"
        best_objective = 4.231335287107440e6 #42.31335287107440
    elseif network_name == "Barcelona"
        network_data_file = "Barcelona_net.txt"
        trip_table_file = "Barcelona_trips.txt"
        best_objective = 1
    elseif network_name =="Chicago Sketch"
        network_data_file = "ChicagoSketch_net.txt"
        trip_table_file = "ChicagoSketch_trips.txt"
        best_objective = 1
        toll_factor = 0.02
        distance_factor = 0.04
    elseif network_name == "Anaheim"
        network_data_file = "Anaheim_net.txt"
        trip_table_file = "Anaheim_trips.txt"
        best_objective = 1
    elseif network_name == "Winnipeg"
        network_data_file = "Winnipeg_net.txt"
        trip_table_file = "Winnipeg_trips.txt"
        best_objective = 1
    end

    network_data_file = joinpath("/home/jzh/Dropbox/Research/Data-driven_estimation_inverse_optimization/TrafficAssignment.jl", "data", network_data_file)
    trip_table_file = joinpath("/home/jzh/Dropbox/Research/Data-driven_estimation_inverse_optimization/TrafficAssignment.jl", "data", trip_table_file)

    ##################################################
    # Network Data
    ##################################################

    number_of_zones = 0
    number_of_links = 0
    number_of_nodes = 0
    first_thru_node = 0

    n = open(network_data_file, "r")

    while (line=readline(n)) != ""
        if contains(line, "<NUMBER OF ZONES>")
            number_of_zones = parseint( line[ search(line, '>')+1 : end-1 ] )
        elseif contains(line, "<NUMBER OF NODES>")
            number_of_nodes = parseint( line[ search(line, '>')+1 : end-1 ] )
        elseif contains(line, "<FIRST THRU NODE>")
            first_thru_node = parseint( line[ search(line, '>')+1 : end-1 ] )
        elseif contains(line, "<NUMBER OF LINKS>")
            number_of_links = parseint( line[ search(line, '>')+1 : end-1 ] )
        elseif contains(line, "<END OF METADATA>")
            break
        end
    end

    @assert number_of_links > 0

    start_node = Array(Int, number_of_links)
    end_node = Array(Int, number_of_links)
    capacity = zeros(number_of_links)
    link_length = zeros(number_of_links)
    free_flow_time = zeros(number_of_links)
    B = zeros(number_of_links)
    power = zeros(number_of_links)
    speed_limit = zeros(number_of_links)
    toll = zeros(number_of_links)
    link_type = Array(Int, number_of_links)

    idx = 1
    while (line=readline(n)) != ""
        if contains(line, "~")
            continue
        end

        if contains(line, ";")
            line = strip(line, '\n')
            line = strip(line, ';')

            numbers = split(line)

            start_node[idx] = parseint(numbers[1])
            end_node[idx] = parseint(numbers[2])
            capacity[idx] = parsefloat(numbers[3])
            link_length[idx] = parsefloat(numbers[4])
            free_flow_time[idx] = parsefloat(numbers[5])
            B[idx] = parsefloat(numbers[6])
            power[idx] = parsefloat(numbers[7])
            speed_limit[idx] = parsefloat(numbers[8])
            toll[idx] = parsefloat(numbers[9])
            link_type[idx] = parseint(numbers[10])

            idx = idx + 1
        end
    end

    ##################################################
    # Trip Table
    ##################################################

    number_of_zones_trip = 0
    total_od_flow = 0

    f = open(trip_table_file, "r")

    while (line=readline(f)) != ""
        if contains(line, "<NUMBER OF ZONES>")
            number_of_zones_trip = parseint( line[ search(line, '>')+1 : end-1 ] )
        elseif contains(line, "<TOTAL OD FLOW>")
            total_od_flow = parsefloat( line[ search(line, '>')+1 : end-1 ] )
        elseif contains(line, "<END OF METADATA>")
            break
        end
    end

    @assert number_of_zones_trip == number_of_zones # Check if number_of_zone is same in both txt files
    @assert total_od_flow > 0

    travel_demand = zeros(number_of_zones, number_of_zones)
    od_pairs = Array((Int64,Int64),0)
    while (line=readline(f)) != ""
        if contains(line, "Origin")
            origin = parseint( split(line)[2] )
        elseif contains(line, ";")
            pairs = split(line, ";")
            for i=1:size(pairs)[1]
                if contains(pairs[i], ":")
                    pair = split(pairs[i], ":")
                    destination = parseint( strip(pair[1]) )
                    od_flow = parsefloat( strip(pair[2]) )
                    travel_demand[origin, destination] = od_flow
                    push!(od_pairs, (origin, destination))
                    # println("origin=$origin, destination=$destination, flow=$od_flow")
                end
            end
        end
    end

    # Preparing data to return
    ta_data = TA_Data(
        network_name,
        number_of_zones,
        number_of_nodes,
        first_thru_node,
        number_of_links,
        start_node,
        end_node,
        capacity,
        link_length,
        free_flow_time,
        B,
        power,
        speed_limit,
        toll,
        link_type,
        total_od_flow,
        travel_demand,
        od_pairs,
        toll_factor,
        distance_factor,
        best_objective)

    return ta_data

end # end of load_network function

max_stars_repo_name: jingzbu/InverseVITraffic

max_issues_repo_name: jingzbu/InverseVITraffic

max_forks_repo_name: jingzbu/InverseVITraffic


text:
function gather_kpts_scfres(scfres::NamedTuple)
    # TODO Rename this to gather_kpts once scfres gets its proper type

    # Need gathering over k-Points:
    kpt_properties = (:ψ, :occupation, :eigenvalues)
    scfdict = Dict{Symbol, Any}()
    for symbol in kpt_properties
        scfdict[symbol] = gather_kpts(getproperty(scfres, symbol), scfres.basis)
    end
    scfdict[:basis] = gather_kpts(scfres.basis)

    if mpi_master()
        (; (symbol => get(scfdict, symbol, getproperty(scfres, symbol))
            for symbol in propertynames(scfres))...)
    else
        nothing
    end
end

"""
    load_scfres(filename)

Load back an `scfres`, which has previously been stored with `save_scfres`.
Note the warning in `save_scfres`.
"""
function load_scfres end

"""
    save_scfres(filename, scfres)

Save an `scfres` obtained from `self_consistent_field` to a JLD2
or VTK file depending on the extension.

!!! warning "No compatibility guarantees"
    No guarantees are made with respect to this function at this point.
    It may change incompatibly between DFTK versions or stop working / be removed
    in the future.
"""
function save_scfres(filename::AbstractString, scfres::NamedTuple; kwargs...)
    _, ext = splitext(filename)
    ext = Symbol(ext[2:end])

    # Whitelist valid extensions
    !(ext in (:jld2, :vts)) && error("Extension '$ext' not supported by DFTK.")

    # Gather scfres data on master MPI process
    scfres = gather_kpts_scfres(scfres)

    # On master MPI process dispatch to individual functions based on extension
    ret = isnothing(scfres) ? nothing : save_scfres_master(filename, scfres, Val(ext); kwargs...)

    # Return from this function synchronously on all processes
    # (to ensure that data has actually been written upon return)
    MPI.Barrier(MPI.COMM_WORLD)
    ret
end

function save_scfres_master(filename::AbstractString, ::NamedTuple, ::Any; kwargs...)
    error("The extension $(last(splitext(filename))) is currently not available. " *
          "A required package (e.g. JLD2 or WriteVTK) is not yet loaded.")
end


text:
module Test_ginv

using DescriptorSystems
using LinearAlgebra
using Polynomials
using Random
using Test

println("Test_ginv")
Random.seed!(2173)
@testset "ginv" begin

fast = true; Ty = Complex{Float64}; #Ty = Float64
disc = true;

# random examples
for Ty in (Float64, Complex{Float64})
      # for fast in (true, false)
      for disc in (false, true)

sys = rdss(0,0,0);
@time sysinv, info = ginv(sys);
@test iszero(sysinv) && info.nrank == 0 && info.nfp == 0

# invertible
sys = rss(4,3,3; T = Ty, disc);
@time sysinv, info = ginv(sys; fast);
@test iszero(sys*sysinv-I,atol=1.e-7) &&
      info.nrank == 3 && info.nfp == 0

sys = rdss(4,3,3; T = Ty, disc);
@time sysinv, info = ginv(sys);
@test iszero(sys*sysinv-I,atol=1.e-7)

# full row rank 1-2 inverse
sys = rss(2,2,3; T = Ty, disc);
@time sysinv, info = ginv(sys,atol=1.e-7);
@test iszero(sys*sysinv*sys-sys,atol=1.e-7) && iszero(sysinv*sys*sysinv-sysinv,atol=1.e-7) &&
      info.nrank == 2 && info.nfp == 2 && order(sysinv) == 2

sys = rss(2,2,3; T = Ty, disc);
@time sysinv, info = ginv(sys, mindeg = true, atol=1.e-7);
@test iszero(sys*sysinv*sys-sys,atol=1.e-7) && iszero(sysinv*sys*sysinv-sysinv,atol=1.e-7) &&
      info.nrank == 2 && info.nfp == 0 && order(sysinv) == 1

sys = rdss(2,2,3; T = Ty, disc);
disc ? poles = [0.5, 0.8] : poles = [-1,-2]
@time sysinv, info = ginv(sys; poles, atol=1.e-7);
@test iszero(sys*sysinv*sys-sys,atol=1.e-7) && iszero(sysinv*sys*sysinv-sysinv,atol=1.e-7) &&
      info.nrank == 2 && info.nfp == 2 && sort(real(gpole(sysinv))) ≈ sort(real(poles)) &&
      sort(abs.(gpole(sysinv))) ≈ sort(abs.(poles))

# full row rank 1-2-3-4 inverse
sys = rdss(2,2,3; T = Ty, disc);
@time sysinv, info = ginv(sys, type = "1-2-3-4", atol=1.e-7);
@test iszero(sys*sysinv*sys-sys,atol=1.e-7) && iszero(sysinv*sys*sysinv-sysinv,atol=1.e-7) &&
      iszero(sys*sysinv-(sys*sysinv)',atol=1.e-7) && iszero(sysinv*sys-(sysinv*sys)',atol=1.e-7) &&
      info.nrank == 2 && info.nfp == 0

sys = rss(0,2,3; T = Ty, disc);
@time sysinv, info = ginv(sys, type = "1-2-3-4", atol=1.e-7);
@test iszero(sys*sysinv*sys-sys,atol=1.e-7) && iszero(sysinv*sys*sysinv-sysinv,atol=1.e-7) &&
      iszero(sys*sysinv-(sys*sysinv)',atol=1.e-7) && iszero(sysinv*sys-(sysinv*sys)',atol=1.e-7) &&
      info.nrank == 2 && info.nfp == 0 && order(sysinv) == 0

# full column rank
sys = rss(2,3,2; T = Ty, disc);
@time sysinv, info = ginv(sys,atol=1.e-7);
@test iszero(sys*sysinv*sys-sys,atol=1.e-7) && iszero(sysinv*sys*sysinv-sysinv,atol=1.e-7) &&
      info.nrank == 2 && info.nfp == 2 && order(sysinv) == 2

sys = rss(0,3,2; T = Ty, disc);
@time sysinv, info = ginv(sys, type = "1-2-3-4", atol=1.e-7);
@test iszero(sys*sysinv*sys-sys,atol=1.e-7) && iszero(sysinv*sys*sysinv-sysinv,atol=1.e-7) &&
      iszero(sys*sysinv-(sys*sysinv)',atol=1.e-7) && iszero(sysinv*sys-(sysinv*sys)',atol=1.e-7) &&
      info.nrank == 2 && info.nfp == 0 && order(sysinv) == 0

# full column rank: minimum order
sys = rdss(2,3,2; T = Ty, disc);
@time sysinv, info = ginv(sys, mindeg = true, atol=1.e-7);
@test iszero(sys*sysinv*sys-sys,atol=1.e-7) && iszero(sysinv*sys*sysinv-sysinv,atol=1.e-7) &&
      info.nrank == 2 && info.nfp == 0 && order(sysinv) == 1

sys = rdss(2,3,2; T = Ty, disc);
disc ? poles = [0.5, 0.8] : poles = [-1,-2]
@time sysinv, info = ginv(sys; poles, atol=1.e-7);
@test iszero(sys*sysinv*sys-sys,atol=1.e-7) && iszero(sysinv*sys*sysinv-sysinv,atol=1.e-7) &&
      info.nrank == 2 && info.nfp == 2 && sort(real(gpole(sysinv))) ≈ sort(real(poles)) &&
      sort(abs.(gpole(sysinv))) ≈ sort(abs.(poles))

# full column rank 1-2-3-4 inverse
sys = rdss(2,3,2; T = Ty, disc);
@time sysinv, info = ginv(sys, type = "1-2-3-4", atol=1.e-7);
@test iszero(sys*sysinv*sys-sys,atol=1.e-7) && iszero(sysinv*sys*sysinv-sysinv,atol=1.e-7) &&
      iszero(sys*sysinv-(sys*sysinv)',atol=1.e-7) && iszero(sysinv*sys-(sysinv*sys)',atol=1.e-7) &&
      info.nrank == 2 && info.nfp == 0

# rank less than minimum dimension
sys = rss(2,3,2; T = Ty, disc)*rss(2,2,3; T = Ty, disc);
@time sysinv, info = ginv(sys,atol=1.e-7);
@test iszero(sys*sysinv*sys-sys,atol=1.e-7) && iszero(sysinv*sys*sysinv-sysinv,atol=1.e-7) &&
      info.nrank == 2 && info.nfp == 4 && order(sysinv) == 4

sys = rss(0,3,2; T = Ty, disc)*rss(0,2,3; T = Ty, disc);
@time sysinv, info = ginv(sys, type = "1-2", atol=1.e-7)
@test iszero(sys*sysinv*sys-sys,atol=1.e-7) && iszero(sysinv*sys*sysinv-sysinv,atol=1.e-7) &&
      iszero(sys*sysinv-(sys*sysinv)',atol=1.e-7) && iszero(sysinv*sys-(sysinv*sys)',atol=1.e-7) &&
      info.nrank == 2 && info.nfp == 0 && order(sysinv) == 0

sys = rss(0,3,2; T = Ty, disc)*rss(0,2,3; T = Ty, disc);
@time sysinv, info = ginv(sys, type = "1-2-3-4", atol=1.e-7)
@test iszero(sys*sysinv*sys-sys,atol=1.e-7) && iszero(sysinv*sys*sysinv-sysinv,atol=1.e-7) &&
      iszero(sys*sysinv-(sys*sysinv)',atol=1.e-7) && iszero(sysinv*sys-(sysinv*sys)',atol=1.e-7) &&
      info.nrank == 2 && info.nfp == 0 && order(sysinv) == 0

sys = rss(2,3,2; T = Ty, disc)*rss(2,2,3; T = Ty, disc);
@time sysinv, info = ginv(sys, mindeg = true, atol=1.e-7);
@test iszero(sys*sysinv*sys-sys,atol=1.e-7) && iszero(sysinv*sys*sysinv-sysinv,atol=1.e-7) &&
      info.nrank == 2 && info.nfp == 0 && order(sysinv) == 2

sys = rdss(2,3,2; T = Ty, disc)*rss(2,2,3; T = Ty, disc);
disc ? poles = [0.5, 0.8] : poles = [-1,-2]
@time sysinv, info = ginv(sys; poles, atol=1.e-7);
@test iszero(sys*sysinv*sys-sys,atol=1.e-7) && iszero(sysinv*sys*sysinv-sysinv,atol=1.e-7) &&
      info.nrank == 2 && info.nfp == 4 &&
      minimum(abs.(gpole(sysinv).-poles[1])) < 1.e-5 && minimum(abs.(gpole(sysinv).-poles[2])) < 1.e-5

sys = rss(2,3,2; T = Ty, disc)*rdss(2,2,3; T = Ty, disc);
disc ? sdeg = 0.8 : sdeg = -1
@time sysinv, info = ginv(sys; sdeg, atol=1.e-7);
@test iszero(sys*sysinv*sys-sys,atol=1.e-7) && iszero(sysinv*sys*sysinv-sysinv,atol=1.e-7) &&
      info.nrank == 2 && info.nfp == 4 &&
      maximum(real(gpole(sysinv)).-sdeg) < 1.e-4

# 1-2-3-4 inverse
sys = rss(2,3,2; T = Ty, disc)*rdss(2,2,3; T = Ty, disc);
@time sysinv, info = ginv(sys, type = "1-2-3-4", atol=1.e-7);
@test iszero(sys*sysinv*sys-sys,atol=1.e-7) && iszero(sysinv*sys*sysinv-sysinv,atol=1.e-7) &&
      iszero(sys*sysinv-(sys*sysinv)',atol=1.e-7) && iszero(sysinv*sys-(sysinv*sys)',atol=1.e-7) &&
      info.nrank == 2 && info.nfp == 0

# 1-2-3 inverse
sys = rss(2,3,2; T = Ty, disc)*rdss(2,2,3; T = Ty, disc);
@time sysinv, info = ginv(sys, type = "1-2-3", atol=1.e-7);
@test iszero(sys*sysinv*sys-sys,atol=1.e-7) && iszero(sysinv*sys*sysinv-sysinv,atol=1.e-7) &&
      iszero(sys*sysinv-(sys*sysinv)',atol=1.e-7) &&
      info.nrank == 2 && info.nfp == 2

# 1-2-4 inverse
sys = rss(2,3,2; T = Ty, disc)*rdss(2,2,3; T = Ty, disc);
@time sysinv, info = ginv(sys, type = "1-2-4", atol=1.e-7);
@test iszero(sys*sysinv*sys-sys,atol=1.e-7) && iszero(sysinv*sys*sysinv-sysinv,atol=1.e-7) &&
       iszero(sysinv*sys-(sysinv*sys)',atol=1.e-7) &&
      info.nrank == 2 && info.nfp == 2

sys = rss(0,3,2; T = Ty, disc)*rss(0,2,3; T = Ty, disc);
@time sysinv, info = ginv(sys, type = "1-2-3-4", atol=1.e-7);
@test iszero(sys*sysinv*sys-sys,atol=1.e-7) && iszero(sysinv*sys*sysinv-sysinv,atol=1.e-7) &&
      iszero(sys*sysinv-(sys*sysinv)',atol=1.e-7) && iszero(sysinv*sys-(sysinv*sys)',atol=1.e-7) &&
      info.nrank == 2 && info.nfp == 0 && order(sysinv) == 0

sys = rss(0,3,2; T = Ty, disc)*rss(0,2,3; T = Ty, disc);
@time sysinv, info = ginv(sys, type = "1-2-3", atol=1.e-7);
@test iszero(sys*sysinv*sys-sys,atol=1.e-7) && iszero(sysinv*sys*sysinv-sysinv,atol=1.e-7) &&
      iszero(sys*sysinv-(sys*sysinv)',atol=1.e-7) &&
      info.nrank == 2 && info.nfp == 0 && order(sysinv) == 0

sys = rss(0,3,2; T = Ty, disc)*rss(0,2,3; T = Ty, disc);
@time sysinv, info = ginv(sys, type = "1-2-4", atol=1.e-7);
@test iszero(sys*sysinv*sys-sys,atol=1.e-7) && iszero(sysinv*sys*sysinv-sysinv,atol=1.e-7) &&
      iszero(sysinv*sys-(sysinv*sys)',atol=1.e-7) &&
      info.nrank == 2 && info.nfp == 0 && order(sysinv) == 0

sys = rss(0,3,2; T = Ty, disc)*rss(0,2,3; T = Ty, disc);
@time sysinv, info = ginv(sys, type = "1-2-4", atol=1.e-7);
@test iszero(sys*sysinv*sys-sys,atol=1.e-7) && iszero(sysinv*sys*sysinv-sysinv,atol=1.e-7) &&
      iszero(sysinv*sys-(sysinv*sys)',atol=1.e-7) &&
      info.nrank == 2 && info.nfp == 0 && order(sysinv) == 0

@test_throws  ErrorException ginv(sys, atol=1.e-7, type = "1")

end # disc
# end # fast
end # Ty
end # ginv

end # module


text:
mutable struct RateTimer
    tlast::Float64
    elapsed::Float64
    rate::Float64
    paused::Bool
    tpaused::Float64
    function RateTimer(rate)
        t = time_ns()
        new(t, 0, rate, true, t)
    end
end
RateTimer() = RateTimer(1)

function Base.time_ns(rt::RateTimer)
    if !rt.paused
        tnow = time_ns()
        elapsed = tnow - rt.tlast
        rt.tlast = tnow
        rt.elapsed += elapsed * rt.rate
    end
    return rt.elapsed
end

Base.time(rt::RateTimer) = time_ns(rt) / 1e9
stop!(rt::RateTimer) = (rt.elapsed = time_ns(rt); rt.paused = true; rt)
start!(rt::RateTimer) = (rt.tlast = time_ns(); rt.paused = false; rt)
setrate!(rt::RateTimer, r) = (rt.rate = r; rt)
LyceumBase.reset!(rt::RateTimer) = (rt.tlast = time_ns(); rt.elapsed = 0; rt)


text:
@coeff function solve_steady!(reactor::Reactor)
  LHS = B_0

  RHS = 1 / sqrt(
    R_0^1.5 * K_BS + R_0^2 * K_CD
  )

  equation = LHS - RHS

  try
    work_R_0 = find_zero(equation, (R_0_min, R_0_max))

    if is_valid(reactor, work_R_0, B_0)
      reactor.R_0 = work_R_0
    else
      reactor.R_0 = NaN
    end
  catch
    reactor.R_0 = NaN
  end
end


text:
# export Circular, Circle

# """
# $(TYPEDEF)

# Circular is a `Manifolds.Circle{ℝ}` mechanization of one rotation, with `theta in [-pi,pi)`.
# """
# @defVariable Circular Circle() [0.0;]

# Base.convert(::Type{<:MB.AbstractManifold}, ::InstanceType{Circular}) = Manifolds.Circle()

#


text:
@use GADT

import Base: convert

struct Fun{T,R}
    fn::Function
end

function (typed_fn::Fun{T,R})(arg::T)::R where {T,R}
    typed_fn.fn(arg)
end

function convert(::Type{Fun{T,R}}, fn::Function) where {T,R}
    Fun{T,R}(fn)
end

function convert(::Type{Fun{T,R}}, fn::Fun{C,D}) where {T,R,C<:T,D<:R}
    Fun{T,R}(fn.fn)
end

⇒(::Type{A}, ::Type{B}) where {A,B} = Fun{A,B}

@data public Exp{T} begin
    Sym{T}::Symbol => Exp{T}
    Val{A}::A => Exp{A}
    App{A,B}::(Exp{Fun{A,B}}, Exp{<:A}) => Exp{B}
    Lam{A,B}::(Symbol, Exp{B}) => Exp{Fun{A,B}}
    If{A}::(Exp{Bool}, Exp{A}, Exp{A}) => Exp{A}
end

function substitute(template::Exp{T}, pair::Tuple{Symbol,Exp{G}}) where {T,G<:T}
    (sym, exp) = pair
    @match template begin
        Sym(&sym) => exp
        Val(_) => template
        App(f, a) => App(substitute(f, pair), substitute(a, pair))
        Lam(&sym, exp) => template
        Lam(sym, exp) => Lam(sym, substitute(exp, pair))
        If(cond, exp1, exp2) => begin
            (cond, exp1, exp2) =
                substitute(cond, pair), substitute(exp1, pair), substitute(exp2, pair)
            If(cond, exp1, exp2)
        end
    end
end

function eval_exp(exp::Exp{T_}, ctx::Dict{Symbol,Any}) where {T_}
    let T = T_ # fix static parameter issues for Julia 1.0 and 1.1
    @match exp begin
        Sym(a) => (ctx[a]::T, ctx)
        Val(a::T) => (a, ctx)
        App{A,T}(f::Exp{Fun{A,T}}, arg::Exp{<:A}) where {A, T} =>
            let (f, ctx) = eval_exp(f, ctx), (arg, ctx) = eval_exp(arg, ctx)
                (f(arg), ctx)
            end
        Lam{A,B}(sym, exp::Exp{B}) where {A,B} =>
            let f(x) = begin
                    eval_exp(substitute(exp, sym => Val(x)), ctx)[1]
                end
                (f, ctx)
            end
        If(cond, exp1, exp2) => let (cond, ctx) = eval_exp(cond, ctx)
            eval_exp(cond ? exp1 : exp2, ctx)
        end
        _ => error(exp)
    end
    end
end

add = Val{Number ⇒ Number ⇒ Number}(x -> y -> x + y)
sub = Val{Number ⇒ Number ⇒ Number}(x -> y -> x - y)
gt = Val{Number ⇒ Number ⇒ Bool}(x -> y -> x > y)
ctx = Dict{Symbol,Any}()

@assert 3 == eval_exp(App(App(add, Val(1)), Val(2)), ctx)[1]
@assert -1 == eval_exp(App(App(sub, Val(1)), Val(2)), ctx)[1]
@assert 1 == eval_exp(
    If(
        App(App(gt, Sym{Int}(:x)), Sym{Int}(:y)),
        App(App(sub, Sym{Int}(:x)), Sym{Int}(:y)),
        App(App(sub, Sym{Int}(:y)), Sym{Int}(:x)),
    ),
    Dict{Symbol,Any}(:x => 1, :y => 2),
)[1]


text:
using DrWatson
@quickactivate "ChaosThroughBilliards"
include(srcdir("style.jl"))

ps = 2.0

bdsta = billiard_stadium(1.0, 1.0)

N = 300

billiard_video(
    videodir("intro$N.mp4"), bdsta, 1.0, 0.8, 0;
    frames = 1000, colors = COLORS[1:2],
    N = N, dx = 0.01,
    plot_particles = false, tailwidth = 1, particle_size = ps, res = MAXRES,
    dt = 0.0001, speed = 250, tail = 20000, # this makes ultra fine temporal resolution
)

println("Done.")

max_stars_repo_name: JuliaDynamics/ChaosThroughBilliards

max_issues_repo_name: JuliaDynamics/ChaosThroughBilliards

max_forks_repo_name: JuliaDynamics/ChaosThroughBilliards


text:
include("header.jl")
include("simulator.jl")
include("visual.jl")
include("updateyaml.jl")

# include("naturalgas.jl")
gas = CreateSolution(mech);
ns = gas.n_species;
nr = gas.n_reactions;

include("expdataset.jl")

include("sensBVP.jl")

npr = 16;
p = zeros(npr * 3);
include("sensBVPpart.jl")

regression_plot(; max=5)

opt = ADAMW(1.e-3, (0.9, 0.999), 1.e-6);

grad_max = 10^(2.5);

include("callback.jl")

ind_sl = 17:nr

epochs = ProgressBar(iter:1000);
l_epoch = ones(n_exp);
grad_norm = ones(n_exp);
for epoch in epochs
    global p
    l_epoch .= 1.0
    grad_norm .= 1.0
    for i_exp in randperm(n_exp)
        T0 = conds[i_exp, 1]
        P = conds[i_exp, 2] * one_atm
        phi = conds[i_exp, 3]
        idt0 = conds[i_exp, 4]
        ts, pred = get_idt(T0, P, phi, p; dT=dT, tfinal=10.0)
        ngfull = length(ts)
        idt = ts[end]
        l_epoch[i_exp] = (log(idt / idt0))^2
        @printf("%d ng: %d (%d) idt %.2e idt0 %.2e \n",
                i_exp, length(ts), ngfull, idt, idt0)

        if (idt < 9.9) & (length(ts) < 2000) & (i_exp in l_exp_train)
            grad = 2 * log(idt / idt0) * sensBVP_mthreadp(ts, pred, p)
            grad_norm[i_exp] = norm(grad, 2)
            if grad_norm[i_exp] > grad_max
                @. grad = grad / grad_norm[i_exp] * grad_max
            end
            # grad[ind_sl] .= 0.0
            if (iter > 1)
                update!(opt, p, grad)
            end
        end
    end
    loss = mean(l_epoch[l_exp_train])
    loss_val = mean(l_epoch[l_exp_val])
    g_norm = mean(grad_norm[l_exp_train])
    set_description(
        epochs,
        string(
            @sprintf("loss %.2e val %.2e pnorm %.2f gnorm %.2f",
                loss, loss_val, norm(p), g_norm
            )
        ),
    )
    cb(p, loss, loss_val, g_norm; doplot=true)
end

# @load string(ckpt_path, "/../p.bson") p

@save string(ckpt_path, "/../p.bson") p

regression_plot(; max=100)

max_stars_repo_name: DENG-MIT/DeepReduction

max_issues_repo_name: DENG-MIT/DeepReduction

max_forks_repo_name: DENG-MIT/DeepReduction