text:
function _precompile_()
  ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
  precompile(
    Tuple{
      Core.kwftype(typeof(ITensors.eigen)),
      NamedTuple{
        (
          :ishermitian,
          :which_decomp,
          :tags,
          :maxdim,
          :mindim,
          :cutoff,
          :eigen_perturbation,
          :ortho,
        ),
        Tuple{Bool,Nothing,TagSet,Int64,Int64,Float64,Nothing,String},
      },
      typeof(eigen),
      ITensor{2},
      IndexSet{1,Index{Int64},Tuple{Index{Int64}}},
      IndexSet{1,Index{Int64},Tuple{Index{Int64}}},
    },
  )
  precompile(
    Tuple{
      Core.kwftype(typeof(ITensors.factorize)),
      NamedTuple{
        (:which_decomp, :tags, :maxdim, :mindim, :cutoff, :eigen_perturbation, :ortho),
        Tuple{Nothing,TagSet,Int64,Int64,Float64,Nothing,String},
      },
      typeof(factorize),
      ITensor{3},
      IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
    },
  )
  precompile(
    Tuple{
      Core.kwftype(typeof(ITensors.replacebond!)),
      NamedTuple{
        (:maxdim, :mindim, :cutoff, :eigen_perturbation, :ortho, :which_decomp),
        Tuple{Int64,Int64,Float64,Nothing,String,Nothing},
      },
      typeof(replacebond!),
      MPS,
      Int64,
      ITensor{3},
    },
  )
  precompile(
    Tuple{
      Core.kwftype(typeof(KrylovKit.eigsolve)),
      NamedTuple{
        (:ishermitian, :tol, :krylovdim, :maxiter),Tuple{Bool,Float64,Int64,Int64}
      },
      typeof(KrylovKit.eigsolve),
      ProjMPO,
      ITensor{3},
      Int64,
      Symbol,
    },
  )
  precompile(
    Tuple{
      Core.kwftype(typeof(KrylovKit.expand!)),
      NamedTuple{(:verbosity,),Tuple{Int64}},
      typeof(KrylovKit.expand!),
      KrylovKit.LanczosIterator{ProjMPO,ITensor{3},KrylovKit.ModifiedGramSchmidt2},
      KrylovKit.LanczosFactorization{ITensor{3},Float64},
    },
  )
  precompile(
    Tuple{
      Core.kwftype(typeof(NDTensors.eigen)),
      NamedTuple{
        (
          :ishermitian,
          :which_decomp,
          :tags,
          :maxdim,
          :mindim,
          :cutoff,
          :eigen_perturbation,
          :ortho,
        ),
        Tuple{Bool,Nothing,TagSet,Int64,Int64,Float64,Nothing,String},
      },
      typeof(eigen),
      LinearAlgebra.Hermitian{
        Float64,
        NDTensors.Tensor{
          Float64,
          2,
          NDTensors.Dense{Float64,Array{Float64,1}},
          IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
        },
      },
    },
  )
  precompile(
    Tuple{
      Core.kwftype(typeof(NDTensors.svd)),
      NamedTuple{(:alg,),Tuple{String}},
      typeof(svd),
      NDTensors.Tensor{
        Complex{Float64},
        2,
        NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
    },
  )
  precompile(
    Tuple{
      Core.kwftype(typeof(NDTensors.svd)),
      NamedTuple{(:alg,),Tuple{String}},
      typeof(svd),
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
    },
  )
  precompile(Tuple{ProjMPO,ITensor{3}})
  precompile(Tuple{Type{MPO},AutoMPO,Array{Index{Int64},1}})
  precompile(Tuple{Type{TagSet},String})
  precompile(
    Tuple{
      typeof(*),
      Float64,
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
    },
  )
  precompile(Tuple{typeof(*),ITensor{2},ITensor{3}})
  precompile(Tuple{typeof(*),ITensor{3},ITensor{1}})
  precompile(Tuple{typeof(*),ITensor{3},ITensor{2}})
  precompile(Tuple{typeof(*),ITensor{3},ITensor{3}})
  precompile(Tuple{typeof(*),ITensor{4},ITensor{1}})
  precompile(Tuple{typeof(*),ITensor{4},ITensor{4}})
  precompile(Tuple{typeof(==),Array{Index{Int64},1},Array{Index{Int64},1}})
  precompile(
    Tuple{
      typeof(==),
      IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
    },
  )
  precompile(Tuple{typeof(==),TagSet,TagSet})
  precompile(
    Tuple{
      typeof(Base.Broadcast.materialize!),
      ITensor{4},
      Base.Broadcast.Broadcasted{
        ITensors.ITensorStyle,Nothing,typeof(+),Tuple{ITensor{4},ITensor{4}}
      },
    },
  )
  precompile(Tuple{typeof(Base.afoldl),typeof(*),ITensor{4},ITensor{3}})
  precompile(
    Tuple{typeof(ITensors.setstore!),ITensor{2},NDTensors.Dense{Float64,Array{Float64,1}}}
  )
  precompile(
    Tuple{
      typeof(KrylovKit.orthogonalize!),ITensor{3},ITensor{3},KrylovKit.ModifiedGramSchmidt2
    },
  )
  precompile(
    Tuple{
      typeof(NDTensors.contract!!),
      NDTensors.Tensor{
        Complex{Float64},
        2,
        NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Complex{Float64},
        2,
        NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Complex{Float64},
        2,
        NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(NDTensors.contract!!),
      NDTensors.Tensor{
        Complex{Float64},
        2,
        NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Complex{Float64},
        2,
        NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Number,
        2,
        NDTensors.Combiner,
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(NDTensors.contract!!),
      NDTensors.Tensor{
        Complex{Float64},
        2,
        NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Diag{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Complex{Float64},
        2,
        NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(NDTensors.contract!!),
      NDTensors.Tensor{
        Float64,0,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{0,Union{},Tuple{}}
      },
      Tuple{},
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(NDTensors.contract!!),
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(NDTensors.contract!!),
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Number,
        2,
        NDTensors.Combiner,
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(NDTensors.contract!!),
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Diag{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(NDTensors.contract!!),
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
      NDTensors.Tensor{
        Float64,
        1,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{1,Index{Int64},Tuple{Index{Int64}}},
      },
      Tuple{Int64},
    },
  )
  precompile(
    Tuple{
      typeof(NDTensors.contract!!),
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(NDTensors.contract!!),
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(NDTensors.contract!!),
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(NDTensors.contract!!),
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
      NDTensors.Tensor{
        Float64,
        4,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}},
      },
      NTuple{4,Int64},
      NDTensors.Tensor{
        Float64,
        1,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{1,Index{Int64},Tuple{Index{Int64}}},
      },
      Tuple{Int64},
    },
  )
  precompile(
    Tuple{
      typeof(NDTensors.contract!!),
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
      NDTensors.Tensor{
        Float64,
        4,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}},
      },
      NTuple{4,Int64},
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(NDTensors.contract!!),
      NDTensors.Tensor{
        Float64,
        4,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}},
      },
      NTuple{4,Int64},
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(NDTensors.contract!!),
      NDTensors.Tensor{
        Float64,
        4,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}},
      },
      NTuple{4,Int64},
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(NDTensors.contract!!),
      NDTensors.Tensor{
        Float64,
        4,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}},
      },
      NTuple{4,Int64},
      NDTensors.Tensor{
        Float64,
        4,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}},
      },
      NTuple{4,Int64},
      NDTensors.Tensor{
        Float64,
        4,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}},
      },
      NTuple{4,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(NDTensors.contraction_output),
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      IndexSet{0,Union{},Tuple{}},
    },
  )
  precompile(
    Tuple{
      typeof(NDTensors.contraction_output),
      NDTensors.Tensor{
        Float64,
        4,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}},
      },
      NDTensors.Tensor{
        Float64,
        1,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{1,Index{Int64},Tuple{Index{Int64}}},
      },
      IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
    },
  )
  precompile(
    Tuple{
      typeof(NDTensors.getperm),
      IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}},
      IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}},
    },
  )
  precompile(
    Tuple{
      typeof(NDTensors.permutedims!!),
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
      Function,
    },
  )
  precompile(
    Tuple{
      typeof(NDTensors.permutedims!!),
      NDTensors.Tensor{
        Float64,
        4,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}},
      },
      NDTensors.Tensor{
        Float64,
        4,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}},
      },
      NTuple{4,Int64},
      Function,
    },
  )
  precompile(
    Tuple{
      typeof(Random.randn!),
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
    },
  )
  precompile(Tuple{typeof(axpy!),Float64,ITensor{3},ITensor{3}})
  precompile(Tuple{typeof(combiner),Index{Int64}})
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Complex{Float64},
        2,
        NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Complex{Float64},
        2,
        NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      Tuple{Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Complex{Float64},
        2,
        NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Complex{Float64},
        2,
        NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Complex{Float64},
        2,
        NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Number,
        2,
        NDTensors.Combiner,
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NTuple{4,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      Tuple{Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
      Tuple{Int64,Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Number,
        2,
        NDTensors.Combiner,
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      Tuple{Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Number,
        2,
        NDTensors.Combiner,
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Diag{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Complex{Float64},
        2,
        NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      Tuple{Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Diag{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Complex{Float64},
        2,
        NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Diag{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      Tuple{Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Diag{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
      NDTensors.Tensor{
        Float64,
        1,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{1,Index{Int64},Tuple{Index{Int64}}},
      },
      Tuple{Int64},
      Tuple{Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
      NDTensors.Tensor{
        Float64,
        1,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{1,Index{Int64},Tuple{Index{Int64}}},
      },
      Tuple{Int64},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
      Tuple{Int64,Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
      NTuple{4,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
      Tuple{Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
      Tuple{},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Float64,
        4,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}},
      },
      NTuple{4,Int64},
      NDTensors.Tensor{
        Float64,
        1,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{1,Index{Int64},Tuple{Index{Int64}}},
      },
      Tuple{Int64},
      Tuple{Int64,Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Float64,
        4,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}},
      },
      NTuple{4,Int64},
      NDTensors.Tensor{
        Float64,
        1,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{1,Index{Int64},Tuple{Index{Int64}}},
      },
      Tuple{Int64},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Float64,
        4,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}},
      },
      NTuple{4,Int64},
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
      Tuple{Int64,Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Float64,
        4,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}},
      },
      NTuple{4,Int64},
      NDTensors.Tensor{
        Float64,
        3,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      },
      Tuple{Int64,Int64,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Float64,
        4,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}},
      },
      NTuple{4,Int64},
      NDTensors.Tensor{
        Float64,
        4,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}},
      },
      NTuple{4,Int64},
      NTuple{4,Int64},
    },
  )
  precompile(
    Tuple{
      typeof(contract),
      NDTensors.Tensor{
        Float64,
        4,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}},
      },
      NTuple{4,Int64},
      NDTensors.Tensor{
        Float64,
        4,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}},
      },
      NTuple{4,Int64},
    },
  )
  precompile(
    Tuple{typeof(hassameinds),ITensor{3},Tuple{Index{Int64},Index{Int64},Index{Int64}}}
  )
  precompile(Tuple{typeof(hastags),TagSet,String})
  precompile(Tuple{typeof(itensor),Array{Float64,2},Index{Int64},Index{Int64}})
  precompile(
    Tuple{
      typeof(itensor),Array{Float64,3},Index{Int64},Index{Int64},Vararg{Index,N} where N
    },
  )
  precompile(Tuple{typeof(mul!),ITensor{3},ITensor{3},Float64})
  precompile(
    Tuple{
      typeof(mul!),
      ITensor{3},
      KrylovKit.OrthonormalBasis{ITensor{3}},
      SubArray{Float64,1,Array{Float64,2},Tuple{UnitRange{Int64},Int64},true},
    },
  )
  precompile(Tuple{typeof(prime),ITensor{3},ITensors.Not{TagSet}})
  precompile(
    Tuple{
      typeof(prime),
      IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      Int64,
      ITensors.Not{IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},
    },
  )
  precompile(Tuple{typeof(randomITensor),Index{Int64},Index{Int64},Index{Int64}})
  precompile(Tuple{typeof(randomMPS),Array{Index{Int64},1},Int64})
  precompile(Tuple{typeof(replaceind!),ITensor{2},Index{Int64},Index{Int64}})
  precompile(Tuple{typeof(replaceind!),ITensor{3},Index{Int64},Index{Int64}})
  precompile(
    Tuple{
      typeof(replaceind),
      IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      Index{Int64},
      Index{Int64},
    },
  )
  precompile(
    Tuple{
      typeof(replaceinds),
      IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      Tuple{Index{Int64},Index{Int64}},
      Tuple{Index{Int64},Index{Int64}},
    },
  )
  precompile(
    Tuple{
      typeof(replaceinds),
      IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},
      IndexSet{1,Index{Int64},Tuple{Index{Int64}}},
      IndexSet{1,Index{Int64},Tuple{Index{Int64}}},
    },
  )
  precompile(
    Tuple{
      typeof(setindex!),
      Dict{Array{ITensors.SiteOp,1},Array{Float64,2}},
      Array{Float64,2},
      Array{ITensors.SiteOp,1},
    },
  )
  precompile(
    Tuple{
      typeof(svd),
      NDTensors.Tensor{
        Complex{Float64},
        2,
        NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
    },
  )
  precompile(
    Tuple{
      typeof(svd),
      NDTensors.Tensor{
        Float64,
        2,
        NDTensors.Dense{Float64,Array{Float64,1}},
        IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},
      },
    },
  )
  return precompile(
    Tuple{
      typeof(zip),
      IndexSet{1,Index{Int64},Tuple{Index{Int64}}},
      Vararg{IndexSet{1,Index{Int64},Tuple{Index{Int64}}},N} where N,
    },
  )
end


text:
#!/usr/bin/env julia

module Jld2Mat

include("Params.jl")
include("Util.jl")

using .Util

printlog("Loading modules")
printlog("HDF5...")
using HDF5
printlog("JLD...")
using JLD
printlog("MAT...")
using MAT

function main()
    if length(ARGS) == 0
        println("Usage: jld2mat file1.jld [file2.jld ...]")
    end
    for file in ARGS
        if ! isfile(file)
            error("'$file' is not a file")
        end
    end

    for file in ARGS
        outfile = "$file.mat"
        if endswith(lowercase(file), ".jld")
            outfile = string(file[1:end-4], ".mat")
        end
        printlog("Opening $file")
        jldopen(file, "r") do f
            printlog("Extracting variables from $file")
            vars = read(f)
            printlog("Writing to $outfile")
            matwrite(outfile, vars)
        end
    end
end

printlog("Running main()")

main()

end # module

ext: jl


text:
module Learnable

include("../src/Jaynes.jl")
using .Jaynes

function learnable_kernel(x::Float64)
    q = learnable(:q)
    m = rand(:m, Normal(x + q[1], exp(q[2])))
    return m
end

function learnable_sim(q::Float64, len::Int)
    z = learnable(:z)
    x = markov(:x, learnable_kernel, len, z + q)
    return x
end

# Initial parameters.
len = 200
ps = parameters([(:x, :q) => [1.0, 1.0],
                 (:z, ) => 3.0])
sel = selection([(:x, i, :m) => 20.0 for i in 1 : len])

# Train.
ps = train(sel, ps, learnable_sim, 0.0, len;
           opt = ADAM(0.05, (0.9, 0.8)),
           iters = 3000)
display(ps)

end # module


text:
using Test
using Wflow
using Dates
using NCDatasets

tomlpath = joinpath(@__DIR__, "sbm_simple.toml")
Wflow.run(tomlpath)

# test whether restarted runs get the same results as continuous ones, i.e. state is captured
tomlpath = joinpath(@__DIR__, "sbm_config.toml")
config = Wflow.Config(tomlpath)

# January at once, cold start
config.starttime = DateTime("2000-01-01T00:00:00")
config.endtime = DateTime("2000-02-01T00:00:00")
config.model.reinit = true  # cold start
# note that this needs to be relative to the tomlpath
config.state.path_output =
    joinpath(dirname(tomlpath), "data/state-test/outstates-moselle-january.nc")
config.output.path =
    joinpath(dirname(tomlpath), "data/state-test/output-moselle-january.nc")
model = Wflow.initialize_sbm_model(config)
Wflow.run(model)

# first half of January, cold start
config.starttime = DateTime("2000-01-01T00:00:00")
config.endtime = DateTime("2000-01-15T00:00:00")
config.model.reinit = true  # cold start
config.state.path_output =
    joinpath(dirname(tomlpath), "data/state-test/outstates-moselle-january-1of2.nc")
config.output.path =
    joinpath(dirname(tomlpath), "data/state-test/output-moselle-january-1of2.nc")
model = Wflow.initialize_sbm_model(config)
Wflow.run(model)

# second half of January, warm start
config.starttime = DateTime("2000-01-16T00:00:00")
config.endtime = DateTime("2000-02-01T00:00:00")
config.model.reinit = false  # warm start
config.fews_run = false
config.state.path_input =
    joinpath(dirname(tomlpath), "data/state-test/outstates-moselle-january-1of2.nc")
config.state.path_output =
    joinpath(dirname(tomlpath), "data/state-test/outstates-moselle-january-2of2.nc")
config.output.path =
    joinpath(dirname(tomlpath), "data/state-test/output-moselle-january-2of2.nc")
model = Wflow.initialize_sbm_model(config)
Wflow.run(model)

# second half of January, warm start, fews_run set to true, and starttime set one day earlier
# to match endtime of part 1
config.starttime = DateTime("2000-01-15T00:00:00")
config.endtime = DateTime("2000-02-01T00:00:00")
config.model.reinit = false  # warm start
config.fews_run = true
config.state.path_input =
    joinpath(dirname(tomlpath), "data/state-test/outstates-moselle-january-1of2.nc")
config.state.path_output = joinpath(
    dirname(tomlpath),
    "data/state-test/outstates-moselle-january-2of2-fews_run.nc",
)
config.output.path =
    joinpath(dirname(tomlpath), "data/state-test/output-moselle-january-2of2-fews_run.nc")
model = Wflow.initialize_sbm_model(config)
Wflow.run(model)

# verify that there are minimal differences in the end state of the two runs
endstate_one_run_path =
    joinpath(dirname(tomlpath), "data/state-test/outstates-moselle-january.nc")
endstate_restart_path =
    joinpath(dirname(tomlpath), "data/state-test/outstates-moselle-january-2of2.nc")
endstate_one_run = NCDataset(endstate_one_run_path)
endstate_restart = NCDataset(endstate_restart_path)

varnames = setdiff(keys(endstate_restart), keys(endstate_restart.dim))
@testset "restart" begin
    @test length(varnames) > 10
    for varname in varnames
        a = endstate_one_run[varname][:]
        b = endstate_restart[varname][:]
        maxdiff = maximum(abs.(skipmissing(b - a)))
        @test maxdiff < 1e-9
    end
end

# the fews_run restart should match the other restart exactly
endstate_fewsrun_path = joinpath(
    dirname(tomlpath),
    "data/state-test/outstates-moselle-january-2of2-fews_run.nc",
)
endstate_restart_path =
    joinpath(dirname(tomlpath), "data/state-test/outstates-moselle-january-2of2.nc")
endstate_fewsrun = NCDataset(endstate_fewsrun_path)
endstate_restart = NCDataset(endstate_restart_path)

varnames = setdiff(keys(endstate_restart), keys(endstate_restart.dim))
@testset "fews_run" begin
    for varname in varnames
        a = endstate_fewsrun[varname][:]
        b = endstate_restart[varname][:]
        @test all(skipmissing(a) .== skipmissing(b))
    end
end


text:
## Foil type
#==========================================================================================#

"""
    Foil(coords)

Airfoil structure consisting of foil coordinates as an array of points. Should be in Selig format for compatibility with other AeroMDAO tools.
"""
struct Foil{T <: Real}
    coords :: Vector{SVector{2,T}}
    name :: String
    Foil(coords :: Vector{SVector{2,T}}, name = "Unnamed") where T <: Real = new{T}(coords, name)
end

"""
    scale_foil(foil :: Foil, chord)

Scale the coordinates of a Foil, usually to some chord length.
"""
scale_foil(foil :: Foil, chord) = chord * foil.coords

"""
    cosine_foil(foil :: Foil, num :: Integer)

Return a Foil with cosine spacing for a given number of points.
"""
cosine_foil(foil :: Foil, num :: Integer) = cosine_foil(foil.coords, num)

"""
    camber_thickness(foil :: Foil, num :: Integer)

Compute the camber-thickness distribution of a Foil with cosine spacing..
"""
camber_thickness(foil :: Foil, num :: Integer) = foil_camthick(cosine_foil(foil.coords), num + 1)

function max_thickness_to_chord_ratio_location(coords)
    xs, thiccs = getindex.(coords, 1), getindex.(coords, 3)
    max_thick_arg = argmax(thiccs)
    xs[max_thick_arg], thiccs[max_thick_arg]
end

## Foil processing
#==========================================================================================#

"""
    read_foil(path :: String; header = true)

Read a '.dat' file consisting of 2D coordinates, for an airfoil as an array of `SVector`s, with an optional argument to skip the header.
"""
function read_foil(path :: String; header = true)
    coords = readdlm(path, skipstart = header ? 1 : 0)
    @views SVector.(coords[:,1], coords[:,2])
end

function split_foil(coords)
    for (i, ((xp, yp), (x, y), (xn, yn))) ∈ (enumerate ∘ adj3)(coords)
        if x < xp && x < xn
            if slope(x, y, xp, yp) >= slope(x, y, xn, yn)
                return splitat(i, coords)
            else
                return splitat(i, coords[end:-1:1,:])
            end
        end
    end
    (coords, [])
end

function paneller(foil :: Foil, num_panels :: Integer)
    coords = cosine_foil(foil.coords, Int(ceil(num_panels / 2)))
    @views Panel2D.(coords[2:end], coords[1:end-1])[end:-1:1]
end

"""
Discretises a foil profile into panels by projecting the x-coordinates of a circle onto the geometry.
"""
function cosine_foil(coords, n :: Integer = 40)
    upper, lower = split_foil(coords)
    n_upper = [ upper    ;
               [lower[1]]] # Append leading edge point from lower to upper
    upper_cos, lower_cos = cosine_interp(n_upper[end:-1:1], n), cosine_interp(lower, n)

    @views [ upper_cos[end:-1:2] ;
             lower_cos           ]
end

## Class shape transformation method
#==========================================================================================#

# Basic shape function
function shape_function(x, basis_func, coeffs, coeff_LE = 0)
    n     = length(coeffs)
    terms = basis_func.(x, n - 1, 0:n-1)
    dot(coeffs, terms) + coeff_LE * (x^0.5) * (1 - x)^(n - 0.5)
end

# Computing coordinates
cst_coords(class_func, basis_func, x, alphas, dz, coeff_LE, args...) = class_func(x) * shape_function(x, basis_func, alphas, coeff_LE) + x * dz

## Bernstein basis
#==========================================================================================#

bernstein_class(x, N1, N2) = x^N1 * (1 - x)^N2
bernstein_basis(x, n, k)   = binomial(n, k) * bernstein_class(x, k, n - k)

"""
    kulfan_CST(alpha_u, alpha_l,
               (dz_u, dz_l) :: NTuple{2, Real},
               coeff_LE = 0,
               n :: Integer = 40)

Define a cosine-spaced airfoil with ``2n`` points using the Class Shape Transformation method on a Bernstein polynomial basis with arrays of coefficients ``(\\alpha_u,~ \\alpha_l)`` for the upper and lower surfaces, trailing-edge spacing values ``(\\Delta z_u,~ \\Delta z_l)``, and support for leading edge modifications.
"""
function kulfan_CST(alpha_u, alpha_l, (dz_u, dz_l), coeff_LE = 0, n :: Integer = 40, N1 = 0.5, N2 = 1.)
    # Cosine spacing for airfoil of unit chord length
    xs = cosine_spacing(0.5, 1, n)

    # λ-function for Bernstein polynomials
    bernie(x, alphas, dz) = cst_coords(y -> bernstein_class(y, N1, N2), bernstein_basis, x, alphas, dz, coeff_LE)

    # Upper and lower surface generation
    upper_surf = [ bernie(x, alpha_u, dz_u) for x ∈ xs ]
    lower_surf = [ bernie(x, alpha_l, dz_l) for x ∈ xs ]

    # Counter-clockwise ordering
    @views [ SVector.(xs, upper_surf)[end:-1:2] ;
             SVector.(xs, lower_surf)           ]
end

function camber_CST(alpha_cam, alpha_thicc, (dz_cam, dz_thicc), coeff_LE = 0, n :: Integer = 40, N1 = 0.5, N2 = 1.)
    # Cosine spacing for airfoil of unit chord length
    xs = cosine_spacing(0.5, 1, n)

    # λ-function for Bernstein polynomials
    bernie(x, alphas, dz) = cst_coords(y -> bernstein_class(y, N1, N2), bernstein_basis, x, alphas, dz, coeff_LE)

    # Upper and lower surface generation
    cam   = [ bernie(x, alpha_cam, dz_cam) for x ∈ xs ]
    thicc = [ bernie(x, alpha_thicc, dz_thicc) for x ∈ xs ]

    camthick_foil(xs, cam, thicc)
end

function coords_to_CST(coords, num_dvs)
    xs 		 = first.(coords)
    S_matrix = reduce(hcat, @. bernstein_class(xs, 0.5, 1.0) * bernstein_basis(xs, num_dvs - 1, i) for i in 0:num_dvs - 1)

    alphas 	 = S_matrix \ last.(coords)

    return alphas
end

function camthick_to_CST(coords, num_dvs)
    xs, camber, thickness = (columns ∘ foil_camthick)(coords)

    alpha_cam  	= coords_to_CST(SVector.(xs, camber), num_dvs)
    alpha_thick = coords_to_CST(SVector.(xs, thickness), num_dvs)

    alpha_cam, alpha_thick
end

## Camber-thickness representation
#==========================================================================================#

"""
Converts an airfoil to its camber-thickness representation in cosine spacing.
"""
function foil_camthick(coords, num :: Integer = 40)
    upper, lower = split_foil(cosine_foil(coords, num))

    xs, y_LE  	 = first.(lower), first(lower)[2]   # Getting abscissa and leading edge ordinate
    y_upper, y_lower = last.(upper[end:-1:1]), last.(lower[2:end]) # Excluding leading edge point

    camber    	 = [ y_LE; (y_upper + y_lower) / 2 ]
    thickness 	 = [ 0; y_upper - y_lower ]

    SVector.(xs, camber, thickness)
end

"""
    camthick_foil(xs, camber, thickness)

Converts the camber-thickness representation to coordinates given the ``x``-locations and their corresponding camber and thickness values.
"""
function camthick_foil(xs, camber, thickness)
    coords = [ [xs camber + thickness / 2][end:-1:2,:];
                xs camber - thickness / 2 			  ]
    @views SVector.(coords[:,1], coords[:,2])
end

camthick_foil(coords) = @views camthick_foil(getindex.(coords, 1), getindex.(coords, 2), getindex.(coords, 3))

"""
    camber_coordinates(coords :: Array{2, <: Real})

Provides the camber coordinates on the ``x``-``z`` plane at ``y = 0``.
"""
camber_coordinates(coords) = @views SVector.(getindex.(coords, 1), 0, getindex.(coords, 2))

"""
    thickness_coordinates(coords :: Array{2, <: Real})

Provides the thickness coordinates on the ``x``-``z`` plane at ``y = 0``.
"""
thickness_coordinates(coords) = @views SVector.(getindex.(coords, 1), 0, getindex.(coords, 3))

## NACA Parametrisation
#==========================================================================================#

# NACA 4-digit parameter functions
naca4_thickness(t_by_c, xc, sharp_trailing_edge :: Bool) = 5 * t_by_c * (0.2969 * √xc - 0.1260 * xc - 0.3516 * xc^2 + 0.2843 * xc^3 - (ifelse(sharp_trailing_edge, 0.1036, 0.1015) * xc^4))
naca4_camberline(pos, cam, xc) = ifelse(xc < pos, (cam / pos^2) * xc * (2 * pos - xc), cam / (1 - pos)^2 * ( (1 - 2 * pos) + 2 * pos * xc - xc^2) )
naca4_gradient(pos, cam, xc) = atan(2 * cam / (ifelse(xc < pos, pos^2, (1 - pos)^2)) * (pos - xc))

"""
    naca4(digits :: NTuple{4, <: Real}, n :: Integer; sharp_trailing_edge :: Bool)

Defines a NACA 4-digit series profile.
"""
function naca4(digits :: NTuple{4, <: Real}, n :: Integer = 40; sharp_trailing_edge = true)
    # Camber
    cam = digits[1] / 100
    # Position
    pos = digits[2] / 10
    # Thickness-to-chord ratio
    t_by_c = (10 * digits[3] + digits[4]) / 100

    # Cosine spacing
    xs = cosine_spacing(0.5, 1.0, n)

    # Thickness distribution
    thickness = naca4_thickness.(Ref(t_by_c), xs, Ref(sharp_trailing_edge))

    if pos == 0 || cam == 0
        x_upper = xs
        y_upper = thickness
        x_lower = xs
        y_lower = -thickness
    else
        # Compute camberline
        camber 	= naca4_camberline.(Ref(pos), Ref(cam), xs)
        # Compute gradients
        grads 	= naca4_gradient.(Ref(pos), Ref(cam), xs)
        # Upper surface
        x_upper = @. xs - thickness * sin(grads)
        y_upper = @. camber + thickness * cos(grads)
        # Lower surface
        x_lower = @. xs + thickness * sin(grads)
        y_lower = @. camber - thickness * cos(grads)
    end
    coords = [ [x_upper y_upper][end:-1:2,:];
                x_lower y_lower             ]

    @views SVector.(coords[:,1], coords[:,2])
end

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
export VAE_Encoder, kl_divergence, bernoulli_expected_loglik, normal_expected_loglik
export mean, hidden, log_std

##################
# VAE_Encoder
##################
struct VAE_Encoder{T, MT<:AbstractMatrix{T}, VT<:AbstractVector{T}}
	Wh::MT
	bh::VT
	Wμ::MT
	bμ::VT
	Wlogσ::MT
	blogσ::VT
	f::Function
end
Flux.@functor VAE_Encoder
function VAE_Encoder(in::Integer, h_dim::Integer, out::Integer,
								f=identity; initW=glorot_uniform, initb=zeros)
	Wh = initW(h_dim, in)
	bh = initb(Float32, h_dim)
	Wμ = initW(out, h_dim)
	bμ = initb(Float32, out)
	Wlogσ = initW(out, h_dim)
	blogσ = initb(Float32, out)
	VAE_Encoder(Wh, bh, Wμ, bμ, Wlogσ, blogσ, f)
end

Base.eltype(::Type{<:VAE_Encoder{T}}) where {T} = T
hidden(d::VAE_Encoder, x) = d.f.(d.Wh*x .+ d.bh)
mean(d::VAE_Encoder, h) = d.Wμ*h .+ d.bμ
log_std(d::VAE_Encoder, h) = d.Wlogσ*h .+ d.blogσ

function (d::VAE_Encoder)(rng::AbstractRNG, x)
	h = hidden(d, x)
	μ = mean(d, h)
	logσ = log_std(d, h)
	z = μ .+ exp.(logσ).*randn(rng, eltype(d), size(logσ))
	z
end
(d::VAE_Encoder)(x) = (d::VAE_Encoder)(GLOBAL_RNG, x)

function Base.show(io::IO, en::VAE_Encoder)
	print(io, "VAE(", size(en.Wh, 2), ", ", size(en.Wh, 1), ", ", size(en.Wμ, 1))
	en.f == identity || print(io, ", ", en.f)
	print(io, ")")
end

#############
# KL
#############
function kl_divergence(d::VAE_Encoder, x)
	h = hidden(d, x)
	μ = mean(d, h)
	logσ = log_std(d, h)
	kl = sum(logσ .+ 0.5*μ.^2 .+ 0.5*exp.(2.0*logσ) .- 1)
	kl
end

#####################
# expected log lik
#####################
function bernoulli_expected_loglik(fp, en::VAE_Encoder, x, sample_z=1)
	loglik = 0.0f0
	for i in 1:sample_z
		z = en(x)
		p = fp(z)
		loglik += sum(logpdf.(Bernoulli.(p), x))
	end
	loglik / sample_z
end

function normal_expected_loglik(fμ, flogσ, en::VAE_Encoder, x, sample_z=1)
	loglik = 0.0f0
	for i in 1:sample_z
		z = en(x)
		μ = fμ(z)
		logσ = flogσ(z)
		σ = exp.(logσ)
		loglik += sum(logpdf.(Normal.(μ, σ), x))
	end
	loglik / sample_z
end

max_stars_repo_name: HamletWantToCode/Generative

max_issues_repo_name: HamletWantToCode/Generative

max_forks_repo_name: HamletWantToCode/Generative


text:
using PyPlot

include("./topoprovision.jl")
include("./topodraw.jl")

## Increment Counter and Get Next Dir
function nextRun(dir, singleplot)
    # dir = joinpath(@__DIR__, "./$dir")
    if isfile("$(dir)/counter.txt") == false
        mkdir("$(dir)")
        writedlm("$(dir)/counter.txt", "1")
    end
    if singleplot
        c = Int(readdlm("$(dir)/counter.txt")[1])
        writedlm("$(dir)/counter.txt", c + 1)
        newdir = "$(dir)/$(c)"
        mkdir(newdir)
    else
        if Int(readdlm("$(dir)/counter.txt")[1]) == 1  # handling case when simulation last longer than a day
            c = Int(readdlm("$(dir)/counter.txt")[1])
            writedlm("$(dir)/counter.txt", c + 1)
            newdir = "$(dir)/$(c)"
            if isdir(newdir) == false
                mkdir(newdir)
            end
        else
            c = Int(readdlm("$(dir)/counter.txt")[1])-1
        end
        newdir = "$(dir)/$(c)"
    end

    return newdir
end

## Compute Weibull Probabilities
function weibullProbs(num; shape=.8, scale=.0001)
    w = Distributions.Weibull(shape, scale)
    probs = []
    for i in 1:num
        push!(probs, rand(w))
    end
    return probs
end

## get optical/IP topo
function ReadCrossLayerTopo(dir, topology, topology_index, verbose, expanded_spectrum; weibull_failure=false, IPfromFile=false, tofile=false)
    topodir =  "./data/topology/$(topology)"
    input_ip_nodes = readdlm("$(topodir)/IP_nodes.txt", header=true)[1]
    input_optical_nodes = readdlm("$(topodir)/optical_nodes.txt", header=true)[1]
    input_topology_o = readdlm("$(topodir)/optical_topo.txt", header=true)[1]
    ignore = ()
    fromNodesOptical = input_topology_o[:,1]
    toNodesOptical = input_topology_o[:,2]
    lengthOptical = input_topology_o[:,3]
    if isdir("$(dir)/$(topology)") == false
        mkdir("$(dir)/$(topology)")
    end

    if weibull_failure
        failureprobOptical = weibullProbs(length(lengthOptical), shape=.8, scale=.02)  # derive fiber cut probability from weibull
        ## draw initial weibull distribution CDF
        PyPlot.clf()
        sorted_failureprobOptical = sort(failureprobOptical)
        cdf = []
        for i in 1:length(sorted_failureprobOptical)
            push!(cdf, i/length(sorted_failureprobOptical))
        end
        PyPlot.plot(sorted_failureprobOptical, cdf, marker="P", linewidth=1, label="shape=.8, scale=.02")
        figname = "$(dir)/$(topology)/00_weibull_distribution.png"
        PyPlot.xlabel("Fiber Failure scenario probability")
        PyPlot.ylabel("CDF")
        PyPlot.legend(loc="best")
        PyPlot.xscale("log")
        PyPlot.savefig(figname)
    else
        failureprobOptical = input_topology_o[:,4]  # if read fiber cut probability from file
        ## draw initial weibull distribution CDF
        PyPlot.clf()
        sorted_failureprobOptical = sort(failureprobOptical)
        cdf = []
        for i in 1:length(sorted_failureprobOptical)
            push!(cdf, i/length(sorted_failureprobOptical))
        end
        PyPlot.plot(sorted_failureprobOptical, cdf, marker="P", linewidth=1, label="failure distribution from file")
        # PyPlot.plot(sorted_failureprobOptical.* 20, cdf, marker="P", linewidth=1, label="shifted")
        figname = "$(dir)/$(topology)/00_filefailure_distribution.png"
        PyPlot.xlabel("Fiber Failure scenario probability")
        PyPlot.ylabel("CDF")
        PyPlot.legend(loc="best")
        PyPlot.xscale("log")
        PyPlot.savefig(figname)
    end

    fiberlinks = []
    bidirect_links = []
    fiberlinkslength = []
    fiberlinksSpectrum = []
    fiberlinksFailure = []
    bidirect_fiberlinksFailure = []

    for i in 1:size(fromNodesOptical, 1)
        if (!(fromNodesOptical[i] in ignore) && !(toNodesOptical[i] in ignore))
            push!(fiberlinks, (Int(fromNodesOptical[i]), Int(toNodesOptical[i])))
            push!(fiberlinkslength, lengthOptical[i])
            push!(fiberlinksSpectrum, [])
            push!(fiberlinksFailure, failureprobOptical[i])
            if !in((Int(fromNodesOptical[i]), Int(toNodesOptical[i])), bidirect_links) && !in((Int(toNodesOptical[i]), Int(fromNodesOptical[i])), bidirect_links)
                push!(bidirect_links, (Int(fromNodesOptical[i]), Int(toNodesOptical[i])))
                push!(bidirect_fiberlinksFailure, failureprobOptical[i])
            end
        end
    end

    parsed_wavelength = []
    parsed_optical_route = []
    links_length = []
    if IPfromFile == true
        input_topology_ip = readdlm("$(topodir)/IP_topo_$(topology_index)/IP_topo_$(topology_index).txt", header=true)[1]
        if verbose println("Read IP layer topo from file: $(topodir)/IP_topo_$(topology_index)/IP_topo_$(topology_index).txt") end
        ignore = ()
        optical_route = input_topology_ip[:,5]
        wavelength = input_topology_ip[:,6]
        for x in 1:length(optical_route)
            current_wavelength = [parse(Int, ss) for ss in split(wavelength[x][2:end-1], ",")]  # example: [1,2,3,4,5,6,7]
            push!(parsed_wavelength, current_wavelength)
            current_optical_route = [parse(Int, ss) for ss in split(optical_route[x][2:end-1], ",")]
            push!(parsed_optical_route, current_optical_route)
            current_length = 0
            for j in 1:length(current_optical_route)
                current_length += lengthOptical[j]
                for t in 1:length(current_wavelength)
                    if current_wavelength[t] in fiberlinksSpectrum[current_optical_route[j]]
                        continue
                    else
                        push!(fiberlinksSpectrum[current_optical_route[j]],current_wavelength[t])
                    end
                end
            end
            push!(links_length, current_length)
        end
    else
        input_topology_ip = ProvisionIPTopology(topology, tofile=tofile)  # generate IP provisioning only no output to file
        optical_route = input_topology_ip[:,5]
        wavelength = input_topology_ip[:,6]
        for x in 1:length(optical_route)
            current_wavelength = [parse(Int, ss) for ss in split(wavelength[x][5:end-1], ", ")]  # example: Any[1, 2, 3, 4, 5, 6, 7]
            push!(parsed_wavelength, current_wavelength)
            current_optical_route = [parse(Int, ss) for ss in split(optical_route[x][5:end-1], ", ")]
            push!(parsed_optical_route, current_optical_route)
            current_length = 0
            for j in 1:length(current_optical_route)
                current_length += lengthOptical[j]
                for t in 1:length(current_wavelength)
                    if current_wavelength[t] in fiberlinksSpectrum[current_optical_route[j]]
                        continue
                    else
                        push!(fiberlinksSpectrum[current_optical_route[j]], current_wavelength[t])
                    end
                end
            end
            push!(links_length, current_length)
        end
    end

    fromNodes = input_topology_ip[:,1]
    toNodes = input_topology_ip[:,2]
    marker = input_topology_ip[:,3]
    capacity = input_topology_ip[:,4] .* 100  # we assume each wavelength is 100 Gbps
    failure_prob = input_topology_ip[:,7]

    links = []
    links_fiberroute = []
    links_wavelength = []
    links_hops = []
    # println("fromNodes ", fromNodes)

    for i in 1:size(fromNodes, 1)
        if (!(fromNodes[i] in ignore) && !(toNodes[i] in ignore))
            push!(links, (Int(fromNodes[i]), Int(toNodes[i]), Int(marker[i])))
            push!(links_fiberroute, parsed_optical_route[i])
            push!(links_wavelength, parsed_wavelength[i])
            push!(links_hops, length(parsed_optical_route[i]))
        end
    end

    # IP links are represented in uni-directional form
    IPTopo = Dict()
    IPTopo["nodes"] = input_ip_nodes
    IPTopo["links"] = links
    IPTopo["capacity"] = capacity
    IPTopo["link_probs"] = failure_prob
    IPTopo["link_fiberroute"] = links_fiberroute
    IPTopo["link_length"] = links_length
    IPTopo["link_wavelength"] = links_wavelength

    # optical fiber links are represented in uni-directional form
    OpticalTopo = Dict()
    OpticalTopo["nodes"] = input_optical_nodes
    OpticalTopo["links"] = fiberlinks  # unidirectional
    OpticalTopo["bidirect_links"] = bidirect_links
    OpticalTopo["fiber_length"] = fiberlinkslength
    OpticalTopo["fiber_spectrum"] = fiberlinksSpectrum
    OpticalTopo["fiber_probs"] = fiberlinksFailure
    OpticalTopo["bidirect_fiber_probs"] = bidirect_fiberlinksFailure
    OpticalTopo["capacity"] = [96+expanded_spectrum-length(x) for x in fiberlinksSpectrum]

    capacitycode = ones(length(fiberlinks), 96+expanded_spectrum)
    for fiber in 1:length(fiberlinksSpectrum)
        for slot in fiberlinksSpectrum[fiber]
            capacitycode[fiber, slot] = 0
        end
    end
    OpticalTopo["capacityCode"] = capacitycode  # 0 means occupied, 1 means available

    drawGraph(topology, IPTopo, topology_index, OpticalTopo, dir)

    return IPTopo, OpticalTopo
end

## different TM may have different sum demand, this function find the proper demand rescale if we want to consider a uniform constant demand time series
function FindDemandScale(topology, nodes_num, AllTraffic)
    demand_upscale = 1
    demand_downscale = 1
    sum_demands = []
    scaling_factor = []
    for traffic_num in 1:length(AllTraffic)
        initial_demand, flows = readDemand("$(topology)/demand", nodes_num, AllTraffic[traffic_num], demand_upscale, demand_downscale, false)  # read initial demand
        # println(sum(initial_demand))
        push!(sum_demands, sum(initial_demand))
    end

    min_demand = minimum(sum_demands)
    for i in 1:length(sum_demands)
        push!(scaling_factor, sum_demands[i]/min_demand)
    end

    return scaling_factor
end

## Read Demand
function readDemand(filename, num_nodes, num_demand, scale, downscale, rescale; matrix=true, sigfigs=1, zeroindex=false)
    filename = "./data/topology/$(filename)"
    input_demand = matrix ? ParseMatrix("$(filename).txt", num_nodes, num_demand) : IgnoreCycles(readdlm("$(filename)/$(num_demand).txt", header=true)[1], zeroindex=zeroindex)
    # println("input_demand ", input_demand)
    fromNodes = input_demand[:,1]
    toNodes = input_demand[:,2]
    flows = map(tup -> (Int(tup[2]), Int(toNodes[tup[1]])), enumerate(fromNodes))
    demand = convert(Array{Float64}, input_demand[:,3] ./ downscale .* scale)

    if rescale
        rescale_set = readdlm("$(filename)_rescale.txt", header=false)
        # println("rescale_set ", rescale_set)
        demand = demand ./ rescale_set[num_demand]
        println("Rescale demand $(num_demand) at $(rescale_set[num_demand]) -> $(sum(demand))")
    end

    return demand, flows
end

## IgnoreCycles
function IgnoreCycles(demand; zeroindex=false)
    z = [0 0 0;]
    for row in 1:size(demand,1)
        if demand[row, 1] != demand[row, 2]
            z = vcat(z, [demand[row,1] + zeroindex demand[row, 2] + zeroindex demand[row,3];])
            # z = vcat(z, transpose(demand[row,:]))
        end
    end
    return z[2:end,:]

end

## Parse Matrix
function ParseMatrix(filename, num_nodes, num_demand)
    ignore = ()
    start_range = 0
    end_range = Inf
    x = readdlm(filename)[num_demand,:]
    # println(x)
    m = zeros(length(x)-Int(sqrt(length(x))), 3)
    # m = zeros(length(x)-round(Int, sqrt(length(x))), 3)
    fromNode = 0
    count = 1
    for i in 0:(num_nodes^2-1)
        toNode = i%num_nodes + 1
        if toNode == 1
            fromNode += 1
        end
        if (fromNode != toNode && i >= start_range && i < end_range && !(fromNode in ignore) && !(toNode in ignore))
            m[count,1] = fromNode
            m[count,2] = toNode
            m[count,3] = x[i+1]
            count += 1
        end
    end
    ret = m
    for row in size(m,1):-1:1
        if (m[row,3] == 0)
            ret = ret[setdiff(1:end, row), :]
        end
    end
    return ret
end

## Tunnels routing
function getTunnels(IPtopo, flows, k, verbose, scenarios; edge_disjoint=1)  # 1 means k shortest path, 2 means IPedge disjoint, 3 means fiber disjoint
    nodes = IPtopo["nodes"]
    edges = IPtopo["links"]
    capacities = IPtopo["capacity"]
    fiberroute = IPtopo["link_fiberroute"]
    num_edges = length(edges)
    # print("num_edges:", num_edges)
    # print(edges)
    num_nodes = length(nodes)
    num_flows = length(flows)
    paral_edges = []
    # count number of paralle links for adding dummy links/nodes
    dummy_nodes = 0
    for l in 1:num_edges
        if edges[l][3] > 1  # parallel IP links
            dummy_nodes += 1
            push!(paral_edges, edges[l])
        end
    end

    graph = LightGraphs.SimpleDiGraph(num_nodes + dummy_nodes)  # considering parallel links
    distances = Inf*ones(num_nodes+dummy_nodes, num_nodes+dummy_nodes)

    # println("paral_edges: ", length(paral_edges), paral_edges)

    for i in 1:num_edges
        if edges[i][3] == 1
            LightGraphs.add_edge!(graph, edges[i][1], edges[i][2])
            if edge_disjoint == 5  # weight KSP
                distances[edges[i][1], edges[i][2]] = 1 / capacities[i]  # weight edge, more capacity means less distance
                distances[edges[i][2], edges[i][1]] = 1 / capacities[i]  # weight edge, more capacity means less distance
            else
                distances[edges[i][1], edges[i][2]] = 1
                distances[edges[i][2], edges[i][1]] = 1
            end
        else
            e = findfirst(x -> x == (edges[i][1], edges[i][2], edges[i][3]), paral_edges)
            # println("parallel link: ", e, paral_edges[e])
            LightGraphs.add_edge!(graph, edges[i][1], num_nodes+e)
            LightGraphs.add_edge!(graph, num_nodes+e, edges[i][2])
            distances[edges[i][1], num_nodes+e] = 0.5
            distances[num_nodes+e, edges[i][2]] = 0.5
            # distances[edges[i][2], num_nodes+e] = 0.5
            # distances[num_nodes+e, edges[i][1]] = 0.5
        end
    end

    T = []
    Tf = []
    ti = 1  # global tunnel index
    max_k = 1

    # KSP routing
    if edge_disjoint == 1
        printstyled("KSP tunnel routing\n", color=:blue)
        for f in 1:num_flows
            tf = []
            curr_k = 0

            state = LightGraphs.yen_k_shortest_paths(graph, flows[f][1], flows[f][2], distances, k)
            paths = state.paths
            edges_used = []
            fiber_used = []
            for i in 1:k
                t = []
                traversed_fiber = []
                path = i <= size(paths,1) ? paths[i] : []
                jump_flag = 0
                for n in 2:size(path,1)
                    if path[n] > num_nodes   # use parallel links
                        e = findfirst(x -> x == (path[n-1], path[n+1], paral_edges[path[n]-num_nodes][3]), edges)
                        jump_flag = 1
                        push!(t, e)
                    elseif jump_flag == 0
                        e = findfirst(x -> x == (path[n-1], path[n], 1), edges)
                        push!(t, e)
                        for b in 1:length(fiberroute[e])
                            if in(fiberroute[e][b], traversed_fiber)
                                continue
                            else
                                push!(traversed_fiber, fiberroute[e][b])
                            end
                        end
                        f = edges[e]
                    elseif jump_flag == 1
                        jump_flag = 0
                    end
                end

                if length(t) == 0 break end

                push!(T, t)
                push!(tf, ti)
                ti += 1
                curr_k += 1

            end
            push!(Tf, tf)
            if verbose println("K-shortest tunnel routing: $(paths)") end
        end

    # IPedge disjoint routing
    elseif edge_disjoint == 2
        printstyled("IP disjoint tunnel routing\n", color=:blue)
        for f in 1:num_flows
            tf = []
            disjoint_graph = deepcopy(graph)
            for i in 1:k
                state = LightGraphs.yen_k_shortest_paths(disjoint_graph, flows[f][1], flows[f][2], distances, 1)
                paths = state.paths
                # println("paths ", paths)
                t = []
                path = 0 < size(paths,1) ? paths[1] : []  # because we use KSP k=1
                # println("path ", path)
                for n in 2:size(path,1)
                    e = findfirst(x -> x == (path[n-1], path[n], 1), edges)
                    push!(t, e)
                    remove_status = LightGraphs.rem_edge!(disjoint_graph, path[n-1], path[n])
                    # println("remove edge $(path[n-1])-$(path[n]), $(remove_status)")
                    distances[path[n-1], path[n]] = Inf
                end

                if length(t) == 0 break end

                push!(T, t)
                push!(tf, ti)
                ti += 1
            end
            push!(Tf, tf)
        end

    # fiber disjoint routing
    elseif edge_disjoint == 3
        printstyled("fiber disjoint tunnel routing\n", color=:blue)
        for f in 1:num_flows
            tf = []
            disjoint_graph = deepcopy(graph)
            for i in 1:k
                state = LightGraphs.yen_k_shortest_paths(disjoint_graph, flows[f][1], flows[f][2], distances, 1)
                paths = state.paths
                # println("paths ", paths)
                t = []
                path = 0 < size(paths,1) ? paths[1] : []  # because we use KSP k=1
                # println("path ", path)
                for n in 2:size(path,1)
                    e = findfirst(x -> x == (path[n-1], path[n], 1), edges)
                    push!(t, e)
                    # remove all IP edges on related b fiber
                    for b in 1:length(fiberroute[e])
                        for ee in 1:num_edges
                            if fiberroute[e][b] in fiberroute[ee]
                                remove_status = LightGraphs.rem_edge!(disjoint_graph, edges[ee][1], edges[ee][2])
                                # println("remove edge $(edges[ee][1])-$(edges[ee][2]), $(remove_status)")
                                distances[edges[ee][1], edges[ee][2]] = Inf
                            end
                        end
                    end
                    # the reverse path fiber should also be removed because fiber are bidirectional
                    re = findfirst(x -> x == (path[n], path[n-1], 1), edges)
                    for b in 1:length(fiberroute[re])
                        for ee in 1:num_edges
                            if fiberroute[re][b] in fiberroute[ee]
                                remove_status = LightGraphs.rem_edge!(disjoint_graph, edges[ee][1], edges[ee][2])
                                # println("remove edge $(edges[ee][1])-$(edges[ee][2]), $(remove_status)")
                                distances[edges[ee][1], edges[ee][2]] = Inf
                            end
                        end
                    end
                end

                if length(t) == 0 break end

                push!(T, t)
                push!(tf, ti)
                ti += 1
            end
            push!(Tf, tf)
        end

    # failure aware routing, in real life tunnels should be properly routed such that failure will disrupt all tunnels of a flow
    elseif edge_disjoint == 4 || edge_disjoint == 5
        printstyled("failure aware tunnel routing\n", color=:blue)
        nscenarios = length(scenarios)
        for f in 1:num_flows
            accumulate_k = 0
            tf = []

            # failure scenario aware tunnel routing
            for s in 1:nscenarios
                disjoint_graph = deepcopy(graph)
                for q in 1:length(scenarios[s])
                    if scenarios[s][q] == 0
                        remove_status = LightGraphs.rem_edge!(disjoint_graph, edges[q][1], edges[q][2])
                        distances[edges[q][1], edges[q][2]] = Inf
                    end
                end
                state = LightGraphs.yen_k_shortest_paths(disjoint_graph, flows[f][1], flows[f][2], distances, 1)
                paths = state.paths
                t = []
                path = 0 < size(paths,1) ? paths[1] : []  # because we use KSP k=1

                for n in 2:size(path,1)
                    e = findfirst(x -> x == (path[n-1], path[n], 1), edges)
                    push!(t, e)
                end

                if length(t) == 0 break end

                if !in(t, T) && accumulate_k < k
                    push!(T, t)
                    push!(tf, ti)
                    ti += 1
                    accumulate_k += 1
                end
            end

            # if tunnel number is smaller than k, then use KSP to fill it
            if accumulate_k < k
                state = LightGraphs.yen_k_shortest_paths(graph, flows[f][1], flows[f][2], distances, k)
                paths = state.paths

                for i in 1:k
                    t = []
                    traversed_fiber = []
                    path = i <= size(paths,1) ? paths[i] : []
                    jump_flag = 0
                    for n in 2:size(path,1)
                        if path[n] > num_nodes   # use parallel links
                            e = findfirst(x -> x == (path[n-1], path[n+1], paral_edges[path[n]-num_nodes][3]), edges)
                            jump_flag = 1
                            push!(t, e)
                        elseif jump_flag == 0
                            e = findfirst(x -> x == (path[n-1], path[n], 1), edges)
                            push!(t, e)
                            for b in 1:length(fiberroute[e])
                                if in(fiberroute[e][b], traversed_fiber)
                                    continue
                                else
                                    push!(traversed_fiber, fiberroute[e][b])
                                end
                            end
                            f = edges[e]
                        elseif jump_flag == 1
                            jump_flag = 0
                        end
                    end

                    if length(t) == 0 break end

                    if !in(t, T) && accumulate_k < k
                        push!(T, t)
                        push!(tf, ti)
                        ti += 1
                        accumulate_k += 1
                    end
                end
            end
            push!(Tf, tf)
        end
    end
    return T, Tf, max_k, graph
end

## parse tunnels from input files (format 1)
function parseTunnels(allocs_tunnels_paths, links)
    println("Reading tunnels from: $(allocs_tunnels_paths)")
    allocs_tunnels = readdlm(allocs_tunnels_paths)

    num_lines = size(allocs_tunnels,1)
    flows = []
    T = [[]]
    Tf = []
    tf = []
    max_t_paths = 0
    num_flows = 0

    for i in 1:num_lines
        # Reading flows
        if allocs_tunnels[i, 1] == "Flow:"
            flow_str = allocs_tunnels[i, 2][2:end-1]
            f1, f2 = split(flow_str, ",")
            f1, f2 = parse(Int64, f1) + 1, parse(Int64, f2) + 1
            flow = (f1, f2)

            push!(flows, flow)
            if num_flows > 0
                push!(Tf, tf)
            end

            num_flows = num_flows + 1
            tf = []

        # Reading tunnels and allocations
        elseif allocs_tunnels[i, 2] == ""
            _hops = split(allocs_tunnels[i, 1], ";")[1:end-1]
            t = []
            for _hop in _hops
                _hop_str = _hop[2:end-1]
                _h1, _h2 = split(_hop_str, ",")
                _h1, _h2 = parse(Int64, _h1) + 1, parse(Int64, _h2) + 1
                _hop_int = (_h1, _h2)

                if _h1 != 0 && _h2 != 0
                    e = findfirst(x -> (x[1], x[2]) == _hop_int, links)
                    if e !== nothing
                        push!(t, e)
                    end
                end
            end

            if !in(t, T)
                push!(T, t)
            end

            _t_idx = findfirst(x -> x == t, T)
            push!(tf, _t_idx)
            max_t_paths = max(max_t_paths, size(tf, 1))
        end
    end

    push!(Tf, tf)

    Tf_adjusted = []
    t_empty = findfirst(x -> x == [], T)

    for f in 1:size(flows, 1)
        tf = Tf[f]
        t = T[f]
        if length(tf) == max_t_paths
            push!(Tf_adjusted, tf)
        else
            missing_tunnels = max_t_paths - length(tf)

            for i in 1:missing_tunnels
                push!(tf, t_empty)
            end
            push!(Tf_adjusted, tf)
        end
    end

    return T, Tf_adjusted, flows
end

## parse tunnels from input files (format 2)
function parsePaths(filename, links, flows; zeroindex=false)
    filename = joinpath(@__DIR__, "./data/$filename")
    nflows = length(flows)
    x = readdlm(filename)
    T = []
    Tf = [[] for i=1:nflows]
    tf = []
    fromNode = 0
    toNode = 0
    num_flow = 0
    tindex = 1
    max_paths = 0
    paths = Matrix(undef, nflows, nflows)
    for row in 1:size(x,1)
        if "->" in x[row,:]
            max_paths = max(max_paths, size(tf, 1))
            if fromNode != 0 && num_flow != 0
                #add tunnel if not first run
                paths[fromNode, toNode] = tf
                Tf[num_flow] = tf
            end
            if x[row,1] isa Number
                fromNode = x[row, 1] + zeroindex
                toNode = x[row, 3] + zeroindex
            else
                fromNode = parse(Int, replace(x[row, 1], "h" => "")) + zeroindex
                toNode = parse(Int, replace(x[row, 3], "h" => "")) + zeroindex
            end
            num_flow = findfirst(x -> x == (fromNode, toNode), flows)
            if num_flow === nothing
                num_flow = 0
            end
            tf = []
        else
            t = []
            #parse each tunnel into edges
            for col in 2:size(x,2)
                if occursin("]", x[row,col])
                    break
                end
                r = replace(x[row,col][2:end-2], "s" => "")
                stringtup = split(r, ",")
                #find edge in edge matrix
                e = (parse(Int,stringtup[1]) + zeroindex, parse(Int,stringtup[2]) + zeroindex)
                index = findfirst(x -> x == e, links)
                push!(t, index)
            end
            #create new tunnel and add index to that flows tunnels
            if num_flow != 0
                push!(T, t)
                push!(tf, tindex)
                tindex += 1
            end
        end
    end

    #add last tunnel
    if num_flow != 0
        Tf[num_flow] = tf
        paths[fromNode, toNode] = tf
    end
    push!(T, [])

    for f in 1:size(Tf,1)
        for t in 1:max_paths
            try Tf[f][t]
            catch
                push!(Tf[f], tindex)
            end
        end
    end
    return T, Tf, max_paths
end

max_stars_repo_name: hipersys-team/arrow

max_issues_repo_name: arrow-bow/arrow

max_forks_repo_name: arrow-bow/arrow


text:
using Gadfly
using DataFrames

include("./utility.jl")
function bench_cpu_regarding_n()
    N = [   1,2,3,4,5,6,7,8,9,
            10,20,30,40,50,60,70,80,90,100,200,300,400,500,750,
            1_000,2_500,5_000,7_500,10_000,100_000,1_000_000]

    cpu = Float64[]
    for n in N
        a = rand(n)
        sqrt(a)
        gc()
        gc_enable(false)
        t = mean([@elapsed sqrt(a) for i=1:100])*(1e9/n)
        # t = mean(ibench(sqrt, (a,),0.1 ; print_eval=false, print_stats=true)[1])*(1e9/n)
        gc_enable(true)
        push!(cpu,t)
        println(round(cpu[end]))
    end

    df = DataFrame()
    df[:N] = N
    df[:CPU] = cpu

    path = Pkg.dir("MKL") * "/benchmark/"
    p = Gadfly.plot(
                    layer(df,x="N",y="CPU",Geom.line),
                    Scale.x_log10,
                    Guide.xlabel("n-element vector"),
                    Guide.ylabel("CPU time in nsec/element"),
                    Guide.title("CPU time for sqrt(X) where X = Float64[] with n elements"))
    draw(PNG(path*"sqrt_cpu(n).png", 20cm, 20cm), p)
    p
end

bench_cpu_regarding_n()
@time bench_cpu_regarding_n()

begin
    N = 1_000_000
    p = convert(Int64,log10(N))
    cpu = []
    for i=0:p
        n = 10^(p-i)
        ne = div(N,n)
        a = [rand(n) for j=1:ne]
        println("ne = $ne | n = $n | nexn = ", ne*n)
        t = mean([
            @elapsed  begin
                        for j=1:ne
                            sqrt(a[j])
                        end
                    end
            for i=1:100])*(1e9/N)
        push!(cpu,t)
    end
end
println(cpu)
p = 6
df = DataFrame()
df[:N] = [10^(i) for i=0:p]
df[:CPU] = reverse(cpu)
df

p = Gadfly.plot(
                layer(df,x="N",y="CPU",Geom.line),
                Scale.x_log10,
                Guide.xlabel("n-element vector"),
                Guide.ylabel("CPU time in nsec/element"),
                Guide.title("CPU time for sqrt(X) where X = Float64[] with n elements"))

path = Pkg.dir("MKL") * "/benchmark/"
draw(PNG(path*"sqrt_frag.png", 20cm, 20cm), p)
p

function main()
    N = [1,2,3,4,5,6,7,8,9,
         10,20,30,40,50,60,70,80,90,100,200,300,400,500,750,
         1_000,2_500,5_000,7_500,10_000,100_000,1_000_000]

    cpu = Vector{Float64}(length(N))
    @inbounds for i in 1:length(N)
        n = N[i]
        a = rand(n)
        dest = rand(n)

        # mkl_sqrt!(n,a,dest)
        map!(sqrt,dest,a)
        # broadcast!(sqrt,dest,a)
        # sqrt(a)

        gc()
        tsum = 0.0
        tsum2 = 0.0

        ncycle = 10_000 ÷ n + 1
        nrep = 100

        for j in 1:nrep
            t = @elapsed for k in 1:ncycle
                # Vector{Float64}(n)
                # mkl_sqrt!(n,a,dest)
                map!(sqrt,dest,a)
                # broadcast!(sqrt,dest,a)
                # sqrt(a)
            end
            t = t / ncycle / n * 1e9
            tsum += t
            # tsum2 += t^2
        end

        tavg = tsum / nrep
        # tavg2 = tsum2 / nrep
        # tunc = sqrt((tavg2 - tavg^2) / (nrep - 1))

        # tavg = mean([@elapsed sqrt(a) for i in 1:100]) * (1e9 / n)
        cpu[i] = tavg
    end
    N, cpu
end

df = DataFrame()
df[:N], df[:CPU] = main()

p = Gadfly.plot(layer(df, x="N", y="CPU", Geom.line),
                Scale.x_log10,
                Guide.xlabel("n-element vector"),
                Guide.ylabel("CPU time in nsec/element"),
                Guide.title("CPU time with n elements")
                )

path = Pkg.dir("MKL") * "/benchmark/"
draw(PNG(path*"sqrt_map!.png", 20cm, 20cm), p)