text:
using Combinatorics
using ITensors
using Test

include("util.jl")

@testset "MPS Basics" begin

  N = 10
  sites = [Index(2,"Site") for n=1:N]
  psi = MPS(sites)
  @test length(psi) == N
  @test length(MPS()) == 0
  @test isnothing(flux(psi))

  str = split(sprint(show, psi), '\n')
  @test str[1] == "MPS"
  @test length(str) == length(psi) + 2

  @test siteind(psi,2) == sites[2]
  @test findfirstsiteind(psi, sites[2]) == 2
  @test findfirstsiteind(psi, sites[4]) == 4
  @test findfirstsiteinds(psi, IndexSet(sites[5])) == 5
  @test hasind(psi[3],linkind(psi,2))
  @test hasind(psi[3],linkind(psi,3))

  @test isnothing(linkind(psi, N))
  @test isnothing(linkind(psi, N+1))
  @test isnothing(linkind(psi, 0))
  @test isnothing(linkind(psi, -1))
  @test linkind(psi, 3) == commonind(psi[3], psi[4])

  psi[1] = ITensor(sites[1])
  @test hasind(psi[1],sites[1])

  @testset "N=1 MPS" begin
    sites1 = [Index(2,"Site,n=1")]
    psi = MPS(sites1)
    @test length(psi) == 1
    @test siteind(psi,1) == sites1[1]
    @test siteinds(psi)[1] == sites1[1]
  end

  @testset "Missing links" begin
    psi = MPS([randomITensor(sites[i]) for i in 1:N])
    @test isnothing(linkind(psi, 1))
    @test isnothing(linkind(psi, 5))
    @test isnothing(linkind(psi, N))
    @test maxlinkdim(psi) == 0
    @test psi ⋅ psi ≈ *(dag(psi)..., psi...)[]
  end

  @testset "productMPS" begin

    @testset "vector of string input" begin
      sites = siteinds("S=1/2", N)
      state = fill("", N)
      for j=1:N
        state[j] = isodd(j) ? "Up" : "Dn"
      end
      psi = productMPS(sites,state)
      for j=1:N
        sign = isodd(j) ? +1.0 : -1.0
        @test (psi[j] * op(sites,"Sz",j) * dag(prime(psi[j],"Site")))[] ≈ sign/2
      end
      @test_throws DimensionMismatch productMPS(sites, fill("", N - 1))
    end

    @testset "String input" begin
      sites = siteinds("S=1/2", N)
      psi = productMPS(sites, "Dn")
      for j=1:N
        sign = -1.0
        @test (psi[j] * op(sites,"Sz",j) * dag(prime(psi[j],"Site")))[] ≈ sign/2
      end
    end

    @testset "Int input" begin
      sites = siteinds("S=1/2", N)
      psi = productMPS(sites, 2)
      for j=1:N
        sign = -1.0
        @test (psi[j] * op(sites,"Sz",j) * dag(prime(psi[j],"Site")))[] ≈ sign/2
      end
    end

    @testset "vector of int input" begin
      sites = siteinds("S=1/2",N)
      state = fill(0,N)
      for j=1:N
        state[j] = isodd(j) ? 1 : 2
      end
      psi = productMPS(sites,state)
      for j=1:N
        sign = isodd(j) ? +1.0 : -1.0
        @test (psi[j]*op(sites,"Sz",j)*dag(prime(psi[j],"Site")))[] ≈ sign/2
      end
    end

    @testset "vector of ivals input" begin
      sites  = siteinds("S=1/2",N)
      states = fill(0,N)
      for j=1:N
        states[j] = isodd(j) ? 1 : 2
      end
      ivals  = [state(sites[n],states[n]) for n=1:length(sites)]
      psi = productMPS(ivals)
      for j=1:N
        sign = isodd(j) ? +1.0 : -1.0
        @test (psi[j]*op(sites,"Sz",j)*dag(prime(psi[j],"Site")))[] ≈ sign/2
      end

      @testset "ComplexF64 eltype" begin
        sites  = siteinds("S=1/2",N)
        psi = productMPS(ComplexF64,sites,fill(1,N))
        for j=1:N
          @test eltype(psi[j]) <: ComplexF64
        end
      end
    end

    @testset "N=1 case" begin
      site = Index(2,"Site,n=1")
      psi = productMPS([site],[1])
      @test psi[1][1] ≈ 1.0
      @test psi[1][2] ≈ 0.0
      psi = productMPS([site],[2])
      @test psi[1][1] ≈ 0.0
      @test psi[1][2] ≈ 1.0
    end
  end

  @testset "randomMPS with chi==1" begin
    phi = randomMPS(sites)

    @test maxlinkdim(phi) == 1

    @test hasind(phi[1],sites[1])
    @test norm(phi[1])≈1.0

    @test hasind(phi[4],sites[4])
    @test norm(phi[4])≈1.0
  end

  @testset "inner different MPS" begin
    phi = randomMPS(sites)
    psi = randomMPS(sites)
    phipsi = dag(phi[1])*psi[1]
    for j = 2:N
      phipsi *= dag(phi[j])*psi[j]
    end
    @test phipsi[] ≈ inner(phi,psi)

    badsites = [Index(2) for n=1:N+1]
    badpsi = randomMPS(badsites)
    @test_throws DimensionMismatch inner(phi,badpsi)
  end

  @testset "inner same MPS" begin
    psi = randomMPS(sites)
    psidag = dag(psi)
    #ITensors.prime_linkinds!(psidag)
    psipsi = psidag[1]*psi[1]
    for j = 2:N
      psipsi *= psidag[j]*psi[j]
    end
    @test psipsi[] ≈ inner(psi,psi)
  end

  @testset "norm MPS" begin
    psi = randomMPS(sites,10)
    psidag = ITensors.sim_linkinds(dag(psi))
    psi² = ITensor(1)
    for j = 1:N
      psi² *= psidag[j] * psi[j]
    end
    @test psi²[] ≈ psi ⋅ psi
    @test sqrt(psi²[]) ≈ norm(psi)
    for j in 1:N
      psi[j] .*= j
    end
    @test norm(psi) ≈ factorial(N)
  end

  @testset "lognorm MPS" begin
    psi = randomMPS(sites,10)
    for j in 1:N
      psi[j] .*= j
    end
    psidag = ITensors.sim_linkinds(dag(psi))
    psi² = ITensor(1)
    for j = 1:N
      psi² *= psidag[j] * psi[j]
    end
    @test psi²[] ≈ psi ⋅ psi
    @test 0.5 * log(psi²[]) ≈ lognorm(psi)
    @test lognorm(psi) ≈ log(factorial(N))
  end

  @testset "scaling MPS" begin
    psi = randomMPS(sites)
    twopsidag = 2.0*dag(psi)
    #ITensors.prime_linkinds!(twopsidag)
    @test inner(twopsidag, psi) ≈ 2.0*inner(psi,psi)
  end

  @testset "flip sign of MPS" begin
    psi = randomMPS(sites)
    minuspsidag = -dag(psi)
    #ITensors.primelinkinds!(minuspsidag)
    @test inner(minuspsidag, psi) ≈ -inner(psi,psi)
  end

  @testset "add MPS" begin
    psi = randomMPS(sites)
    phi = deepcopy(psi)
    xi = add(psi, phi)
    @test inner(xi, xi) ≈ 4.0 * inner(psi, psi)
    # sum of many MPSs
    Ks = [randomMPS(sites) for i in 1:3]
    K12  = add(Ks[1], Ks[2])
    K123 = add(K12, Ks[3])
    @test inner(sum(Ks), K123) ≈ inner(K123,K123)
  end

  @testset "+ MPS" begin
    psi = randomMPS(sites)
    phi = deepcopy(psi)
    xi = psi + phi
    @test inner(xi, xi) ≈ 4.0 * inner(psi, psi)
    # sum of many MPSs
    Ks = [randomMPS(sites) for i in 1:3]
    K12  = Ks[1] + Ks[2]
    K123 = K12 + Ks[3]
    @test inner(sum(Ks), K123) ≈ inner(K123,K123)
  end

  sites = siteinds(2,N)
  psi = MPS(sites)
  @test length(psi) == N # just make sure this works
  @test length(siteinds(psi)) == N

  psi = randomMPS(sites)
  l0s = linkinds(psi)
  orthogonalize!(psi, N-1)
  ls = linkinds(psi)
  for (l0,l) in zip(l0s,ls)
    @test tags(l0) == tags(l)
  end
  @test ITensors.leftlim(psi) == N-2
  @test ITensors.rightlim(psi) == N
  orthogonalize!(psi, 2)
  @test ITensors.leftlim(psi) == 1
  @test ITensors.rightlim(psi) == 3
  psi = randomMPS(sites)
  ITensors.setrightlim!(psi, N+1) # do this to test qr
                                  # from rightmost tensor
  orthogonalize!(psi, div(N, 2))
  @test ITensors.leftlim(psi) == div(N, 2) - 1
  @test ITensors.rightlim(psi) == div(N, 2) + 1

  @test isnothing(linkind(MPS(fill(ITensor(), N), 0, N + 1), 1))

  @testset "replacebond!" begin
  # make sure factorization preserves the bond index tags
    psi = randomMPS(sites)
    phi = psi[1]*psi[2]
    bondindtags = tags(linkind(psi,1))
    replacebond!(psi,1,phi)
    @test tags(linkind(psi,1)) == bondindtags

    # check that replacebond! updates llim and rlim properly
    orthogonalize!(psi,5)
    phi = psi[5]*psi[6]
    replacebond!(psi, 5, phi; ortho = "left")
    @test ITensors.leftlim(psi) == 5
    @test ITensors.rightlim(psi) == 7

    phi = psi[5]*psi[6]
    replacebond!(psi, 5, phi; ortho = "right")
    @test ITensors.leftlim(psi) == 4
    @test ITensors.rightlim(psi)==6

    ITensors.setleftlim!(psi, 3)
    ITensors.setrightlim!(psi, 7)
    phi = psi[5]*psi[6]
    replacebond!(psi, 5, phi; ortho = "left")
    @test ITensors.leftlim(psi) == 3
    @test ITensors.rightlim(psi) == 7
  end

end

@testset "orthogonalize! with QNs" begin
  N = 8
  sites = siteinds("S=1/2",N, conserve_qns=true)
  init_state = [isodd(n) ? "Up" : "Dn" for n=1:N]
  psi0 = productMPS(sites,init_state)
  orthogonalize!(psi0,4)
  @test ITensors.leftlim(psi0) == 3
  @test ITensors.rightlim(psi0) == 5
end

# Helper function for making MPS
function basicRandomMPS(N::Int;dim=4)
  sites = [Index(2,"Site") for n=1:N]
  M = MPS(sites)
  links = [Index(dim,"n=$(n-1),Link") for n=1:N+1]
  for n=1:N
    M[n] = randomITensor(links[n],sites[n],links[n+1])
  end
  M[1] *= delta(links[1])
  M[N] *= delta(links[N+1])
  M[1] /= sqrt(inner(M,M))
  return M
end

@testset "MPS gauging and truncation" begin

  N = 30

  @testset "orthogonalize! method" begin
    c = 12
    M = basicRandomMPS(N)
    orthogonalize!(M,c)

    @test ITensors.leftlim(M) == c-1
    @test ITensors.rightlim(M) == c+1

    # Test for left-orthogonality
    L = M[1]*prime(M[1],"Link")
    l = linkind(M,1)
    @test norm(L-delta(l,l')) < 1E-12
    for j=2:c-1
      L = L*M[j]*prime(M[j],"Link")
      l = linkind(M,j)
      @test norm(L-delta(l,l')) < 1E-12
    end

    # Test for right-orthogonality
    R = M[N]*prime(M[N],"Link")
    r = linkind(M,N-1)
    @test norm(R-delta(r,r')) < 1E-12
    for j in reverse(c+1:N-1)
      R = R*M[j]*prime(M[j],"Link")
      r = linkind(M,j-1)
      @test norm(R-delta(r,r')) < 1E-12
    end

    @test norm(M[c]) ≈ 1.0
  end

  @testset "truncate! method" begin
    M = basicRandomMPS(N;dim=10)
    M0 = copy(M)
    truncate!(M;maxdim=5)

    @test ITensors.rightlim(M) == 2

    # Test for right-orthogonality
    R = M[N]*prime(M[N],"Link")
    r = linkind(M,N-1)
    @test norm(R-delta(r,r')) < 1E-12
    for j in reverse(2:N-1)
      R = R*M[j]*prime(M[j],"Link")
      r = linkind(M,j-1)
      @test norm(R-delta(r,r')) < 1E-12
    end

    @test inner(M,M0) > 0.1
  end

end

@testset "Other MPS methods" begin

  @testset "sample! method" begin
    N = 10
    sites = [Index(3,"Site,n=$n") for n=1:N]
    psi = makeRandomMPS(sites,chi=3)
    nrm2 = inner(psi,psi)
    psi[1] *= (1.0/sqrt(nrm2))

    s = sample!(psi)

    @test length(s) == N
    for n=1:N
      @test 1 <= s[n] <= 3
    end

    # Throws becase not orthogonalized to site 1:
    orthogonalize!(psi,3)
    @test_throws ErrorException sample(psi)

    # Throws becase not normalized
    orthogonalize!(psi,1)
    psi[1] *= (5.0/norm(psi[1]))
    @test_throws ErrorException sample(psi)

    # Works when ortho & normalized:
    orthogonalize!(psi,1)
    psi[1] *= (1.0/norm(psi[1]))
    s = sample(psi)
    @test length(s) == N
  end

  @testset "randomMPS with chi > 1" begin
    N = 20
    chi = 8
    sites = siteinds(2,N)
    M = randomMPS(sites,chi)

    @test ITensors.leftlim(M) == 0
    @test ITensors.rightlim(M) == 2

    @test norm(M[1]) ≈ 1.0

    @test maxlinkdim(M) == chi

    # Test for right-orthogonality
    R = M[N]*prime(M[N],"Link")
    r = linkind(M,N-1)
    @test norm(R-delta(r,r')) < 1E-10
    for j in reverse(2:N-1)
      R = R*M[j]*prime(M[j],"Link")
      r = linkind(M,j-1)
      @test norm(R-delta(r,r')) < 1E-10
    end
  end

  @testset "randomMPS from initial state (QN case)" begin
    N = 20
    chi = 8
    sites = siteinds("S=1/2",N;conserve_qns=true)

    # Make flux-zero random MPS
    state = [isodd(n) ? 1 : 2 for n=1:N]
    M = randomMPS(sites,state,chi)
    @test flux(M) == QN("Sz",0)

    @test ITensors.leftlim(M) == 0
    @test ITensors.rightlim(M) == 2

    @test norm(M[1]) ≈ 1.0
    @test inner(M,M) ≈ 1.0

    @test maxlinkdim(M) == chi

    # Test making random MPS with different flux
    state[1] = 2
    M = randomMPS(sites,state,chi)
    @test flux(M) == QN("Sz",-2)
    state[3] = 2
    M = randomMPS(sites,state,chi)
    @test flux(M) == QN("Sz",-4)
  end

  @testset "swapbondsites" begin
    N = 5
    sites = siteinds("S=1/2", N)
    ψ0 = randomMPS(sites)
    ψ = replacebond(ψ0, 3, ψ0[3] * ψ0[4];
                    swapsites = true,
                    cutoff = 1e-15)
    @test siteind(ψ, 1) == siteind(ψ0, 1)
    @test siteind(ψ, 2) == siteind(ψ0, 2)
    @test siteind(ψ, 4) == siteind(ψ0, 3)
    @test siteind(ψ, 3) == siteind(ψ0, 4)
    @test siteind(ψ, 5) == siteind(ψ0, 5)
    @test prod(ψ) ≈ prod(ψ0)
    @test maxlinkdim(ψ) == 1

    ψ = swapbondsites(ψ0, 4;
                      cutoff = 1e-15)
    @test siteind(ψ, 1) == siteind(ψ0, 1)
    @test siteind(ψ, 2) == siteind(ψ0, 2)
    @test siteind(ψ, 3) == siteind(ψ0, 3)
    @test siteind(ψ, 5) == siteind(ψ0, 4)
    @test siteind(ψ, 4) == siteind(ψ0, 5)
    @test prod(ψ) ≈ prod(ψ0)
    @test maxlinkdim(ψ) == 1
  end

  @testset "map!" begin
    N = 5
    s = siteinds("S=½", N)
    M0 = productMPS(s, "↑")

    # Test map! with limits getting set
    M = orthogonalize(M0, 1)
    @test ITensors.leftlim(M) == 0
    @test ITensors.rightlim(M) == 2
    map!(prime, M)
    @test ITensors.leftlim(M) == 0
    @test ITensors.rightlim(M) == N+1

    # Test map! without limits getting set
    M = orthogonalize(M0, 1)
    map!(prime, M, set_limits = false)
    @test ITensors.leftlim(M) == 0
    @test ITensors.rightlim(M) == 2

    # Test prime! with limits getting set
    M = orthogonalize(M0, 1)
    @test ITensors.leftlim(M) == 0
    @test ITensors.rightlim(M) == 2
    prime!(M, set_limits = true)
    @test ITensors.leftlim(M) == 0
    @test ITensors.rightlim(M) == N+1

    # Test prime! without limits getting set
    M = orthogonalize(M0, 1)
    prime!(M)
    @test ITensors.leftlim(M) == 0
    @test ITensors.rightlim(M) == 2
  end

  @testset "setindex!(::MPS, _, ::Colon)" begin
    N = 4
    s = siteinds("S=½", N)
    ψ = randomMPS(s)
    ϕ = productMPS(s, "↑")
    orthogonalize!(ϕ, 1)
    ψ[:] = ϕ
    @test ITensors.orthocenter(ψ) == 1
    @test inner(ψ, ϕ) ≈ 1

    ψ = randomMPS(s)
    ϕ = productMPS(s, "↑")
    orthogonalize!(ϕ, 1)
    ψ[:] = ITensors.data(ϕ)
    @test ITensors.leftlim(ψ) == 0
    @test ITensors.rightlim(ψ) == N+1
    @test inner(ψ, ϕ) ≈ 1
  end

  @testset "findsite[s](::MPS/MPO, is)" begin
    s = siteinds("S=1/2", 5)
    ψ = randomMPS(s)
    l = linkinds(ψ)

    A = randomITensor(s[4]', s[2]', dag(s[4]), dag(s[2]))

    @test findsite(ψ, s[3]) == 3
    @test findsite(ψ, (s[3], s[5])) == 3
    @test findsite(ψ, l[2]) == 2
    @test findsite(ψ, A) == 2

    @test findsites(ψ, s[3]) == [3]
    @test findsites(ψ, (s[4], s[1])) == [1, 4]
    @test findsites(ψ, l[2]) == [2, 3]
    @test findsites(ψ, (l[2], l[3])) == [2, 3, 4]
    @test findsites(ψ, A) == [2, 4]

    M = randomMPO(s)
    lM = linkinds(M)

    @test findsite(M, s[4]) == 4
    @test findsite(M, s[4]') == 4
    @test findsite(M, (s[4]', s[4])) == 4
    @test findsite(M, (s[4]', s[3])) == 3
    @test findsite(M, lM[2]) == 2
    @test findsite(M, A) == 2

    @test findsites(M, s[4]) == [4]
    @test findsites(M, s[4]') == [4]
    @test findsites(M, (s[4]', s[4])) == [4]
    @test findsites(M, (s[4]', s[3])) == [3, 4]
    @test findsites(M, (lM[2], lM[3])) == [2, 3, 4]
    @test findsites(M, A) == [2, 4]
  end

  @testset "[first]siteind[s](::MPS/MPO, j::Int)" begin
    s = siteinds("S=1/2", 5)
    ψ = randomMPS(s)
    @test firstsiteind(ψ, 3) == s[3]
    @test siteind(ψ, 4) == s[4]
    @test isnothing(siteind(ψ, 4; plev = 1))
    @test siteinds(ψ, 3) == IndexSet(s[3])
    @test siteinds(ψ, 3; plev = 1) == IndexSet()

    M = randomMPO(s)
    @test noprime(firstsiteind(M, 4)) == s[4]
    @test firstsiteind(M, 4; plev = 0) == s[4]
    @test firstsiteind(M, 4; plev = 1) == s[4]'
    @test siteind(M, 4) == s[4]
    @test siteind(M, 4; plev = 0) == s[4]
    @test siteind(M, 4; plev = 1) == s[4]'
    @test isnothing(siteind(M, 4; plev = 2))
    @test siteinds(M, 3) == IndexSet(s[3], s[3]')
    @test siteinds(M, 3; plev = 1) == IndexSet(s[3]')
    @test siteinds(M, 3; plev = 0) == IndexSet(s[3])
    @test siteinds(M, 3; tags = "n=2") == IndexSet()
  end

  @testset "movesites $N sites" for N in 1:7
    s0 = siteinds("S=1/2", N)
    ψ0 = productMPS(s0, "↑")
    for perm in permutations(1:N)
      s = s0[perm]
      ψ = productMPS(s, rand(("↑", "↓"), N))
      ns′ = [findsite(ψ0, i) for i in s]
      @test ns′ == perm
      ψ′ = movesites(ψ, 1:N .=> ns′)
      for n in 1:N
        @test siteind(ψ0, n) == siteind(ψ′, n)
      end
      @test prod(ψ) ≈ prod(ψ′)
    end
  end

  @testset "Construct MPS from ITensor" begin
    N = 5
    s = siteinds("S=1/2", N)
    l = [Index(3, "left_$n") for n in 1:2]
    r = [Index(3, "right_$n") for n in 1:2]

    #
    # MPS
    #

    A = randomITensor(s...)
    ψ = MPS(A, s)
    @test prod(ψ) ≈ A
    @test ITensors.orthocenter(ψ) == N
    @test maxlinkdim(ψ) == 4

    ψ0 = productMPS(s, "↑")
    A = prod(ψ0)
    ψ = MPS(A, s; cutoff = 1e-15)
    @test prod(ψ) ≈ A
    @test ITensors.orthocenter(ψ) == N
    @test maxlinkdim(ψ) == 1

    ψ0 = randomMPS(s, 2)
    A = prod(ψ0)
    ψ = MPS(A, s; cutoff = 1e-15, orthocenter = 2)
    @test prod(ψ) ≈ A
    @test ITensors.orthocenter(ψ) == 2
    @test maxlinkdim(ψ) == 2

    A = randomITensor(s..., l[1], r[1])
    ψ = MPS(A, s, leftinds = l[1], orthocenter = 3)
    ls = linkinds(ψ)
    @test hassameinds(ψ[1], (l[1], s[1], ls[1]))
    @test hassameinds(ψ[N], (r[1], s[N], ls[N - 1]))
    @test prod(ψ) ≈ A
    @test ITensors.orthocenter(ψ) == 3
    @test maxlinkdim(ψ) == 12

    A = randomITensor(s..., l..., r...)
    ψ = MPS(A, s, leftinds = l)
    ls = linkinds(ψ)
    @test hassameinds(ψ[1], (l..., s[1], ls[1]))
    @test hassameinds(ψ[N], (r..., s[N], ls[N - 1]))
    @test prod(ψ) ≈ A
    @test ITensors.orthocenter(ψ) == N
    @test maxlinkdim(ψ) == 36
  end

  @testset "Set range of MPS tensors" begin
    N = 5
    s = siteinds("S=1/2", N)
    ψ0 = randomMPS(s, 3)

    ψ = orthogonalize(ψ0, 2)
    A = prod(ITensors.data(ψ)[2:N-1])
    randn!(A)
    ϕ = MPS(A, s[2:N-1], orthocenter = 1)
    ψ[2:N-1] = ϕ
    @test prod(ψ) ≈ ψ[1] * A * ψ[N]
    @test maxlinkdim(ψ) == 4
    @test ITensors.orthocenter(ψ) == 2

    ψ = orthogonalize(ψ0, 1)
    A = prod(ITensors.data(ψ)[2:N-1])
    randn!(A)
    @test_throws AssertionError ψ[2:N-1] = A

    ψ = orthogonalize(ψ0, 2)
    A = prod(ITensors.data(ψ)[2:N-1])
    randn!(A)
    ψ[2:N-1, orthocenter = 3] = A
    @test prod(ψ) ≈ ψ[1] * A * ψ[N]
    @test maxlinkdim(ψ) == 4
    @test ITensors.orthocenter(ψ) == 3
  end

  @testset "movesites reverse sites" begin
    N = 6
    s = siteinds("S=1/2", N)
    ψ0 = randomMPS(s)
    ψ = movesites(ψ0, 1:N .=> reverse(1:N))
    for n in 1:N
      @test siteind(ψ, n) == s[N-n+1]
    end
  end

end

nothing


text:
# [![Visitors](https://hits.seeyoufarm.com/api/count/incr/badge.svg?url=https%3A%2F%2Fgithub.com%2FRoyiAvital%2FStackExchangeCodes&count_bg=%2379C83D&title_bg=%23555555&icon=&icon_color=%23E7E7E7&title=Visitors+%28Daily+%2F+Total%29&edge_flat=false)](https://github.com/RoyiAvital/Julia100Exercises)
#
# # 100 Julia Exercises with Solutions
#
# A set of introductory exercises for Julia. Based on [100 NumPy Exercises](https://github.com/rougier/numpy-100).
#
# In order to generate this file:
# 1. Clone the repository (Or download).
# 2. Activate the project.
# 3. Run:
# ```Julia
# using Literate;
# Literate.markdown("Julia100Exercises.jl", name = "README", execute = true, flavor = Literate.CommonMarkFlavor());
# ```
#
# **Remark**: Tested with Julia `1.7.2`.
#
# **To Do**:
# 1. Reevaluate the difficulty level of each question.

using Literate;
using LinearAlgebra;
using Statistics;
using Dates;
using DelimitedFiles;
using UnicodePlots;
using Random;
using Tullio;
using StaticKernels;

# ## Question 001
# Import the `LinearAlgebra` package under the name `LA`. (★☆☆)

import LinearAlgebra as LA;

# ## Question 002
# Print the version of Julia. (★☆☆)

println(VERSION);

# ## Question 003
# Create a non initialized vector of size 10 of `Float64`. (★☆☆)

vA = Vector{Float64}(undef, 10)

# Which is equivalent of

vA = Array{Float64, 1}(undef, 10)

# ## Question 004
# Find the memory size of any array. (★☆☆)

sizeof(vA)

# ## Question 005
# Show the documentation of the `+` (Add) method. (★☆☆)

@doc +

# ## Question 006
# Create a vector of zeros of size 10 but the fifth value which is 1. (★☆☆)

vA = zeros(10);
vA[5] = 1.0;
vA

# ## Question 007
# Create a vector with values ranging from 7 to 12. (★☆☆)

vA = 7:12

# The above is efficient type. In order to explicitly create a vector:

vA = collect(7:12)

# ## Question 008
# Reverse a vector (first element becomes last). (★☆☆)

vA = collect(1:3);
vB = vA[end:-1:1];
vB

# Alternative 001:

vB = reverse(vA);

# Alternative 002 (In place):

reverse!(vA);
vA

# ## Question 009
# Create a `3x3` matrix with values ranging from 0 to 8. (★☆☆)

mA = reshape(0:8, 3, 3)

# Another way would be:

mA = Matrix{Float64}(undef, 3, 3);
mA[:] = 0:8;

# ## Question 010
# Find indices of non zero elements from `[1, 2, 0, 0, 4, 0]`. (★☆☆)

findall(!iszero, [1, 2, 0, 0, 4, 0])

# ## Question 011
# Create a 3x3 identity matrix. (★☆☆)

mA = I(3)

# An alternative method (Explicit matrix) would be:

mA = Matrix(I, 3, 3) #<! For Float64: Matrix{Float64}(I, 3, 3)

# ## Question 012
# Create a `2x2x2` array with random values. (★☆☆)

mA = randn(2, 2, 2)

# ## Question 013
# Create a `5x5` array with random values and find the minimum and maximum values. (★☆☆)

mA = rand(5, 5);
minVal = minimum(mA)
#+
maxVal = maximum(mA)

# Using `extrema()` one could get both values at once:

minVal, maxVal = extrema(mA);

# ## Question 014
# Create a random vector of size 30 and find the mean value. (★☆☆)

meanVal = mean(randn(30))

# ## Question 015
# Create a 2d array with 1 on the border and 0 inside. (★☆☆)

mA = zeros(4, 4);
mA[:, [1, end]] .= 1;
mA[[1, end], :] .= 1;
mA

# An alternative way:

mA = ones(4, 5);
mA[2:(end - 1), 2:(end - 1)] .= 0;

# Using one line code:

mA = zeros(4, 5);
mA[[LinearIndices(mA)[cartIdx] for cartIdx in CartesianIndices(mA) if (any(cartIdx.I .== 1) || cartIdx.I[1] == size(mA, 1) || cartIdx.I[2] == size(mA, 2))]] .= 1;

# By Tomer Arnon (https://github.com/tomerarnon):

numRows = 4;
numCols = 5;
mA = Int[ii ∈ (1, numRows) || jj ∈ (1, numCols) for ii in 1:numRows, jj in 1:numCols];

# ## Question 016
# Add a border of zeros around the array. (★☆☆)

mB = zeros(size(mA) .+ 2);
mB[2:(end - 1), 2:(end - 1)] = mA;
mB

# ## Question 017
# Evaluate the following expressions. (★☆☆)

0 * NaN
#+
NaN == NaN
#+
Inf > NaN
#+
NaN - NaN
#+
NaN in [NaN]
#+
0.3 == 3 * 0.1

# ## Question 018
# Create a `5x5` matrix with values `[1, 2, 3, 4]` just below the diagonal. (★☆☆)

mA = diagm(5, 5, -1 => 1:4)

# ## Question 019
# Create a `8x8` matrix and fill it with a checkerboard pattern. (★☆☆)

mA = zeros(8, 8);
mA[2:2:end, 1:2:end] .= 1;
mA[1:2:end, 2:2:end] .= 1;
mA

# By Tomer Arnon (https://github.com/tomerarnon):

mA = Int[isodd(ii + jj) for ii in 1:8, jj in 1:8];

# ## Question 020
# Convert the linear index 100 to a _Cartesian Index_ of a size `(6,7,8)`. (★☆☆)

mA = rand(6, 7, 8);
cartIdx = CartesianIndices(mA)[100]; #<! See https://discourse.julialang.org/t/14666
mA[cartIdx] == mA[100]

# ## Question 021
# Create a checkerboard `8x8` matrix using the `repeat()` function. (★☆☆)

mA = repeat([0 1; 1 0], 4, 4)

# ## Question 022
# Normalize a `4x4` random matrix. (★☆☆)

mA = rand(4, 4);
mA .= (mA .- mean(mA)) ./ std(mA) #<! Pay attention that `@.` will yield error (`std()` and `mean()`)

# ## Question 023
# Create a custom type that describes a color as four unsigned bytes (`RGBA`). (★☆☆)

struct sColor
    R::UInt8;
    G::UInt8;
    B::UInt8;
    A::UInt8;
end

sMyColor = sColor(rand(UInt8, 4)...)

# ## Question 024
# Multiply a `2x4` matrix by a `4x3` matrix. (★☆☆)

mA = rand(2, 4) * randn(4, 3)

# ## Question 025
# Given a 1D array, negate all elements which are between 3 and 8, in place. (★☆☆)

vA = rand(1:10, 8);
map!(x -> ((x > 3) && (x < 8)) ? -x : x, vA, vA)

# Julia allows Math like notation as well (See `Q0027`):

vA = rand(1:10, 8);
map!(x -> 3 < x < 8 ? -x : x, vA, vA)

# Using logical indices one could use:

vA[3 .< vA .< 8] .*= -1;

# ## Question 026
# Sum the array `1:4` with initial value of -10. (★☆☆)

sum(1:4, init = -10)

# ## Question 027
# Consider an integer vector `vZ` validate the following expressions. (★☆☆)
# ```julia
# vZ .^ vZ
# 2 << vZ >> 2
# vZ <- vZ
# 1im * vZ
# vZ / 1 / 1
# vZ < Z > Z
# ```

vZ = rand(1:10, 3);
#+
vZ .^ vZ
#+
try
    2 << vZ >> 2
catch e
    println(e)
end
#+
vZ <- vZ
#+
1im * vZ
#+
vZ / 1 / 1
#+
vZ < vZ > vZ

# ## Question 028
# Evaluate the following expressions. (★☆☆)

[0] ./ [0]
#+
try
    [0] .÷ [0]
catch e
    println(e)
end
#+
try
    convert(Float, convert(Int, NaN))
catch e
    println(e)
end

# ## Question 029
# Round away from zero a float array. (★☆☆)

vA = randn(10);
map(x -> x > 0 ? ceil(x) : floor(x), vA)

# ## Question 030
# Find common values between two arrays. (★☆☆)

vA = rand(1:10, 6);
vB = rand(1:10, 6);

vA[findall(in(vB), vA)]

# ## Question 031
# Suppress Julia's warnings. (★☆☆)

# One could use [Suppressor.jl](https://github.com/JuliaIO/Suppressor.jl).

# ## Question 032
# Compare `sqrt(-1)` and `sqrt(-1 + 0im)`. (★☆☆)

try
    sqrt(-1)
catch e
    println(e)
end
#+
sqrt(-1 + 0im)

# ## Question 033
# Display yesterday, today and tomorrow's date. (★☆☆)

println("Yesterday: $(today() - Day(1))");
println("Today: $(today())");
println("Tomorrow: $(today() + Day(1))");

# ## Question 034
# Display all the dates corresponding to the month of July 2016. (★★☆)

collect(Date(2016,7,1):Day(1):Date(2016,7,31))

# ## Question 035
# Compute `((mA + mB) * (-mA / 2))` in place. (★★☆)

mA = rand(2, 2);
mB = rand(2, 2);
mA .= ((mA .+ mB) .* (.-mA ./ 2))

# Using the dot macro:

@. mA = ((mA + mB) * (-mA / 2));

# ## Question 036
# Extract the integer part of a random array of positive numbers using 4 different methods. (★★☆)

mA = 5 * rand(3, 3);

# Option 1:
floor.(mA)

# Option 2:
round.(mA .- 0.5) #<! Generates -0.0 for numbers smaller than 0.5

# Option 3:
mA .÷ 1

# Option 4:
mA .- rem.(mA, 1)

# ## Question 037
# Create a `5x5` matrix with row values ranging from 0 to 4. (★★☆)

mA = repeat(reshape(0:4, 1, 5), 5, 1)

# One could also generate _row like_ range using tranpose:
mA = repeat((0:4)', 5, 1);

# ## Question 038
# Generate an array using a generator of 10 numbers. (★☆☆)

vA = collect(x for x in 1:10)

# In Julia the result of collect can be achieved directly using _Array Comprehension_:

vA = [x for x in 1:10];

# ## Question 039
# Create a vector of size 10 with values ranging from 0 to 1, both excluded. (★★☆)

vA = LinRange(0, 1, 12)[2:(end - 1)]

# ## Question 040
# Create a random vector of size 10 and sort it. (★★☆)

vA = rand(1:10, 10);
sort(vA) #<! Use `sort!()` for inplace sorting

# ## Question 041
# Implement the `sum()` function manually. (★★☆)

vA = rand(100);

function MySum(vA::Vector{T}) where{T <: Number}
sumVal = vA[1];
for ii in 2:length(vA)
    sumVal += vA[ii];
end

return sumVal;

end

MySum(vA)

# ## Question 042
# Check for equality of 2 arrays. (★★☆)

vA = rand(10);
vB = rand(10);

all(vA .== vB)

# ## Question 043
# Make an array immutable (Read only). (★★☆)

# This is a work in progress for Julia at in [Issue 31630](https://github.com/JuliaLang/julia/pull/31630).

# ## Question 044
# Consider a random `10x2` matrix representing cartesian coordinates, convert them to polar coordinates. (★★☆)

mA = rand(10, 2);

ConvToPolar = vX -> [hypot(vX[1], vX[2]), atan(vX[2], vX[1])]

mB = [ConvToPolar(vX) for vX in eachrow(mA)]

# In order to have the same output size:

mC = reduce(hcat, mB)';

# ## Question 045
# Create random vector of size 10 and replace the maximum value by 0. (★★☆)

vA = randn(10);

# In case of a single maximum or all different values:
vA[argmax(vA)] = 0;
vA

# General solution:

maxVal = maximum(vA);
vA .= (valA == maxVal ? 0 : valA for valA in vA); #<! Non allocating generator by using `.=`

# ## Question 046
# Create a a grid of `x` and `y` coordinates covering the `[0, 1] x [0, 1]` area. (★★☆)

numGridPts = 5;
vX = LinRange(0, 1, numGridPts);
vY = LinRange(0, 1, numGridPts);
MeshGrid = (vX, vY) -> ([x for _ in vY, x in vX], [y for y in vY, _ in vX]);

mX, mY = MeshGrid(vX, vY); #<! See https://discourse.julialang.org/t/48679
@show mX
#+
@show mY

# By Tomer Arnon (https://github.com/tomerarnon):

mXY = [(ii, jj) for ii in 0:0.25:1, jj in 0:0.25:1]; #<! Also `tuple.(0:0.25:1, (0:0.25:1)')`

# ## Question 047
# Given two vectors, `vX` and `vY`, construct the Cauchy matrix `mC`: `(Cij = 1 / (xi - yj))`. (★★☆)

vX = rand(5);
vY = rand(5);

mC = 1 ./ (vX .- vY')

# ## Question 048
# Print the minimum and maximum representable value for each Julia scalar type. (★★☆)

vT = [UInt8 UInt16 UInt32 UInt64 Int8 Int16 Int32 Int64 Float16 Float32 Float64]

for juliaType in vT
    println(typemin(juliaType));
    println(typemax(juliaType));
end

# ## Question 049
# Print all the values of an array. (★★☆)

mA = rand(3, 3);
print(mA);

# ## Question 050
# Find the closest value to a given scalar in a vector. (★★☆)

inputVal = 0.5;
vA = rand(10);

vA[argmin(abs.(vA .- inputVal))]

# By Tomer Arnon (https://github.com/tomerarnon):

function ClosestValue(vA::Vector{T}, inputVal::T) where{T <: Number}
    return argmin(y -> abs(y - inputVal), vA);
end

ClosestValue(vA, inputVal)

# ## Question 051
# Create a structured array representing a position `(x, y)` and a color `(r, g, b)`. (★★☆)

struct sPosColor
    x::Int
    y::Int
    R::UInt8;
    G::UInt8;
    B::UInt8;
    A::UInt8;
end

numPixels   = 10;
maxVal      = typemax(UInt32);
vMyColor    = [sPosColor(rand(1:maxVal, 2)..., rand(UInt8, 4)...) for _ in 1:numPixels];

# ## Question 052
# Consider a random vector with shape `(5, 2)` representing coordinates, find the distances matrix `mD`: $ {D}_{i, j} = {\left\| {x}_{i} - {x}_{j} \right\|}_{2} $. (★★☆)

mX = rand(5, 2);
vSumSqr = sum(vX -> vX .^ 2, mX, dims = 2);
mD = vSumSqr .+ vSumSqr' - 2 * (mX * mX');
mD #<! Apply `sqrt.()` for the actual norm

# ## Question 053
# Convert a float (32 bits) array into an integer (32 bits) in place. (★★☆)

vA = 9999 .* rand(Float32, 5);
vB = reinterpret(Int32, vA); #<! Creates a view
@. vB = trunc(Int32, vA) #<! Updates the byes in th view (Inplace for `vA`)
# The above is equivalent of:
# ```julia
# for ii in eachindex(vB)
#     vB[ii] = trunc(Int32, vA[ii]);
# end
# vB
# ```

# ## Question 054
# Read the following file (`Q0054.txt`). (★★☆)
# ```
# 1, 2, 3, 4, 5
# 6,  ,  , 7, 8
#  ,  , 9,10,11
# ```

mA = readdlm("Q0054.txt", ',')

# ## Question 055
# Enumerate array in a loop. (★★☆)

mA = rand(3, 3);

for (elmIdx, elmVal) in enumerate(mA) #<! See https://discourse.julialang.org/t/48877
    println(elmIdx);
    println(elmVal);
end

# ## Question 056
# Generate a generic 2D Gaussian like array with `μ = 0`, `σ = 1` and indices over `{-5, -4, ..., 0, 1, ..., 5}`. (★★☆)

vA = -5:5;
μ = 0;
σ = 1;
mG = [(1 / (2 * pi * σ)) * exp(-0.5 * ((([x, y] .- μ)' * ([x, y] .- μ)) / (σ * σ))) for x in vA, y in vA];

heatmap(mG)

# Using the separability of the Gaussian function:

vG = (1 / (sqrt(2 * pi) * σ)) .* exp.(-0.5 .* (((vA .- μ) .^ 2) / (σ * σ)));
mG = vG * vG';

# ## Question 057
# Place `5` elements in a `5x5` array randomly. (★★☆)

mA = rand(5, 5);
mA[rand(1:25, 5)] = rand(5);

# Another option which avoids setting into the same indices:
mA[randperm(25)[1:5]] = rand(5);

# ## Question 058
# Subtract the mean of each row of a matrix. (★★☆)

mA = rand(3, 3);
mA .-= mean(mA, dims = 2);
mean(mA, dims = 1)

# ## Question 059
# Sort an array by a column. (★★☆)

colIdx = 2;

mA = rand(3, 3);
mA[sortperm(mA[:, colIdx]), :]

# Using `sortslices()`:
sortslices(mA, dims = 1, by = x -> x[colIdx]);

# ## Question 060
# Tell if a given 2D array has null (All zeros) columns. (★★☆)

mA = rand(0:1, 3, 9);
any(all(iszero.(mA), dims = 1))

# ## Question 061
# Find the 2nd nearest value from a given value in an array. (★★☆)

inputVal = 0.5;
vA = rand(10);

vA[sortperm(abs.(vA .- inputVal))[2]]

# Alternative way (More efficient)

closeFirst  = Inf;
closeSecond = Inf;
closeFirstIdx  = 0;
closeSecondIdx = 0;

## Using `global` for scope in Literate
for (elmIdx, elmVal) in enumerate(abs.(vA .- inputVal))
    if (elmVal < closeFirst)
        global closeSecond = closeFirst;
        global closeFirst = elmVal;
        global closeSecondIdx  = closeFirstIdx;
        global closeFirstIdx   = elmIdx;
    elseif (elmVal < closeSecond)
        global closeSecond = elmVal;
        global closeSecondIdx = elmIdx;
    end
end

vA[closeSecondIdx] == vA[sortperm(abs.(vA .- inputVal))[2]]

# By Tomer Arnon (https://github.com/tomerarnon):

vA[partialsortperm(abs.(vA .- inputVal), 2)]

# ## Question 062
# Considering two arrays with shape `(1, 3)` and `(3, 1)`, Compute their sum using an iterator. (★★☆)

vA = rand(1, 3);
vB = rand(3, 1);

sum(aVal + bVal for aVal in vA, bVal in vB)

# ## Question 063
# Create an array class that has a name attribute. (★★☆)
#+
# One could use `NamedArrays.jl` or `AxisArrays.jl`.

# ## Question 064
# Given a vector, add `1` to each element indexed by a second vector (Be careful with repeated indices). (★★★)

vA = rand(1:10, 5);
vB = rand(1:5, 3);

println(vA);

## Julia is very efficient with loops
for bIdx in vB
    vA[bIdx] += 1;
end

println(vA);

# ## Question 065
# Accumulate elements of a vector `X` to an array `F` based on an index list `I`. (★★★)

vX = rand(1:5, 10);
vI = rand(1:15, 10);

numElements = maximum(vI);
vF = zeros(numElements);

for (ii, iIdx) in enumerate(vI)
    vF[iIdx] += vX[ii];
end

println("vX: $vX");
println("vI: $vI");
println("vF: $vF");

# One could also use `counts()` from `StatsBase.jl`.

# ## Question 066
# Considering an image of size `w x h x 3` image of type `UInt8`, compute the number of unique colors. (★★☆)

mI = rand(UInt8, 1000, 1000, 3);

numColors = length(unique([reinterpret(UInt32, [iPx[1], iPx[2], iPx[3], 0x00])[1] for iPx in eachrow(reshape(mI, :, 3))])); #<! Reshaping as at the moment `eachslice()` doesn't support multiple `dims`.
print("Number of Unique Colors: $numColors");

# Another option:

numColors = length(unique([UInt32(iPx[1]) + UInt32(iPx[2]) << 8 + UInt32(iPx[3]) << 16 for iPx in eachrow(reshape(mI, :, 3))]));
print("Number of Unique Colors: $numColors");

# Simpler way to slice a pixel:

numColors = length(unique([UInt32(mI[ii, jj, 1]) + UInt32(mI[ii, jj, 2]) << 8 + UInt32(mI[ii, jj, 3]) << 16 for ii in 1:size(mI, 1), jj in 1:size(mI, 2)]));
print("Number of Unique Colors: $numColors");

# ## Question 067
# Considering a four dimensions array, get sum over the last two axis at once. (★★★)

mA = rand(2, 2, 2, 2);
sum(reshape(mA, (2, 2, :)), dims = 3)

# ## Question 068
# Considering a one dimensional vector `vA`, how to compute means of subsets of `vA` using a vector `vS` of same size describing subset indices. (★★★)

## Bascially extending `Q0065` with another vector of number of additions.

vX = rand(1:5, 10);
vI = rand(1:15, 10);

numElements = maximum(vI);
vF = zeros(numElements);
vN = zeros(Int, numElements);

for (ii, iIdx) in enumerate(vI)
    vF[iIdx] += vX[ii];
    vN[iIdx] += 1;
end

## We only divide the mean if the number of elements accumulated is bigger than 1
for ii in 1:numElements
    vF[ii] = ifelse(vN[ii] > 1, vF[ii] / vN[ii], vF[ii]);
end

println("vX: $vX");
println("vI: $vI");
println("vF: $vF");

# ## Question 069
# Get the diagonal of a matrix product. (★★★)

mA = rand(5, 7);
mB = rand(7, 4);

numDiagElements = min(size(mA, 1), size(mB, 2));
vD = [dot(mA[ii, :], mB[:, ii]) for ii in 1:numDiagElements]

# Alternative way:

vD = reshape(sum(mA[1:numDiagElements, :]' .* mB[:, 1:numDiagElements], dims = 1), numDiagElements)

# ## Question 070
# Consider the vector `[1, 2, 3, 4, 5]`, build a new vector with 3 consecutive zeros interleaved between each value. (★★★)

vA = 1:5;

## Since Julia is fast with loops, it would be the easiest choice

numElements = (4 * length(vA)) - 3;
vB = zeros(Int, numElements);

for (ii, bIdx) in enumerate(1:4:numElements)
    vB[bIdx] = vA[ii];
end
println(vB);

## Alternative (MATLAB style) way:

mB = [reshape(collect(vA), 1, :); zeros(Int, 3, length(vA))];
vB = reshape(mB[1:(end - 3)], :);
println(vB);

# ## Question 071
# Consider an array of dimension `5 x 5 x 3`, mulitply it by an array with dimensions `5 x 5` using broadcasting. (★★★)

mA = rand(5, 5, 3);
mB = rand(5, 5);

mA .* mB #<! Very easy in Julia

# ## Question 072
# Swap two rows of a 2D array. (★★★)

mA = rand(UInt8, 3, 2);
println(mA);
mA[[1, 2], :] .= mA[[2, 1], :];
println(mA);

# ## Question 073
# Consider a set of 10 triplets describing 10 triangles (with shared vertices), find the set of unique line segments composing all the triangles. (★★★)

mA = rand(0:100, 10, 3); #<! Each row composes 3 veritces ([1] -> [2], [2] -> [3], [3] -> [1])
mC = [sort!([vC[mod1(ii, end)], vC[mod1(ii + 1, end)]]) for ii in 1:(size(mA, 2) + 1), vC in eachrow(mA)][:] #<! Sorted combinations of vertices
mC = unique(mC)

# ## Question 074
# Given a sorted array `vC` that corresponds to a bincount, produce an array `vA` such that `bincount(vA) == vC`. (★★★)

vC = rand(0:7, 5);
numElements = sum(vC);
vA = zeros(Int, numElements);

elmIdx = 1;
## Using `global` for scope in Literate
for (ii, binCount) in enumerate(vC)
    for jj in 1:binCount
        vA[elmIdx] = ii;
        global elmIdx += 1;
    end
end

# ## Question 075
# Compute averages using a sliding window over an array. (★★★)

numElements = 10;
winRadius   = 1;
winReach    = 2 * winRadius;
winLength   = 1 + winReach;

vA = rand(0:3, numElements);
vB = zeros(numElements - (2 * winRadius));

aIdx = 1 + winRadius;
## Using `global` for scope in Literate
for ii in 1:length(vB)
    vB[ii] = mean(vA[(aIdx - winRadius):(aIdx + winRadius)]); #<! Using integral / running sum it would be faster.
    global aIdx += 1;
end

# Another method using running sum:

vC = zeros(numElements - winReach);

jj = 1;
sumVal = sum(vA[1:winLength]);
vC[jj] = sumVal / winLength;
jj += 1;

## Using `global` for scope in Literate
for ii in 2:(numElements - winReach)
    global sumVal += vA[ii + winReach] - vA[ii - 1];
    vC[jj] = sumVal / winLength;
    global jj += 1;
end

maximum(abs.(vC - vB)) < 1e-8

# ## Question 076
# Consider a one dimensional array `vA`, build a two dimensional array whose first row is `[ vA[0], vA[1], vA[2] ]`  and each subsequent row is shifted by 1. (★★★)

vA = rand(10);
numCols = 3;

numRows = length(vA) - numCols + 1;
mA = zeros(numRows, numCols);

for ii in 1:numRows
    mA[ii, :] = vA[ii:(ii + numCols - 1)]; #<! One could optimize the `-1` out
end

# ## Question 077
# Negate a boolean or to change the sign of a float inplace. (★★★)

vA = rand(Bool, 10);
vA .= .!vA;

vA = randn(10);
vA .*= -1;

# ## Question 078
# Consider 2 sets of points `mP1`, `mP2` describing lines (2d) and a point `vP`, how to compute distance from the point `vP` to each line `i`: `[mP1[i, :], mP2[i, :]`. (★★★)

## See distance of a point from a line in Wikipedia (https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line).
## Specifically _Line Defined by Two Points_.

numLines = 10;
mP1 = randn(numLines, 2);
mP2 = randn(numLines, 2);
vP  = randn(2);

vD = [(abs(((vP2[1] - vP1[1]) * (vP1[2] - vP[2])) - ((vP1[1] - vP[1]) * (vP2[2] - vP1[2]))) / hypot((vP2 - vP1)...)) for (vP1, vP2) in zip(eachrow(mP1), eachrow(mP2))];
minDist = minimum(vD);

println("Min Distance: $minDist");

# ## Question 079
# Consider 2 sets of points `mP1`, `mP2` describing lines (2d) and a set of points `mP`, how to compute distance from the point `vP = mP[j, :]` to each line `i`: `[mP1[i, :], mP2[i, :]`. (★★★)

numLines = 5;
mP1 = randn(numLines, 2);
mP2 = randn(numLines, 2);
mP  = randn(numLines, 2);

mD = [(abs(((vP2[1] - vP1[1]) * (vP1[2] - vP[2])) - ((vP1[1] - vP[1]) * (vP2[2] - vP1[2]))) / hypot((vP2 - vP1)...)) for (vP1, vP2) in zip(eachrow(mP1), eachrow(mP2)), vP in eachrow(mP)];

for jj in 1:numLines
    minDist = minimum(mD[jj, :]);
    println("The minimum distance from the $jj -th point: $minDist");
end

# ## Question 080
# Consider an arbitrary 2D array, write a function that extract a subpart with a fixed shape and centered on a given element (Handel out of bounds). (★★★)

## One could use `PaddedViews.jl` to easily solve this.

arrayLength = 10;
winRadius   = 3;
vWinCenter  = [7, 9];

mA = rand(arrayLength, arrayLength);
winLength = (2 * winRadius) + 1;
mB = zeros(winLength, winLength);

verShift = -winRadius;
## Using `global` for scope in Literate
for ii in 1:winLength
    horShift = -winRadius;
    for jj in 1:winLength
        mB[ii, jj] = mA[min(max(vWinCenter[1] + verShift, 1), arrayLength), min(max(vWinCenter[2] + horShift, 1), arrayLength)]; #<! Nearest neighbor extrapolation
        horShift += 1;
    end
    global verShift += 1;
end

# ## Question 081
# Consider an array `vA = [1, 2, 3, ..., 13, 14]`, generate an array `vB = [[1, 2, 3, 4], [2, 3, 4, 5], [3, 4, 5, 6], ..., [11, 12, 13, 14]]`. (★★★)

vA = collect(1:14);

winNumElements  = 4;
winReach        = winNumElements - 1;

vB = [vA[ii:(ii + winReach)] for ii in 1:(length(vA) - winReach)]

# ## Question 082
# Compute a matrix rank. (★★★)

numRows = 5;
numCols = 4;
mA = randn(numRows, numCols);
rank(mA)

# ## Question 083
# Find the most frequent value in an array. (★★★)

vA = rand(1:5, 15);

# MATLAB Style (Manual loop might be faster)

vB = unique(vA);
## vB[argmax(sum(vA .== vB', dims = 1)[:])] #<! The input to `argmax()` is a `1 x n` vector, hence squeezed so `argmax()` won't return Cartesian Index.
vB[argmax(dropdims(sum(vA .== vB', dims = 1), dims = 1))] #<! The input to `argmax()` is a `1 x n` vector, hence squeezed so `argmax()` won't return Cartesian Index.

# Comparing bits:

# One could convert at the bits level to integers and then use something like `counts()` from `StatsBase.jl`.
# Support to 1:4 bytes of data:
# ```julia
# numBytes = sizeof(vA[1]);
# if (sizeof(vA[1]) == 1)
#     vB = reinterpret(UInt8, vA);
# elseif (sizeof(vA[1]) == 2)
#     vB = reinterpret(UInt16, vA);
# elseif (sizeof(vA[1]) == 4)
#     vB = reinterpret(UInt32, vA);
# elseif (sizeof(vA[1]) == 8)
#     vB = reinterpret(UInt64, vA);
# end
# ```

# ## Question 084
# Extract all the contiguous `3x3` blocks from a random `5x5` matrix. (★★★)

numRows = 5;
numCols = 5;

mA = rand(1:9, numRows, numCols);

winRadius   = 1;
winReach    = 2 * winRadius;
winLength   = winReach + 1;

mB = [mA[ii:(ii + winReach), jj:(jj + winReach)] for ii in 1:(numRows - winReach), jj in 1:(numCols - winReach)]

# ## Question 085
# Create a 2D array struct such that `mA[i, j] == mA[j, i]` (Symmetric matrix). (★★★)

struct SymmetricMatrix{T <: Number} <: AbstractArray{T, 2}
    numRows::Int
    data::Matrix{T}

    function SymmetricMatrix(mA::Matrix{T}) where {T <: Number}
        size(mA, 1) == size(mA, 2) || throw(ArgumentError("Input matrix must be square"))
        new{T}(size(mA, 1), Matrix(Symmetric(mA)));
    end

end

function Base.size(mA::SymmetricMatrix)
    (mA.numRows, mA.numRows);
end
function Base.getindex(mA::SymmetricMatrix, ii::Int)
    mA.data[ii];
end
function Base.getindex(mA::SymmetricMatrix, ii::Int, jj::Int)
    mA.data[ii, jj];
end
function Base.setindex!(mA::SymmetricMatrix, v, ii::Int, jj::Int)
    setindex!(mA.data, v, ii, jj);
    setindex!(mA.data, v, jj, ii);
end

mA = SymmetricMatrix(zeros(Int, 2, 2));
mA[1, 2] = 5;
mA

# ## Question 086
# Consider a set of `p` matrices of shape `nxn` and a set of `p` vectors with length `n`. Compute the sum of of the `p` matrix vector products at once (Result is a vector of length `n`). (★★★)

## One could use `TensorOperations.jl` or `Einsum.jl` for a more elegant solution.

numRows = 5;
numMat  = 3;

tP = [randn(numRows, numRows) for _ in 1:numMat];
mP = [randn(numRows) for _ in 1:numMat];

vA = reduce(+, (mP * vP for (mP, vP) in zip(tP, mP)));

#+ Vanilla solution

vB = zeros(numRows);
for ii in 1:numMat
    vB .+= tP[ii] * mP[ii];
end

vA == vB

# ## Question 087
# Consider a `16x16` array, calculate the block sum (Block size is `4x4`). (★★★)

# We solve a more general case for any size of blocks.

numRows = 16;
numCols = 8;

vBlockSize = [2, 4]; #<! [numRows, numCols] ./ vBlockSize == integer

mA = rand(numRows, numCols);

numBlocksVert   = numRows ÷ vBlockSize[1];
numBlocksHori   = numCols ÷ vBlockSize[2];
numBlocks       = numBlocksVert * numBlocksHori;

mA = reshape(mA, vBlockSize[1], :);

# We number the blocks column wise and create their block index per column of the reshaped `mA`.
vBlockIdx = 1:numBlocks;
mBlockIdx = reshape(vBlockIdx, numBlocksVert, numBlocksHori);
mBlockIdx = repeat(mBlockIdx, 1, 1, vBlockSize[2]);
mBlockIdx = permutedims(mBlockIdx, (1, 3, 2));
vBlockIdx = mBlockIdx[:]; #<! Matches the block index per column of the reshaped `mA`.

vA = dropdims(sum(mA, dims = 1), dims = 1);
vB = zeros(numBlocks);

for ii = 1:(numBlocks * vBlockSize[2])
    vB[vBlockIdx[ii]] += vA[ii];
end
vB

# ## Question 088
# Implement the simulation _Game of Life_ using arrays. (★★★)

numRows = 20;
numCols = 20;

gofKernel           = @kernel w -> w[-1, -1] + w[-1, 0] + w[-1, 1] + w[0, -1] + w[0, 1] + w[1, -1] + w[1, 0] + w[1, 1];
gofNumLives         = round(Int, 0.05 * numRows * numCols);
gofNumGenerations   = 50;

vI = randperm(numRows * numCols)[1:gofNumLives];

mG = zeros(UInt8, numRows, numCols);
mG[vI] .= UInt8(1);
mB = similar(mG);

heatmap(mG) #<! Initialization

#+

for ii in 1:numGridPts
    map!(gofKernel, mB, extend(mG, StaticKernels.ExtensionConstant(0))); #<! One may use `ExtensionCircular`
    for ii in eachindex(mB)
        mG[ii] = UInt8((mB[ii] >= 3) || ((mB[ii] > 0) && (mB[ii] == 2))); #<! Could be done with broadcasting
    end
end

heatmap(mG) #<! Final state

# TODO: Use O(1) implementation for Box Blur.

# ## Question 089
# Get the `n` largest values of an array. (★★★)

vA = rand(10);
numValues = 3;

vA[partialsortperm(vA, 1:numValues, rev = true)]

# ## Question 090
# Given an arbitrary number of vectors, build the _Cartesian Product_ (Every combinations of every item). (★★★)

function CartesianProduct(tupleX)
    return collect(Iterators.product(tupleX...))[:];
end

vA = 1:3;
vB = 8:9;
vC = 4:5;

CartesianProduct((vA, vB, vC))

# ## Question 091
# Create an array which can be accessed like a _record array_ in _NumPy_. (★★★)

# One could use `StructArrays.jl`.

# ## Question 092
# Consider a large vector `vA`, compute `vA` to the power of 3 using 3 different methods. (★★★)

vA = rand(1000);

# Method 001:

vB = vA .^ 3;

# Method 002:

vC = [valA ^ 3 for valA in vA];

# Method 003:

vD = zeros(length(vA));
for (ii, valA) in enumerate(vA)
    vD[ii] = valA * valA * valA;
end

#+

vB ≈ vC ≈ vD

# ## Question 093
# Consider two arrays `mA` and `mB` of shape `8x3` and `2x2`. Find rows of `mA` that contain elements of each row of `mB` regardless of the order of the elements in `mB`. (★★★)

# The way I interpret the question is rows in `mA` which contain at least 1 element from each row of `mB`.

mA = rand(0:4, 8, 3);
mB = rand(0:4, 2, 2);
mC = [any(vA .== vB') for vB in eachrow(mB), vA in eachrow(mA)]; #<! General solution, will work for any size of `mA` and `mB`
vD = [all(vC) for vC in eachcol(mC)]

# In order to have a solution without the intermediate array `mC`

function Iterate2(iA; iterState = missing)
    if(ismissing(iterState))
        valA, iterState = iterate(iA);
    else
        valA, iterState = iterate(iA, iterState);
    end
    valB, iterState = iterate(iA, iterState);
    return (valA, valB), iterState
end

tT = (any(vA .== vB') for vB in eachrow(mB), vA in eachrow(mA));

iterState = missing;

vE = zeros(Bool, size(mA, 1));

for ii = 1:length(vD)
    global iterState;
    (valA, valB), iterState = Iterate2(tT; iterState = iterState);
    vE[ii] = valA && valB;
end

vD == vE

# ## Question 094
# Considering a `10x3` matrix, extract rows with unequal values. (★★★)

mA = rand(1:3, 10, 3);
vD = [maximum(vA) != minimum(vA) for vA in eachrow(mA)]

# ## Question 095
# Convert a vector of ints into a matrix binary representation. (★★★)

vA = rand(UInt8, 10);
mB = zeros(Bool, length(vA), 8);

## See https://discourse.julialang.org/t/26663
for ii in 1:length(vA)
    vS = bitstring(vA[ii]);
    for jj in 1:size(mB, 2)
        mB[ii, jj] = vS[jj] == '1';
    end
end

mB

# By Tomer Arnon (https://github.com/tomerarnon):

mB = reverse!(reduce(hcat, digits.(vA, base = 2, pad = 8))', dims = 2);

# ## Question 096
# Given a two dimensional array, extract unique rows. (★★★)

mA = UInt8.(rand(1:3, 10, 3));

vS = [reduce(*, bitstring(valA) for valA in vA) for vA in eachrow(mA)]; #<! Supports any array!
vU = unique(vS);
vI = [findfirst(valU .== vS) for valU in vU];

# An alternative way:
vB = indexin(vU, vS);
vB == vI

# ## Question 097
# Considering 2 vectors `vA` and `vB`, write the einsum equivalent (Using `Einsum.jl`) of inner, outer, sum, and mul function. (★★★)

vA = rand(5);
vB = rand(5);

# Inner Product
@tullio tullioVal = vA[ii] * vB[ii];
tullioVal ≈ dot(vA, vB) #<! Inner product

# Outer  Product
@tullio mTullio[ii, jj] := vA[ii] * vB[jj]; #<! Outer product
mTullio ≈ vA * vB'

# Sum
@tullio tullioVal = vA[ii];
tullioVal ≈ sum(vA) #<! Sum

# Multiplication
@tullio vTullio[ii] := vA[ii] * vB[ii];
vTullio ≈ vA .* vB #<! Multiplication

# ## Question 098
# Considering a path described by two vectors `vX` and `vY`, sample it using equidistant samples. (★★★)

# The way I interpreted the question is to create sub segments of the same length.

numPts      = 100;
numSegments = 1000;

vX = sort(10 * rand(numPts));
vY = sort(10 * rand(numPts));

vR = cumsum(hypot.(diff(vX), diff(vY)));
pushfirst!(vR, 0.0);
vRSegment = LinRange(0.0, vR[end], numSegments);

struct LinearInterpolator1D{T <: Number} <: AbstractArray{T, 1}
    vX::Vector{T};
    vY::Vector{T};

    function LinearInterpolator1D(vX::Vector{T}, vY::Vector{T}) where {T <: Number}
        length(vX) == length(vX) || throw(ArgumentError("Input vectors must have the same length"));
        new{T}(vX, vY);
    end
end

function Base.size(vA::LinearInterpolator1D)
    size(vA.vX);
end
function Base.getindex(vA::LinearInterpolator1D, ii::Number)
    if (ii >= vA.vX[end])
        return vA.vY[end];
    end
    if (ii <= vA.vX[1])
        return vA.vY[1];
    end

    rightIdx = findfirst(vA.vX .>= ii);
    leftIdx = rightIdx - 1;

    tt = (ii - vA.vX[leftIdx]) / (vA.vX[rightIdx] - vA.vX[leftIdx]);

    return ((1 - tt) * vA.vY[leftIdx]) + (tt * vA.vY[rightIdx]);

end
function Base.setindex!(vA::LinearInterpolator1D, valX, valY, ii::Int, jj::Int)
    setindex!(sLinInterp.vX, valX, ii);
    setindex!(sLinInterp.vY, valY, ii);
end

vXInt = LinearInterpolator1D(vR, vX);
vYInt = LinearInterpolator1D(vR, vY);

vXSegment = [vXInt[intIdx] for intIdx in vRSegment];
vYSegment = [vYInt[intIdx] for intIdx in vRSegment];

hP = lineplot(vX, vY, canvas = DotCanvas, name = "Samples");
lineplot!(hP, vXSegment, vYSegment, name = "Interpolated");
hP

# **Remark**: In order to be mathematically accurate the resolution of `vRSegment` must be high enough to be an integer factor of each segment.
# It can be done if the resolution is `1 / (prod(vR))` which is easily infeasible with `Float64`. So this is a good enough approximation.

# ## Question 099
# Given an integer `n` and a 2D array `mA`, find the rows which can be interpreted as draws from a multinomial distribution with `n` (Rows which only contain integers and which sum to `n`). (★★★)

mA = rand([0, 0.5, 1, 2, 3], 15, 3);
sumVal = 4;
vI = [all(vA .== round.(vA)) && sum(vA) == sumVal for vA in eachrow(mA)];

# ## Question 100
# Compute bootstrapped `95%` confidence intervals for the mean of a 1D array `vA`. Namely, resample the elements of an array with replacement `N` times, compute the mean of each sample and then compute percentiles over the means. (★★★)

numTrials   = 1000;
numSamples  = 100;
μ           = 0.5;

vA = μ .+ randn(numSamples);
tM = (mean(vA[rand(1:numSamples, numSamples)]) for _ in 1:numTrials);
quantile(tM, [0.025, 0.975])

ext: jl


text:
using Documenter, DocumenterTools
include("DocumenterShowcase.jl")

# The DOCSARGS environment variable can be used to pass additional arguments to make.jl.
# This is useful on CI, if you need to change the behavior of the build slightly but you
# can not change the .travis.yml or make.jl scripts any more (e.g. for a tag build).
if haskey(ENV, "DOCSARGS")
    for arg in split(ENV["DOCSARGS"])
        (arg in ARGS) || push!(ARGS, arg)
    end
end

makedocs(
    modules = [Documenter, DocumenterTools, DocumenterShowcase],
    format = Documenter.HTML(
        # Use clean URLs, unless built as a "local" build
        prettyurls = !("local" in ARGS),
        canonical = "https://juliadocs.github.io/Documenter.jl/stable/",
        assets = ["assets/favicon.ico"],
        analytics = "UA-136089579-2",
        highlights = ["yaml"],
    ),
    clean = false,
    sitename = "Documenter.jl",
    authors = "Michael Hatherly, Morten Piibeleht, and contributors.",
    linkcheck = !("skiplinks" in ARGS),
    pages = [
        "Home" => "index.md",
        "Manual" => Any[
            "Guide" => "man/guide.md",
            "man/examples.md",
            "man/syntax.md",
            "man/doctests.md",
            "man/latex.md",
            hide("man/hosting.md", [
                "man/hosting/walkthrough.md"
            ]),
            "man/other-formats.md",
        ],
        "showcase.md",
        "Library" => Any[
            "Public" => "lib/public.md",
            "Internals" => map(
                s -> "lib/internals/$(s)",
                sort(readdir(joinpath(@__DIR__, "src/lib/internals")))
            ),
        ],
        "contributing.md",
    ],
    strict = !("strict=false" in ARGS),
    doctest = ("doctest=only" in ARGS) ? :only : true,
)

deploydocs(
    repo = "github.com/JuliaDocs/Documenter.jl.git",
    target = "build",
)


text:
module ITensorsGPU

using CuArrays
using CuArrays.CUTENSOR
using CuArrays.CUBLAS
using CuArrays.CUSOLVER
using LinearAlgebra
using Random
using TimerOutputs
using StaticArrays
using ITensors
import CuArrays: CuArray, CuMatrix, CuVector
import CuArrays.CUTENSOR: cutensorContractionPlan_t, cutensorAlgo_t
import ITensors: randn!, compute_contraction_labels,
                 plussers, eigen, similar_type, tensor,
                 scale!, unioninds, array, matrix, vector,
                 polar, tensors, truncate!, leftlim, rightlim
import ITensors.NDTensors: ContractionProperties, contract!!, _contract!!, _contract!, contract!, contract,
                         contraction_output, UniformDiagTensor, CombinerTensor, contraction_output_type,
                         UniformDiag, Diag, DiagTensor, NonuniformDiag, NonuniformDiagTensor, zero_contraction_output,
                         outer!, outer!!, is_trivial_permutation, ind, permutedims!!, Dense, DenseTensor, Combiner,
                         Tensor, data, permute, getperm
import Base.*, Base.permutedims!
include("tensor/cudense.jl")
include("tensor/culinearalgebra.jl")
include("tensor/cutruncate.jl")
include("tensor/cucombiner.jl")
include("tensor/cudiag.jl")
include("cuitensor.jl")
include("mps/cumps.jl")
include("mps/cumpo.jl")

#const ContractionPlans = Dict{String, Tuple{cutensorAlgo_t, cutensorContractionPlan_t}}()
const ContractionPlans = Dict{String, cutensorAlgo_t}()

export cuITensor,
       randomCuITensor,
       cuMPS,
       randomCuMPS,
       productCuMPS,
       randomCuMPO,
       cuMPO
end #module


text:
using TestFunctionRunner
TestFunctionRunner.@run(paths = ["../benchmark/ConcurrentUtilsBenchmarks"])


text:
"""
    Rogue.usein(downpath; dryrun, from, rev, commit, push)

Update `(Julia)Manifest.toml` file(s) in a downstream project at
`downpath` to use the current version of the upstream project.

* Make sure that the downstream project has no un-committed changes.

* If `manifests` is `nothing` (default), find all `JuliaManifest.toml` and
  `Manifest.toml` file(s) that are tracked by git and have the
  upstream project as a dependency.

* Update the downstream manifest files (update `git-tree-sha1`).

* Resolve dependencies.

    * Note that this makes sure that the current version of the
      upstream project is available in the repository referenced by
      the downstream manifest files (i.e., `git push`ed).

* Commit the changed manifest files with a git commit message
  generated from the upstream commit.  In particular, it contains a
  URL to the commit page of the VCS hosting service used by the
  upstream project.

# Arguments
- `downpath :: AbstractString`: Path to the downstream project.

# Keyword Arguments
- `dryrun :: Bool = false`: If `true`, only print the operations that
  would be performed.

- `from :: AbstractString = "."`: Specify the location of the upstream
  project.

- `rev :: AbstractString = "HEAD"`: Revision of the upstream project.

- `commit :: Union{Bool, Cmd} = true`: If it is a `Bool`, it
  determines if the change should be committed.  If it is a `Cmd`, the
  change is committed and this is passed as options to the `git
  commit` command.

- `push :: Union{Bool, Cmd} = false`: Similar to `commit` but for `git
  push`.
"""
function usein(
    downpath::AbstractString;
    dryrun::Bool = false,
    from::AbstractString = ".",
    rev::AbstractString = "HEAD",
    commit::Union{Bool, Cmd} = true,
    push::Union{Bool, Cmd} = false,
    committitle::Union{Nothing, AbstractString} = nothing,
    commitcomment::Union{Nothing, AbstractString} = nothing,
)

    uppkgid = pkgat(from)
    fullrev = strip(read(git_cmd(`rev-parse $rev`, from), String))
    treesha1 = strip(read(git_cmd(`rev-parse "$fullrev^{tree}"`, from), String))
    manifests = find_downstream_manifests(downpath, uppkgid)
    if isempty(manifests)
        @error "No git-tracked manifest files containing $uppkgid found in $downpath"
        return
    end
    if !git_is_clean(downpath)
        error("Git repository at `$downpath` has un-committed files.")
    end

    # Preparing commit message now so that it works nicely with
    # `dryrun=true`:
    downroot = strip(read(git_cmd(`rev-parse --show-toplevel`, downpath), String))
    commitfiles = relpath.(manifests, downroot)
    msg = commitmessage(fullrev, uppkgid, from, committitle, commitcomment)
    commitargs = `commit --message $msg`
    if !(committitle === nothing && commitcomment === nothing)
        commitargs = `$commitargs --allow-empty`
    end
    if commit isa Cmd
        commitargs = `$commitargs $commit`
    end
    pushargs = `push`
    if push isa Cmd
        pushargs = `$pushargs $push`
    end

    if dryrun
        for path in manifests
            @info "(dry-run) $path would be updated."
        end
        # TODO: somehow do this inside `update_manifest`.

        commit == false && return
        # TODO: parse `commitargs` to extract the correct msg
        @info """
        (dry-run) Commit manifest files with message:

        $msg
        """

        push === false && return
        @info """
        (dry-run) Pushing changes to remote...
        Execute: $(`git $pushargs`)
        """
        return
    end

    for path in manifests
        @info "Updating $path"
        update_manifest(path, uppkgid, treesha1)
    end
    if git_is_clean(downpath) && committitle === nothing && commitcomment === nothing
        @info "No updates were required in `$downpath`."
        return
    end
    if commit === false
        @info "`commit=false` is specified. Skipping commit."
        return
    end

    run(git_cmd(`add -- $commitfiles`, downroot))
    run(git_cmd(`$commitargs -- $commitfiles`, downroot))

    @info "Successfully update manifest files:\n$(join(manifests, "\n"))"

    push === false && return
    @info "Pushing changes to remote..."
    run(git_cmd(pushargs, downroot))
    @info "Pushing changes to remote... DONE"
end

"""
    Rogue.updateto(upstream; kwargs...)

A shortcut for `Rogue.usein(".", from=upstream)`.
"""
function updateto(upstream; downpath=".", from=nothing, kwargs...)
    if from !== nothing
        throw(ArgumentError("Use `usein` to specify `from` as a keyword argument."))
    end
    usein(downpath; from=upstream, kwargs...)
end

"""
    Rogue.add(name; project)

Install an unregistered package checked out at `~/.julia/dev/\$name`.
Its unregistered dependencies are installed using information stored
in `Manifest.toml` file checked in its repository (e.g.,
`~/.julia/dev/\$name/test/Manifest.toml`).

# Arguments
- `name::AbstractString`: Name of the package to be installed.

# Keyword Arguments
- `project::Union{Nothing, AbstractString} = nothing`: Project
  in which the package is installed.  `nothing` (default) means
  the current activated project.
- `prefer_https::Bool = false`: Prefer HTTPS repository URL rather
  than the one used in Git repository.
"""
function add(
    name::AbstractString;
    project::Union{Nothing, AbstractString} = nothing,
    prefer_https::Bool = false,
)
    if isdir(name)
        path = name
        name = basename(path)
    else
        path = joinpath(expanduser("~/.julia/dev"), name)  # TODO: don't
    end
    @assert isdir(path)

    function sortkey(p)
        r = relpath(p, path)
        rank = get(Dict(
            "" => 0,
            "test" => 1,
            "docs" => 2,
        ), dirname(r), typemax(Int))
        return (rank, r)
    end
    manifestfile, = sort(find_manifests(path), by=sortkey)
    spec = pkgspecof(path; prefer_https=prefer_https)
    deps = keys(get(
        TOML.parsefile(projecttomlpath(path)),
        "deps",
        Dict{String,String}(),
    ))

    @info "Installing private package from: $(spec.repo.url)"
    @info "Installing private dependencies from: $manifestfile"
    local private
    temporaryactivating(project) do
        private = add_private_projects_from(deps, manifestfile)
        Pkg.add(spec)
    end

    if isempty(private.deps)
        private_deps = ""
    else
        private_deps = """
        with private dependencies:
            $(join(private.deps, "\n    "))
        """
    end
    @info """
    Added:
        $name
    $private_deps
    """

    return
end

"""
    Rogue.addin(project; package=".", kwargs...)

A shortcut for `Rogue.add(package; project=project, kwargs...)`.
"""
addin(
    project::AbstractString,
    package::AbstractString = ".";
    kwargs...,
) = add(package, project=project; kwargs...)

"""
    Rogue.status()
    Rogue.status(project::String)
"""
function status(
    project::Union{Nothing, AbstractString} = nothing,
)

    if project === nothing
        project = dirname(Base.active_project())
    end
    manifest = TOML.parsefile(manifesttomlpath(project))
    foreach(private_projects_in_manifest(), manifest) do (name, entry)
        haskey(entry, "git-tree-sha1") || return  # possible?

        url = entry["repo-url"]
        treesha1 = entry["git-tree-sha1"]
        uuid = entry["uuid"]

        spec = Pkg.PackageSpec(
            name = name,
            uuid = uuid,
            url = url,
        )
        path = Pkg.Types.fresh_clone(spec)
        try
            open(git_cmd(`log --format=format:'%H %T'`, path)) do stream
                while !eof(stream)
                    ln = readline(stream)
                    logcommit, logtree = split(ln, " ", limit=2)
                    if logtree == treesha1
                        run(git_cmd(`log --max-count=1 --oneline $logcommit`, path))
                        return
                    end
                end
                @warn "Tree SHA `$treesha1` not found for package `$name` [$uuid]"
            end
        finally
            rm(path; force=true, recursive=true)
        end
    end
end


text:
export run_acdcopf_bf

""
function run_acdcopf_bf(file::String, model_constructor::AbstractBFModel, solver; kwargs...)
    data = PowerModels.parse_file(file)
    PowerModelsACDC.process_additional_data!(data)
    return run_acdcopf_bf(data, model_constructor, solver; kwargs...)
end

""
function run_acdcopf_bf(data::Dict{String,Any}, model_constructor::AbstractBFModel, solver; kwargs...)
    pm = PowerModels.build_model(data, model_constructor, post_acdcopf_bf; kwargs...)
    return PowerModels.optimize_model!(pm, solver; solution_builder = get_solution_acdc)
end

function post_acdcopf_bf(pm::AbstractPowerModel)
    add_ref_dcgrid!(pm)
    PowerModels.variable_voltage(pm)
    PowerModels.variable_generation(pm)
    PowerModels.variable_branch_flow(pm)
    PowerModels.variable_branch_current(pm)

    variable_active_dcbranch_flow(pm)
    variable_dcbranch_current(pm)
    variable_dc_converter(pm)
    variable_dcgrid_voltage_magnitude(pm)

    PowerModels.objective_min_fuel_cost(pm)

    PowerModels.constraint_model_current(pm)
    constraint_voltage_dc(pm)

    for i in PowerModels.ids(pm, :ref_buses)
        PowerModels.constraint_theta_ref(pm, i)
    end

    for i in PowerModels.ids(pm, :bus)
        constraint_kcl_shunt(pm, i)
    end

    for i in PowerModels.ids(pm, :branch)

        PowerModels.constraint_flow_losses(pm, i)
        PowerModels.constraint_voltage_magnitude_difference(pm, i)

        PowerModels.constraint_voltage_angle_difference(pm, i)
        PowerModels.constraint_thermal_limit_from(pm, i)
        PowerModels.constraint_thermal_limit_to(pm, i)
    end
    for i in PowerModels.ids(pm, :busdc)
        constraint_kcl_shunt_dcgrid(pm, i)
    end
    for i in PowerModels.ids(pm, :branchdc)
        constraint_ohms_dc_branch(pm, i)
        constraint_dc_branch_current(pm, i)
    end
    for i in PowerModels.ids(pm, :convdc)
        constraint_converter_losses(pm, i)
        constraint_converter_current(pm, i)
        constraint_conv_transformer(pm, i)
        constraint_conv_reactor(pm, i)
        constraint_conv_filter(pm, i)
        if pm.ref[:nw][pm.cnw][:convdc][i]["islcc"] == 1
            constraint_conv_firing_angle(pm, i)
        end
    end
end


text:
# Description: see script_mse_N2 but now we vary N and N2 while keeping
# their ratio fixed s.t. N2 = 1/alpha *N where alpha is the subsampling
# factor.
#
#
# Parameters:
# Nrange: vector of number of measurements
# N2: ambient dimensions
# alpha: subsampling factor
# K: number of sinusoids
# L: number of snapshots
# Niters: number of Monte Carlo iterations
# SNR_db: signal-to-noise-ratio [dB] of the AWGN model
# f_true: location of the frequencies of the sinusoids

include("MUSIC.jl")
include("rootMUSIC.jl")
include("ESPRIT.jl")
include("TLSesprit.jl")
include("MPM.jl")
include("ANM.jl")
include("EMAC.jl")
include("periodogram.jl")
include("findMinWrapDist.jl")
include("findMinWrapDistUnique.jl")
include("WrapDist.jl")

using Plots, Random

Nrange = ceil.(Int,range(6,30,length=5))
L = 1
SNR_db = 30
Niters = 5
alpha = 0.5
f_true = [0.3 0.32 0.34]
K = length(f_true)

MSE = zeros(length(Nrange),3)

A(f,N) = exp.(2*pi*1im*kron(f,(0:N-1)))
SNR = 10 .^(SNR_db./10)

for ii = 1:length(Nrange)
    N = Nrange[ii]
    N2::Int = ceil(1/alpha*N)
    for jj = 1:Niters

        noise = (randn(N,L)+1im*randn(N,L))/ (sqrt(2)*SNR)
        omega = sort(randperm(N2)[1:N])
        S = exp.(2*pi*rand(K,L)*1im)

        Y1 = A(f_true,N)*S + noise
        Y2 = A(f_true,N2)[omega,:]*S + noise

        t0 = time_ns()
        f_MUSIC = MPM(Y1,K)
        t1 = (time_ns()-t0)/1e9
        MinWrapDist = findMinWrapDistUnique(f_true,f_MUSIC)
        MSE[ii,1] += norm(MinWrapDist)^2

        f_EMAC = EMAC(Y2,N2,K,omega,1/SNR)
        MinWrapDist = findMinWrapDistUnique(f_true,f_EMAC)
        MSE[ii,2] += norm(MinWrapDist)^2

        f_ANM = ANM(Y2,N2,K,omega,1/SNR)
        MinWrapDist = findMinWrapDistUnique(f_true,f_ANM)
        MSE[ii,3] += norm(MinWrapDist)^2

    end
end

MSE /= Niters

plot(Nrange, MSE[:,1],yaxis=:log,tickfontsize = 6,
label="Matrix Pencil Method",
#xlims = (0,1),
#xticks = round.(f_range.-f1; digits=2),
xlabel = "N",
ylabel = "MSE"
)
plot!(Nrange,MSE[:,2],label="EMAC")
plot!(Nrange,MSE[:,3],label="ANM")
#savefig("/home/thomas/atom-amd64/usr/share/atom/example/mse_N_N2.png")


text:
# Deprecated functionality

using Base: @deprecate_binding

@deprecate_binding BLAS CUBLAS
@deprecate_binding FFT CUFFT

@deprecate cuzeros CuArrays.zeros
@deprecate cuones CuArrays.ones
@deprecate cufill CuArrays.fill


text:
# function [ret]=metadprimepm(h,f,Hp,Fp,Hm,Fm,sigma)
#
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# %  Code for computing meta-dprime-balance
# %
# %  Adam Barrett August 2012
# %
# % INPUTS:       h: type I HR
# %               f: type I FAR
# %               Hp: type II HR for positive type I responses
# %               Fp: type II FAR for positive type I responses
# %               Hm: type II HR for negative type I responses
# %               Fm: type II FAR for negative type I responses
# %               sigma: assumed ratio of signal standard deviations for stimulus B versus stimulus A (default value 1)
# %
# % OUTPUTS:      ret.dprime
# %               ret.metadprimebalance
# %               ret.stable: indicates stability according to equation (29),
# %                           taking the value 1 if this stability criterion is satisfied
# %                           and the value 0 if it is not satisfied
# %               ret.htp : meta type I HR arising from the type II data H_+ and F_+
# %               ret.ftp : meta type I FAR arising from the type II data H_+ and F_+
# %               ret.htm : meta type I HR arising from the type II data H_- and F_-
# %               ret.ftm : meta type I FAR arising from the type II data H_- and F_-
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#
# % Compute d' and relative type I threshold
# zf=icdf('Normal',f,0,1);
# zh=icdf('Normal',h,0,1);
# theta=-zf;
# dprime=sigma*zh-zf;
# thetaprime=theta/dprime;
#
# global th hp fp hm fm ssigma
# th=thetaprime;
# hp=Hp;
# fp=Fp;
# hm=Hm;
# fm=Fm;
# ssigma=sigma;
#
# % Solve equations for meta-d'_plus and meta-d'_minus
# x0=[theta;dprime];
# ep=fsolve(@eqformetadplus,x0);
# metadprimeplus=ep(2)
#
# x0=[thetaprime;dprime];
# em=fsolve(@eqformetadminus,x0);
# metadprimeminus=em(2)
#
# % Compute meta type I data
# htp=1-cdf('Normal',th*ep(2),ep(2),sigma);
# ftp=1-cdf('Normal',th*ep(2),0,1);
# htm=1-cdf('Normal',th*em(2),em(2),sigma);
# ftm=1-cdf('Normal',th*em(2),0,1);
#
# % Infer stability of output based on equation (29)
# stable=1;
# if htp>0.95 || ftp>0.95 || htm>0.95 || ftm>0.95 || htp<0.05 || ftp<0.05 || htm<0.05 || ftm<0.05
#     disp('Warning, unstable data');
#     stable=0;
# end
#
# % Compute proportion of positive type I responses for the weighting in
# % formula for meta-d'-balance
# r=(h+f)/2;
#
# ret.dprime=dprime;
# ret.metadprimebalance=r*metadprimeplus+(1-r)*metadprimeminus;
# ret.stable=stable;
# ret.htp=htp;
# ret.ftp=ftp;
# ret.htm=htm;
# ret.ftm=ftm;


text:
module Flatten

export flatten

"""
    flatten(x)

Turns a high-dimensional array (e.g., a batch of feature maps) into a 2-d array,
linearizing all except the last (batch) dimension.
"""
flatten(x) = reshape(x, :, size(x, ndims(x)))

end # module Flatten


text:
abstract type BuildTargets end
struct JuliaTarget <: BuildTargets end
struct StanTarget <: BuildTargets end
struct CTarget <: BuildTargets end
struct MATLABTarget <: BuildTargets end

abstract type ParallelForm end
struct SerialForm <: ParallelForm end
struct MultithreadedForm <: ParallelForm end
struct DistributedForm <: ParallelForm end
struct DaggerForm <: ParallelForm end

"""
`build_function`

Generates a numerically-usable function from a ModelingToolkit `Num`.

```julia
build_function(ex, args...;
               expression = Val{true},
               target = JuliaTarget(),
               kwargs...)
```

Arguments:

- `ex`: The `Num` to compile
- `args`: The arguments of the function
- `expression`: Whether to generate code or whether to generate the compiled form.
  By default, `expression = Val{true}`, which means that the code for the
  function is returned. If `Val{false}`, then the returned value is compiled.

Keyword Arguments:

- `target`: The output target of the compilation process. Possible options are:
    - `JuliaTarget`: Generates a Julia function
    - `CTarget`: Generates a C function
    - `StanTarget`: Generates a function for compiling with the Stan probabilistic
      programming language
    - `MATLABTarget`: Generates an anonymous function for use in MATLAB and Octave
      environments
- `fname`: Used by some targets for the name of the function in the target space.

Note that not all build targets support the full compilation interface. Check the
individual target documentation for details.
"""
function build_function(args...;target = JuliaTarget(),kwargs...)
  _build_function(target,args...;kwargs...)
end

function addheader(ex, fargs, iip; X=gensym(:MTIIPVar))
  if iip
    wrappedex = :(
        ($X,$(fargs.args...)) -> begin
        $ex
        nothing
      end
    )
  else
    wrappedex = :(
      ($(fargs.args...),) -> begin
        $ex
      end
    )
  end
  wrappedex
end

function add_integrator_header(ex, fargs, iip; X=gensym(:MTIIPVar))
  integrator = gensym(:MTKIntegrator)
  if iip
    wrappedex = :(
        $integrator -> begin
        ($X,$(fargs.args...)) = (($integrator).u,($integrator).u,($integrator).p,($integrator).t)
        $ex
        nothing
      end
    )
  else
    wrappedex = :(
      $integrator -> begin
      ($(fargs.args...),) = (($integrator).u,($integrator).p,($integrator).t)
        $ex
      end
    )
  end
  wrappedex
end

function unflatten_long_ops(op, N=4)
    rule1 = @rule((+)((~~x)) => length(~~x) > N ?
                 +(+((~~x)[1:N]...) + (+)((~~x)[N+1:end]...)) : nothing)
    rule2 = @rule((*)((~~x)) => length(~~x) > N ?
                 *(*((~~x)[1:N]...) * (*)((~~x)[N+1:end]...)) : nothing)

    op = value(op)
    Rewriters.Fixpoint(Rewriters.Postwalk(Rewriters.Chain([rule1, rule2])))(op)
end

struct Let
    eqs::Vector
    body
end

function observed_let(eqs)
    process -> ex -> begin
        isempty(eqs) && return ex

        assignments = map(eq -> :($(process(eq.lhs)) = $(process(eq.rhs))), eqs)
        letexpr = :(let $(assignments...)
                    end)
        # avoid a superfluous `begin ... end` block
        letexpr.args[2] = ex
        return letexpr
    end
end

function _build_function(target::JuliaTarget, op::Let, args...; conv=toexpr, kw...)
    _build_function(target, op.body, args...;
                    inner_let = observed_let(op.eqs), kw...)
end

# Scalar output
function _build_function(target::JuliaTarget, op, args...;
                         conv = toexpr, expression = Val{true},
                         checkbounds = false,
                         inner_let = nothing,
                         linenumbers = true, headerfun=addheader)

    argnames = [gensym(:MTKArg) for i in 1:length(args)]
    symsdict = Dict()
    arg_pairs = map((x,y)->vars_to_pairs(x,y, symsdict), argnames, args)
    process = unflatten_long_ops∘(x->substitute(x, symsdict, fold=false))
    ls = reduce(vcat,conv.(first.(arg_pairs)))
    rs = reduce(vcat,last.(arg_pairs))
    var_eqs = Expr(:(=), build_expr(:tuple, ls), build_expr(:tuple, conv.(process.(rs))))

    fname = gensym(:ModelingToolkitFunction)
    op = process(op)
    out_expr = conv(substitute(op, symsdict, fold=false))

    if inner_let !== nothing
        inner_let_expr = inner_let(conv ∘ process)
        out_expr = inner_let_expr(out_expr)
    end

    let_expr = Expr(:let, var_eqs, Expr(:block, out_expr))
    bounds_block = checkbounds ? let_expr : :(@inbounds begin $let_expr end)

    fargs = Expr(:tuple,argnames...)
    oop_ex = headerfun(bounds_block, fargs, false)

    if !linenumbers
        oop_ex = striplines(oop_ex)
    end

    if expression == Val{true}
        return ModelingToolkit.inject_registered_module_functions(oop_ex)
    else
        _build_and_inject_function(@__MODULE__, oop_ex)
    end
end

function _build_and_inject_function(mod::Module, ex)
    @RuntimeGeneratedFunction(ModelingToolkit.inject_registered_module_functions(ex))
end

# Detect heterogeneous element types of "arrays of matrices/sparce matrices"
function is_array_matrix(F)
    return isa(F, AbstractVector) && all(x->isa(x, AbstractArray), F)
end
function is_array_sparse_matrix(F)
    return isa(F, AbstractVector) && all(x->isa(x, AbstractSparseMatrix), F)
end
# Detect heterogeneous element types of "arrays of arrays of matrices/sparce matrices"
function is_array_array_matrix(F)
    return isa(F, AbstractVector) && all(x->isa(x, AbstractArray{<:AbstractMatrix}), F)
end
function is_array_array_sparse_matrix(F)
    return isa(F, AbstractVector) && all(x->isa(x, AbstractArray{<:AbstractSparseMatrix}), F)
end

function fill_array_with_zero!(x::AbstractArray)
    if eltype(x) <: AbstractArray
        foreach(fill_array_with_zero!, x)
    else
        fill!(x, false)
    end
    return x
end

"""
Build function target: JuliaTarget

```julia
function _build_function(target::JuliaTarget, rhss, args...;
                         conv = toexpr, expression = Val{true},
                         checkbounds = false,
                         linenumbers = false, multithread=nothing,
                         headerfun = addheader, outputidxs=nothing,
                         convert_oop = true, force_SA = false,
                         skipzeros = outputidxs===nothing,
                         fillzeros = skipzeros && !(typeof(rhss)<:SparseMatrixCSC),
                         parallel=SerialForm(), kwargs...)
```

Generates a Julia function which can then be utilized for further evaluations.
If expression=Val{false}, the return is a Julia function which utilizes
RuntimeGeneratedFunctions.jl in order to be free of world-age issues.

If the `Num` is an `Operation`, the generated function is a function
with a scalar output, otherwise if it's an `AbstractArray{Operation}`, the output
is two functions, one for out-of-place AbstractArray output and a second which
is a mutating function. The outputted functions match the given argument order,
i.e., f(u,p,args...) for the out-of-place and scalar functions and
`f!(du,u,p,args..)` for the in-place version.

Special Keyword Argumnets:

- `parallel`: The kind of parallelism to use in the generated function. Defaults
  to `SerialForm()`, i.e. no parallelism. Note that the parallel forms are not
  exported and thus need to be chosen like `ModelingToolkit.SerialForm()`.
  The choices are:
  - `SerialForm()`: Serial execution.
  - `MultithreadedForm()`: Multithreaded execution with a static split, evenly
    splitting the number of expressions per thread.
  - `DistributedForm()`: Multiprocessing using Julia's Distributed with a static
    schedule, evenly splitting the number of expressions per process.
  - `DaggerForm()`: Multithreading and multiprocessing using Julia's Dagger.jl
    for dynamic scheduling and load balancing.
- `conv`: The conversion function of the Operation to Expr. By default this uses
  the `toexpr` function.
- `checkbounds`: For whether to enable bounds checking inside of the generated
  function. Defaults to false, meaning that `@inbounds` is applied.
- `linenumbers`: Determines whether the generated function expression retains
  the line numbers. Defaults to true.
- `convert_oop`: Determines whether the OOP version should try to convert
  the output to match the type of the first input. This is useful for
  cases like LabelledArrays or other array types that carry extra
  information. Defaults to true.
- `force_SA`: Forces the output of the OOP version to be a StaticArray.
  Defaults to `false`, and outputs a static array when the first argument
  is a static array.
- `skipzeros`: Whether to skip filling zeros in the in-place version if the
  filling function is 0.
- `fillzeros`: Whether to perform `fill(out,0)` before the calculations to ensure
  safety with `skipzeros`.
  """
  function _build_function(target::JuliaTarget, rhss::AbstractArray, args...;
                           conv = toexpr, expression = Val{true},
                           inner_let = nothing,
                           checkbounds = false,
                           linenumbers = false, multithread=nothing,
                           headerfun = addheader, outputidxs=nothing,
                           convert_oop = true, force_SA = false,
                           skipzeros = outputidxs===nothing,
                           fillzeros = skipzeros && !(typeof(rhss)<:SparseMatrixCSC),
                           parallel=SerialForm(), kwargs...)
      if multithread isa Bool
          @warn("multithraded is deprecated for the parallel argument. See the documentation.")
          parallel = multithread ? MultithreadedForm() : SerialForm()
      end

      argnames = [gensym(:MTKArg) for i in 1:length(args)]
      symsdict = Dict()
      arg_pairs = map((x,y)->vars_to_pairs(x,y, symsdict), argnames, args)
      process = unflatten_long_ops∘(x->substitute(x, symsdict, fold=false))

      ls = reduce(vcat,conv.(first.(arg_pairs)))
      rs = reduce(vcat,last.(arg_pairs))
      var_eqs = Expr(:(=), ModelingToolkit.build_expr(:tuple, ls), ModelingToolkit.build_expr(:tuple, conv.(process.(rs))))

      fname = gensym(:ModelingToolkitFunction)
      fargs = Expr(:tuple,argnames...)

      oidx = isnothing(outputidxs) ? (i -> i) : (i -> outputidxs[i])
      X = gensym(:MTIIPVar)

      if rhss isa SparseMatrixCSC
          rhs_length = length(rhss.nzval)
          rhss = SparseMatrixCSC(rhss.m, rhss.m, rhss.colptr, rhss.rowval, map(process, rhss.nzval))
      else
          rhs_length = length(rhss)
          rhss = [process(r) for r in rhss]
      end

      if parallel isa DistributedForm
          numworks = Distributed.nworkers()
          reducevars = [gensym(:MTReduceVar) for i in 1:numworks]
          lens = Int(ceil(rhs_length/numworks))
          finalsize = rhs_length - (numworks-1)*lens
          _rhss = vcat(reduce(vcat,[[Variable(reducevars[i],j) for j in 1:lens] for i in 1:numworks-1],init=Expr[]),
                       [Variable(reducevars[end],j) for j in 1:finalsize])

      elseif parallel isa DaggerForm
          computevars = [gensym(:MTComputeVar) for i in axes(rhss,1)]
          reducevar = Variable(gensym(:MTReduceVar))
          _rhss = [Variable(reducevar,i) for i in axes(rhss,1)]
      elseif rhss isa SparseMatrixCSC
          _rhss = rhss.nzval
      else
          _rhss = rhss
      end

    ip_sys_exprs = Expr[]
    # we cannot reliably fill the array with the presence of index translation
    if is_array_array_sparse_matrix(rhss) # Array of arrays of sparse matrices
        for (i, rhsel) ∈ enumerate(_rhss)
            for (j, rhsel2) ∈ enumerate(rhsel)
                for (k, rhs) ∈ enumerate(rhsel2.nzval)
                    rhs′ = conv(rhs)
                    (skipzeros && rhs′ isa Number && iszero(rhs′)) && continue
                    push!(ip_sys_exprs, :($X[$i][$j].nzval[$k] = $rhs′))
                end
            end
        end
    elseif is_array_array_matrix(rhss) # Array of arrays of arrays
        for (i, rhsel) ∈ enumerate(_rhss)
            for (j, rhsel2) ∈ enumerate(rhsel)
                for (k, rhs) ∈ enumerate(rhsel2)
                    rhs′ = conv(rhs)
                    (skipzeros && rhs′ isa Number && iszero(rhs′)) && continue
                    push!(ip_sys_exprs, :($X[$i][$j][$k] = $rhs′))
                end
            end
        end
    elseif is_array_sparse_matrix(rhss) # Array of sparse matrices
        for (i, rhsel) ∈ enumerate(_rhss)
            for (j, rhs) ∈ enumerate(rhsel.nzval)
                rhs′ = conv(rhs)
                (skipzeros && rhs′ isa Number && iszero(rhs′)) && continue
                push!(ip_sys_exprs, :($X[$i].nzval[$j] = $rhs′))
            end
        end
    elseif is_array_matrix(rhss) # Array of arrays
        for (i, rhsel) ∈ enumerate(_rhss)
            for (j, rhs) ∈ enumerate(rhsel)
                rhs′ = conv(rhs)
                (skipzeros && rhs′ isa Number && iszero(rhs′)) && continue
                push!(ip_sys_exprs, :($X[$i][$j] = $rhs′))
            end
        end
    elseif rhss isa SparseMatrixCSC
        for (i, rhs) ∈ enumerate(_rhss)
            rhs′ = conv(rhs)
            (skipzeros && rhs′ isa Number && iszero(rhs′)) && continue
            push!(ip_sys_exprs, :($X.nzval[$i] = $rhs′))
        end
    else
        for (i, rhs) ∈ enumerate(_rhss)
            rhs′ = conv(rhs)
            (skipzeros && rhs′ isa Number && iszero(rhs′)) && continue
            push!(ip_sys_exprs, :($X[$(oidx(i))] = $rhs′))
        end
    end

    ip_let_expr = Expr(:let, var_eqs, build_expr(:block, ip_sys_exprs))

    if parallel isa MultithreadedForm
        lens = Int(ceil(length(ip_let_expr.args[2].args)/Threads.nthreads()))
        threaded_exprs = vcat([quote
           Threads.@spawn begin
              $(ip_let_expr.args[2].args[((i-1)*lens+1):i*lens]...)
           end
        end for i in 1:Threads.nthreads()-1],
           quote
              Threads.@spawn begin
                  $(ip_let_expr.args[2].args[((Threads.nthreads()-1)*lens+1):end]...)
              end
          end)
        ip_let_expr.args[2] =  ModelingToolkit.build_expr(:block, threaded_exprs)
        ip_let_expr = :(@sync begin $ip_let_expr end)
    elseif parallel isa DistributedForm
        numworks = Distributed.nworkers()
        lens = Int(ceil(length(ip_let_expr.args[2].args)/numworks))
        spawnvars = [gensym(:MTSpawnVar) for i in 1:numworks]
        rhss_flat = rhss isa SparseMatrixCSC ? rhss.nzval : rhss
        spawnvectors = vcat(
                            [build_expr(:vect, [conv(rhs) for rhs ∈ rhss_flat[((i-1)*lens+1):i*lens]]) for i in 1:numworks-1],
                            build_expr(:vect, [conv(rhs) for rhs ∈ rhss_flat[((numworks-1)*lens+1):end]]))

        spawn_exprs = [quote
                           $(spawnvars[i]) = ModelingToolkit.Distributed.remotecall($(i+1)) do
                               $(spawnvectors[i])
                           end
                       end for i in 1:numworks]
        spawn_exprs = ModelingToolkit.build_expr(:block, spawn_exprs)
        resunpack_exprs = [:($(Symbol(reducevars[iter])) = fetch($(spawnvars[iter]))) for iter in 1:numworks]

        ip_let_expr.args[2] = quote
            @sync begin
                $spawn_exprs
                $(resunpack_exprs...)
                $(ip_let_expr.args[2])
            end
        end
    elseif parallel isa DaggerForm
        @assert HAS_DAGGER[] "Dagger.jl is not loaded; please do `using Dagger`"
        dagwrap(x) = x
        dagwrap(ex::Expr) = dagwrap(ex, Val(ex.head))
        dagwrap(ex::Expr, ::Val) = ex
        dagwrap(ex::Expr, ::Val{:call}) = :(Dagger.delayed($(ex.args[1]))($(dagwrap.(ex.args[2:end])...)))
        new_rhss = dagwrap.(conv.(rhss))
        delayed_exprs = build_expr(:block, [:($(Symbol(computevars[i])) = Dagger.delayed(identity)($(new_rhss[i]))) for i in axes(computevars,1)])
        # TODO: treereduce?
        reduce_expr = quote
            $(Symbol(reducevar)) = collect(Dagger.delayed(vcat)($(computevars...)))
        end
        ip_let_expr.args[2] = quote
            @sync begin
                $delayed_exprs
                $reduce_expr
                $(ip_let_expr.args[2])
            end
        end
    end

    if rhss isa SparseMatrixCSC
        rhss′ = map(conv∘process, rhss.nzval)
    else
        rhss′ = [conv(process(r)) for r in rhss]
    end

    tuple_sys_expr = build_expr(:tuple, rhss′)

    if rhss isa Matrix
        arr_sys_expr = build_expr(:vcat, [build_expr(:row,[conv(rhs) for rhs ∈ rhss[i,:]]) for i in 1:size(rhss,1)])
    elseif typeof(rhss) <: Array && !(typeof(rhss) <: Vector)
        vector_form = build_expr(:vect, [conv(rhs) for rhs ∈ rhss])
        arr_sys_expr = :(reshape($vector_form,$(size(rhss)...)))
    elseif rhss isa SparseMatrixCSC
        vector_form = build_expr(:vect, [conv(rhs) for rhs ∈ nonzeros(rhss)])
        arr_sys_expr = :(SparseMatrixCSC{eltype($(first(argnames))),Int}($(size(rhss)...), $(rhss.colptr), $(rhss.rowval), $vector_form))
    else # Vector
        arr_sys_expr = build_expr(:vect, [conv(rhs) for rhs ∈ rhss])
    end

    xname = gensym(:MTK)

    arr_sys_expr = (typeof(rhss) <: Vector || typeof(rhss) <: Matrix) && !(eltype(rhss) <: AbstractArray) ? quote
        if $force_SA || typeof($(fargs.args[1])) <: Union{ModelingToolkit.StaticArrays.SArray,ModelingToolkit.LabelledArrays.SLArray}
            $xname = ModelingToolkit.StaticArrays.@SArray $arr_sys_expr
            if $convert_oop && !(typeof($(fargs.args[1])) <: Number) && $(typeof(rhss) <: Vector) # Only try converting if it should match `u`
                return similar_type($(fargs.args[1]),eltype($xname))($xname)
            else
                return $xname
            end
        else
            $xname = $arr_sys_expr
            if $convert_oop && $(typeof(rhss) <: Vector)
                if !(typeof($(fargs.args[1])) <: Array) && !(typeof($(fargs.args[1])) <: Number) && eltype($(fargs.args[1])) <: eltype($xname)
                    # Last condition: avoid known error because this doesn't change eltypes!
                    return convert(typeof($(fargs.args[1])),$xname)
                elseif typeof($(fargs.args[1])) <: ModelingToolkit.LabelledArrays.LArray
                    # LArray just needs to add the names back!
                    return ModelingToolkit.LabelledArrays.LArray{ModelingToolkit.LabelledArrays.symnames(typeof($(fargs.args[1])))}($xname)
                else
                    return $xname
                end
            else
                return $xname
            end
        end
    end : arr_sys_expr

    if inner_let !== nothing
        inner_let_expr = inner_let(conv ∘ process)
        arr_sys_expr = inner_let_expr(arr_sys_expr)
        ip_let_expr.args[2] = inner_let_expr(ip_let_expr.args[2])
    end

    if fillzeros && outputidxs === nothing
        ip_let_expr = quote
            $fill_array_with_zero!($X)
            $ip_let_expr
        end
    end

    arr_let_expr = Expr(:let, var_eqs, arr_sys_expr)

    oop_bounds_block = checkbounds ? arr_let_expr : :(@inbounds $arr_let_expr)
    ip_bounds_block = checkbounds ? ip_let_expr : :(@inbounds $ip_let_expr)

    oop_ex = headerfun(oop_bounds_block, fargs, false)
    iip_ex = headerfun(ip_bounds_block, fargs, true; X=X)

    if !linenumbers
        oop_ex = striplines(oop_ex)
        iip_ex = striplines(iip_ex)
    end

    if expression == Val{true}
        return ModelingToolkit.inject_registered_module_functions(oop_ex), ModelingToolkit.inject_registered_module_functions(iip_ex)
    else
        return _build_and_inject_function(@__MODULE__, oop_ex), _build_and_inject_function(@__MODULE__, iip_ex)
    end
end

function vars_to_pairs(name,vs::Union{Tuple, AbstractArray}, symsdict=Dict())
    vs_names = tosymbol.(vs)
    for (v,k) in zip(vs_names, vs)
        symsdict[k] = Sym{symtype(k)}(v)
    end
    exs = [:($name[$i]) for (i, u) ∈ enumerate(vs)]
    vs_names,exs
end
function vars_to_pairs(name,vs, symsdict)
    symsdict[vs] = Sym{symtype(vs)}(tosymbol(vs))
    [tosymbol(vs)], [name]
end

get_varnumber(varop, vars::Vector) =  findfirst(x->isequal(x,varop),vars)

function numbered_expr(O::Symbolic,args...;varordering = args[1],offset = 0,
                       lhsname=gensym("du"),rhsnames=[gensym("MTK") for i in 1:length(args)])
    O = value(O)
    if O isa Sym || isa(operation(O), Sym)
        for j in 1:length(args)
            i = get_varnumber(O,args[j])
            if i !== nothing
                return :($(rhsnames[j])[$(i+offset)])
            end
        end
    end
  return Expr(:call, O isa Sym ? tosymbol(O, escape=false) : Symbol(operation(O)),
         [numbered_expr(x,args...;offset=offset,lhsname=lhsname,
                        rhsnames=rhsnames,varordering=varordering) for x in arguments(O)]...)
end

function numbered_expr(de::ModelingToolkit.Equation,args...;varordering = args[1],
                       lhsname=gensym("du"),rhsnames=[gensym("MTK") for i in 1:length(args)],offset=0)

    varordering = value.(args[1])
    var = var_from_nested_derivative(de.lhs)[1]
    i = findfirst(x->isequal(tosymbol(x isa Sym ? x : operation(x), escape=false), tosymbol(var, escape=false)),varordering)
    :($lhsname[$(i+offset)] = $(numbered_expr(de.rhs,args...;offset=offset,
                                              varordering = varordering,
                                              lhsname = lhsname,
                                              rhsnames = rhsnames)))
end
numbered_expr(c,args...;kwargs...) = c
numbered_expr(c::Num,args...;kwargs...) = error("Num found")

"""
Build function target: CTarget

```julia
function _build_function(target::CTarget, eqs::Array{<:Equation}, args...;
                         conv = toexpr, expression = Val{true},
                         fname = :diffeqf,
                         lhsname=:du,rhsnames=[Symbol("RHS\$i") for i in 1:length(args)],
                         libpath=tempname(),compiler=:gcc)
```

This builds an in-place C function. Only works on arrays of equations. If
`expression == Val{false}`, then this builds a function in C, compiles it,
and returns a lambda to that compiled function. These special keyword arguments
control the compilation:

- libpath: the path to store the binary. Defaults to a temporary path.
- compiler: which C compiler to use. Defaults to :gcc, which is currently the
  only available option.
"""
function _build_function(target::CTarget, eqs::Array{<:Equation}, args...;
                         conv = toexpr, expression = Val{true},
                         fname = :diffeqf,
                         lhsname=:du,rhsnames=[Symbol("RHS$i") for i in 1:length(args)],
                         libpath=tempname(),compiler=:gcc)

    differential_equation = string(join([numbered_expr(eq,args...,lhsname=lhsname,
                                  rhsnames=rhsnames,offset=-1) for
                                  (i, eq) ∈ enumerate(eqs)],";\n  "),";")

    argstrs = join(vcat("double* $(lhsname)",[typeof(args[i])<:Array ? "double* $(rhsnames[i])" : "double $(rhsnames[i])" for i in 1:length(args)]),", ")
    ex = """
    void $fname($(argstrs...)) {
      $differential_equation
    }
    """

    if expression == Val{true}
        return ex
    else
        @assert compiler == :gcc
        ex = build_function(eqs,args...;target=ModelingToolkit.CTarget())
        open(`gcc -fPIC -O3 -msse3 -xc -shared -o $(libpath * "." * Libdl.dlext) -`, "w") do f
            print(f, ex)
        end
        @RuntimeGeneratedFunction(:((du::Array{Float64},u::Array{Float64},p::Array{Float64},t::Float64) -> ccall(("diffeqf", $libpath), Cvoid, (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Float64), du, u, p, t)))
    end
end

"""
Build function target: StanTarget

```julia
function _build_function(target::StanTarget, eqs::Array{<:Equation}, vs, ps, iv;
                         conv = toexpr, expression = Val{true},
                         fname = :diffeqf, lhsname=:internal_var___du,
                         rhsnames=[:internal_var___u,:internal_var___p,:internal_var___t])
```

This builds an in-place Stan function compatible with the Stan differential equation solvers.
Unlike other build targets, this one requestions (vs, ps, iv) as the function arguments.
Only allowed on arrays of equations.
"""
function _build_function(target::StanTarget, eqs::Array{<:Equation}, vs, ps, iv;
                         conv = toexpr, expression = Val{true},
                         fname = :diffeqf, lhsname=:internal_var___du,
                         rhsnames=[:internal_var___u,:internal_var___p,:internal_var___t])
    @assert expression == Val{true}
    differential_equation = string(join([numbered_expr(eq,vs,ps,lhsname=lhsname,
                                   rhsnames=rhsnames) for
                                   (i, eq) ∈ enumerate(eqs)],";\n  "),";")
    """
    real[] $fname(real $(conv(iv)),real[] $(rhsnames[1]),real[] $(rhsnames[2]),real[] x_r,int[] x_i) {
      real $lhsname[$(length(eqs))];
      $differential_equation
      return $lhsname;
    }
    """
end

"""
Build function target: MATLABTarget

```julia
function _build_function(target::MATLABTarget, eqs::Array{<:Equation}, args...;
                         conv = toexpr, expression = Val{true},
                         lhsname=:internal_var___du,
                         rhsnames=[:internal_var___u,:internal_var___p,:internal_var___t])
```

This builds an out of place anonymous function @(t,rhsnames[1]) to be used in MATLAB.
Compatible with the MATLAB differential equation solvers. Only allowed on arrays
of equations.
"""
function _build_function(target::MATLABTarget, eqs::Array{<:Equation}, args...;
                         conv = toexpr, expression = Val{true},
                         fname = :diffeqf, lhsname=:internal_var___du,
                         rhsnames=[:internal_var___u,:internal_var___p,:internal_var___t])
    @assert expression == Val{true}
    matstr = join([numbered_expr(eq.rhs,args...,lhsname=lhsname,
                                  rhsnames=rhsnames) for
                                  (i, eq) ∈ enumerate(eqs)],"; ")

    matstr = replace(matstr,"["=>"(")
    matstr = replace(matstr,"]"=>")")
    matstr = "$fname = @(t,$(rhsnames[1])) ["*matstr*"];"
    matstr
end


text:
using CUDA
using IndexFunArrays

function f(arr, x)
    arr .= arr.^2 .+ arr.^3 .+ sqrt.(x)
    return arr
end

function f2(arr)
    arr .= arr.^2 .+ arr.^3 .+ sqrt.(arr .^2)
    return arr
end

function f3(arr)
    arr = arr.^2 .+ arr.^3 .+ sqrt.(arr .^2)
    return arr
end

function test(s)
    arr = randn(Float32, s)
    arr_c = CuArray(arr)
    x = rr2(Float32, s)
    x_c = CuArray(rr(Float32, s))

    @info "f on CPU"
    @time f(arr, x);
    @time f(arr, x);
    @info "f2 on CPU"
    @time f2(arr);
    @time f2(arr);
    @info "f on GPU"
    CUDA.@time f(arr_c, x);
    CUDA.@time f(arr_c, x);
    @info "f2 on GPU"
    CUDA.@time f2(arr_c);
    CUDA.@time f2(arr_c);
    @info "f3 on GPU"
    CUDA.@time f3(arr_c);
    CUDA.@time f3(arr_c);

    return
end

function MandelbrotIterations(c, Niter=500)
    x = c
    res = 0
    for n in 1:Niter
        x = x.^2 .+ c
        #if abs(x) > 4.0 && (res == 0)
        #    res = n
        #end
        res += (res == 0)*(abs(x)> 4.0)*n  # ensure that only the first breakage counts but avoid if
    end
    return (res % 127)
end

function MandelbrotIterations(idx,start,scale, Niter=100)
    c = Complex(start[1].+idx[1].*scale[1],start[2].+idx[2].*scale[2])
    return MandelbrotIterations(c..., Niter)
end

function testMandelbrotArray(s=(1024,1024), Niter=100)
    arr = zeros(Int64, s)
    scale = (3.0,3.0)./s;
    start=(-2.3,-1.5);
    @info "Test MandelBrot Array CPU"
    @time arr .= MandelbrotIterations.(CartesianIndices(arr), [scale], [start])
    @time arr .= MandelbrotIterations.(CartesianIndices(arr), [scale], [start])
end

function testMandelbrotIFA(s=(1024,1024), Niter=100)
    arr = zeros(Int64, s)
    scale = (3.0,3.0)./s;
    start=(-2.3,-1.5);
    Mandelbrot = IndexFunArray(x -> MandelbrotIterations(x,start, scale,Niter), size(arr));
    @info "Test MandelBrot IndexFunArray CPU"
    @time arr .= Mandelbrot;
    @time arr .= Mandelbrot;

    #arr_c = CuArray(arr)
    #CUDA.@time arr_c .= Mandelbrot
    #CUDA.@time arr_c .= Mandelbrot
    #CUDA.@time arr_c .= MandelbrotIterations.(arr_c,start,scale,Niter)
    return arr
end

#using Napari
#napari.view_image(testMandelbrotIFA())