text:
using ITensors,
      LinearAlgebra, # For tr()
      Random,        # To set a seed
      Test

Random.seed!(12345)

digits(::Type{T},x...) where {T} = T(sum([x[length(x)-k+1]*10^(k-1) for k=1:length(x)]))

@testset "Dense ITensor basic functionality" begin

@testset "ITensor constructors" begin
  i = Index(2,"i")
  j = Index(2,"j")
  k = Index(2,"k")

  @testset "Default" begin
    A = ITensor()
    @test store(A) isa Dense{Nothing}
    @test isNull(A)
  end

  @testset "Undef with index" begin
    A = ITensor(undef, i)
    @test store(A) isa Dense{Float64}
    @test !isNull(A)
  end

  @testset "Default with indices" begin
    A = ITensor(i,j)
    @test store(A) isa Dense{Float64}
    @test !isNull(A)
  end

  @testset "Random" begin
    A = randomITensor(i,j)
    @test store(A) isa Dense{Float64}

    @test ndims(A) == order(A) == 2 == length(inds(A))
    @test size(A) == dims(A) == (2,2)

    @test !isNull(A)

    B = randomITensor(IndexSet(i,j))
    @test store(B) isa Dense{Float64}
    @test ndims(B) == order(B) == 2 == length(inds(B))
    @test size(B) == dims(B) == (2,2)
    @test !isNull(B)
  end

  @testset "From matrix" begin
    M = [1 2; 3 4]
    A = ITensor(M,i,j)
    @test store(A) isa Dense{Float64}

    @test M ≈ matrix(A,i,j)
    @test M' ≈ matrix(A,j,i)
    @test_throws DimensionMismatch vector(A)

    @test size(A,1) == size(M,1) == 2
    @test size(A,3) == size(M,3) == 1
    @test_throws ErrorException size(A,0)
    @test_throws ErrorException size(M,0)

    M = [1 2 3; 4 5 6]
    @test_throws DimensionMismatch ITensor(M,i,j)
  end

  @testset "To Matrix" begin
    TM = randomITensor(i,j)

    M1 = matrix(TM)
    for ni in i, nj in j
      @test M1[ni,nj] ≈ TM[i(ni),j(nj)]
    end

    M2 = matrix(TM,j,i)
    for ni in i, nj in j
      @test M2[nj,ni] ≈ TM[i(ni),j(nj)]
    end

    T3 = randomITensor(i,j,k)
    @test_throws DimensionMismatch matrix(T3,i,j)
  end

  @testset "To Vector" begin
    TV = randomITensor(i)

    V = vector(TV)
    for ni in i
      @test V[ni] ≈ TV[i(ni)]
    end

    T2 = randomITensor(i,j)
    @test_throws DimensionMismatch vector(T2)
  end

  @testset "Complex" begin
    A = ITensor(Complex,i,j)
    @test store(A) isa Dense{ComplexF64}
  end

  @testset "Random complex" begin
    A = randomITensor(Complex,i,j)
    @test store(A) isa Dense{ComplexF64}
  end

  @testset "From complex matrix" begin
    M = [1+2im 2; 3 4]
    A = ITensor(M,i,j)
    @test store(A) isa Dense{ComplexF64}
  end

end

@testset "Convert to complex" begin
  i = Index(2,"i")
  j = Index(2,"j")
  A = randomITensor(i,j)
  B = complex(A)
  for ii ∈ dim(i), jj ∈ dim(j)
    @test complex(A[i(ii),j(jj)]) == B[i(ii),j(jj)]
  end
end

@testset "similar" begin
  i = Index(2,"i")
  j = Index(2,"j")
  A = randomITensor(i,j)
  B = similar(A)
  @test inds(B) == inds(A)
  Ac = similar(A, ComplexF32)
  @test store(Ac) isa Dense{ComplexF32}
end

@testset "fill!" begin
  i = Index(2,"i")
  j = Index(2,"j")
  A = randomITensor(i,j)
  fill!(A, 1.0)
  @test all(data(store(A)) .== 1.0)
end

@testset "copyto!" begin
  i = Index(2,"i")
  j = Index(2,"j")
  M = [1 2; 3 4]
  A = ITensor(M,i,j)
  N = 2*M
  B = ITensor(N,i,j)
  copyto!(A, B)
  @test A == B
  @test data(store(A)) == vec(N)
  A = ITensor(M,i,j)
  B = ITensor(N,j,i)
  copyto!(A, B)
  @test A == B
  @test data(store(A)) == vec(transpose(N))
end

@testset "Unary -" begin
  i = Index(2,"i")
  j = Index(2,"j")
  M = [1 2; 3 4]
  A = ITensor(M,i,j)
  @test -A == ITensor(-M, i, j)
end

@testset "dot" begin
  i = Index(2,"i")
  a = [1.0; 2.0]
  b = [3.0; 4.0]
  A = ITensor(a,i)
  B = ITensor(b,i)
  @test dot(A, B) == 11.0
end

@testset "exponentiate" begin
  s1 = Index(2,"s1")
  s2 = Index(2,"s2")
  i1 = Index(2,"i1")
  i2 = Index(2,"i2")
  Amat = rand(2,2,2,2)
  A = ITensor(Amat,i1,i2,s1,s2)

  Aexp = exp(A,(i1,i2),(s1,s2))
  Amatexp = reshape(exp(reshape(Amat,4,4)),2,2,2,2)
  Aexp_from_mat = ITensor(Amatexp,i1,i2,s1,s2)
  @test Aexp ≈ Aexp_from_mat

  #test that exponentiation works when indices need to be permuted
  Aexp = exp(A,(s1,s2),(i1,i2))
  Amatexp = Matrix(exp(reshape(Amat,4,4))')
  Aexp_from_mat = ITensor(reshape(Amatexp,2,2,2,2),s1,s2,i1,i2)
  @test Aexp ≈ Aexp_from_mat

  #test exponentiation when hermitian=true is used
  Amat = reshape(Amat, 4,4)
  Amat = reshape(Amat+Amat'+randn(4,4)*1e-10,2,2,2,2)
  A = ITensor(Amat,i1,i2,s1,s2)
  Aexp = exp(A,(i1,i2),(s1,s2),ishermitian=true)
  Amatexp = reshape(parent(exp(Hermitian(reshape(Amat,4,4)))),
                    2,2,2,2)
  Aexp_from_mat = ITensor(Amatexp,i1,i2,s1,s2)
  @test Aexp ≈ Aexp_from_mat
end

@testset "add and axpy" begin
  i = Index(2,"i")
  a = [1.0; 2.0]
  b = [3.0; 4.0]
  A = ITensor(a,i)
  B = ITensor(b,i)
  c = [5.0; 8.0]
  @test A + B == ITensor([4.0; 6.0], i)
  @test axpy!(2.0, A, B) == ITensor(c, i)
  a = [1.0; 2.0]
  b = [3.0; 4.0]
  A = ITensor(a,i)
  B = ITensor(b,i)
  c = [8.0; 12.0]
  @test add!(A, 2.0, 2.0, B) == ITensor(c, i)

end

@testset "mul! and rmul!" begin
  i = Index(2,"i")
  a = [1.0; 2.0]
  b = [2.0; 4.0]
  A = ITensor(a,i)
  A2, A3 = copy(A), copy(A)
  B = ITensor(b,i)
  @test mul!(A2, A, 2.0) == B == ITensors.add!(A2, 0, 2, A)
  @test rmul!(A, 2.0) == B == ITensors.scale!(A3, 2)

  i = Index(2,"i")
  j = Index(2,"j")
  M = [1 2; 3 4]
  A = ITensor(M,i,j)
  N = 2*M
  B = ITensor(N,j,i)
  @test data(store(mul!(B, A, 2.0))) == 2.0*vec(transpose(M))
end

@testset "show" begin
  i = Index(2,"i")
  a = [1.0; 2.0]
  A = ITensor(a,i)
  s = split(sprint(show, A), '\n')
  @test s[1] == "ITensor ord=1 " * sprint(show, i) * " "
  @test s[2] == "Dense{Float64}"
  @test s[3] == " 1.0"
  @test s[4] == " 2.0"
end

@testset "Test isapprox for ITensors" begin
  m,n = rand(0:20,2)
  i = Index(m)
  j = Index(n)
  realData = rand(m,n)
  complexData = complex(realData)
  A = ITensor(realData, i,j)
  B = ITensor(complexData, i,j)
  @test A≈B
  @test B≈A
  A = permute(A,j,i)
  @test A≈B
  @test B≈A
end

@testset "ITensor tagging and priming" begin
  s1 = Index(2,"Site,s=1")
  s2 = Index(2,"Site,s=2")
  l = Index(3,"Link")
  A1 = randomITensor(s1,l,l')
  A2 = randomITensor(s2,l',l'')
  @testset "findindex(::ITensor,::String)" begin
    @test s1==findindex(A1,"Site")
    @test s1==findindex(A1,"s=1")
    @test s1==findindex(A1,"s=1,Site")
    @test l==findindex(A1,("Link",0))
    @test l'==findindex(A1,("",1))
    @test l'==findindex(A1,("Link",1))
    @test s2==findindex(A2,"Site")
    @test s2==findindex(A2,"s=2")
    @test s2==findindex(A2,"Site")
    @test s2==findindex(A2,("",0))
    @test s2==findindex(A2,("s=2",0))
    @test s2==findindex(A2,("Site",0))
    @test s2==findindex(A2,("s=2,Site",0))
    @test l'==findindex(A2,("",1))
    @test l'==findindex(A2,("Link",1))
    @test l''==findindex(A2,("",2))
    @test l''==findindex(A2,("Link",2))
  end
  @testset "addtags(::ITensor,::String,::String)" begin
    s1u = addtags(s1,"u")
    lu = addtags(l,"u")

    A1u = addtags(A1,"u")
    @test hasinds(A1u,s1u,lu,lu')

    A1u = addtags(A1,"u","Link")
    @test hasinds(A1u,s1,lu,lu')

    A1u = addtags(A1,"u",("",0))
    @test hasinds(A1u,s1u,lu,l')

    A1u = addtags(A1,"u",("Link",0))
    @test hasinds(A1u,s1,lu,l')

    A1u = addtags(A1,"u",("Link",1))
    @test hasinds(A1u,s1,l,lu')
  end
  @testset "removetags(::ITensor,::String,::String)" begin
    A2r = removetags(A2,"Site")
    @test hasinds(A2r,removetags(s2,"Site"),l',l'')

    A2r = removetags(A2,"Link",("",1))
    @test hasinds(A2r,s2,removetags(l,"Link")',l'')
  end
  @testset "replacetags(::ITensor,::String,::String)" begin
    s2tmp = replacetags(s2,"Site","Temp")
    ltmp = replacetags(l,"Link","Temp")

    A2r = replacetags(A2,"Site","Temp")
    @test hasinds(A2r,s2tmp,l',l'')

    A2r = replacetags(A2,"Link","Temp")
    @test hasinds(A2r,s2,ltmp',ltmp'')

    A2r = replacetags(A2,"Link","Temp",("",1))
    @test hasinds(A2r,s2,ltmp',l'')

    A2r = replacetags(A2,("Link",2),("Temp",3))
    @test hasinds(A2r,s2,l',ltmp''')

    A2r = replacetags(A2,("",1),("",5))
    @test hasinds(A2r,s2,prime(l,5),l'')

    #In-place version
    cA2 = copy(A2)
    replacetags!(cA2,("",1),("",5))
    @test hasinds(cA2,s2,prime(l,5),l'')
  end
  @testset "prime(::ITensor,::String)" begin
    A2p = prime(A2)
    @test A2p==A2'
    @test hasinds(A2p,s2',l'',l''')

    A2p = prime(A2,2)
    A2p = A2''
    @test hasinds(A2p,s2'',l''',l'''')

    A2p = prime(A2,"s=2")
    @test hasinds(A2p,s2',l',l'')
  end

  @testset "mapprime" begin
    @test hasinds(mapprime(A2,1,7),s2,l^7,l'')
    @test hasinds(mapprime(A2,0,1),s2',l',l'')
  end
  @testset "setprime" begin
    @test hasinds(setprime(A2,2,s2),s2'',l',l'')
    @test hasinds(setprime(A2,0,l''),s2,l',l)
  end
  @testset "swapprime" begin
    @test hasinds(swapprime(A2,1,3),l''',s2,l'')
  end
end

@testset "Converting Real and Complex Storage" begin

  @testset "Add Real and Complex" begin
    i = Index(2,"i")
    j = Index(2,"j")
    TC = randomITensor(ComplexF64,i,j)
    TR = randomITensor(Float64,i,j)

    S1 = TC+TR
    S2 = TR+TC
    @test typeof(S1.store) == Dense{ComplexF64}
    @test typeof(S2.store) == Dense{ComplexF64}
    for ii=1:dim(i),jj=1:dim(j)
      @test S1[i(ii),j(jj)] ≈ TC[i(ii),j(jj)]+TR[i(ii),j(jj)]
      @test S2[i(ii),j(jj)] ≈ TC[i(ii),j(jj)]+TR[i(ii),j(jj)]
    end
  end

end

@testset "ITensor, Dense{$SType} storage" for SType ∈ (Float64,ComplexF64)
  mi,mj,mk,ml,mα = 2,3,4,5,6,7
  i = Index(mi,"i")
  j = Index(mj,"j")
  k = Index(mk,"k")
  l = Index(ml,"l")
  α = Index(mα,"alpha")
  @testset "Set and get values with IndexVals" begin
    A = ITensor(SType,i,j,k)
    for ii ∈ 1:dim(i), jj ∈ 1:dim(j), kk ∈ 1:dim(k)
      A[k(kk),j(jj),i(ii)] = digits(SType,ii,jj,kk)
    end
    for ii ∈ 1:dim(i), jj ∈ 1:dim(j), kk ∈ 1:dim(k)
      @test A[j(jj),k(kk),i(ii)]==digits(SType,ii,jj,kk)
    end
    @test_throws MethodError A[1]
  end
  @testset "Test permute(ITensor,Index...)" begin
    A = randomITensor(SType,i,k,j)
    permA = permute(A,k,j,i)
    @test k==inds(permA)[1]
    @test j==inds(permA)[2]
    @test i==inds(permA)[3]
    for ii ∈ 1:dim(i), jj ∈ 1:dim(j), kk ∈ 1:dim(k)
      @test A[k(kk),i(ii),j(jj)]==permA[i(ii),j(jj),k(kk)]
    end
    for ii ∈ 1:dim(i), jj ∈ 1:dim(j), kk ∈ 1:dim(k)
      @test A[k(kk),i(ii),j(jj)]==permA[i(ii),j(jj),k(kk)]
    end
    # TODO: I think this was doing slicing, but what is the output
    # of slicing an ITensor?
    #@testset "getindex and setindex with vector of IndexVals" begin
    #    k_inds = [k(kk) for kk ∈ 1:dim(k)]
    #    for ii ∈ 1:dim(i), jj ∈ 1:dim(j)
    #      @test A[k_inds,i(ii),j(jj)]==permA[i(ii),j(jj),k_inds...]
    #    end
    #    for ii ∈ 1:dim(i), jj ∈ 1:dim(j)
    #        A[k_inds,i(ii),j(jj)]=collect(1:length(k_inds))
    #    end
    #    permA = permute(A,k,j,i)
    #    for ii ∈ 1:dim(i), jj ∈ 1:dim(j)
    #      @test A[k_inds,i(ii),j(jj)]==permA[i(ii),j(jj),k_inds...]
    #    end
    #end
  end
  @testset "Set and get values with Ints" begin
    A = ITensor(SType,i,j,k)
    A = permute(A,k,i,j)
    for ii ∈ 1:dim(i), jj ∈ 1:dim(j), kk ∈ 1:dim(k)
      A[kk,ii,jj] = digits(SType,ii,jj,kk)
    end
    A = permute(A,i,j,k)
    for ii ∈ 1:dim(i), jj ∈ 1:dim(j), kk ∈ 1:dim(k)
      @test A[ii,jj,kk]==digits(SType,ii,jj,kk)
    end
  end
  @testset "Test scalar(ITensor)" begin
    x = SType(34)
    A = ITensor(x)
    @test x==scalar(A)
    A = ITensor(SType,i,j,k)
    @test_throws BoundsError scalar(A)
  end
  @testset "Test norm(ITensor)" begin
    A = randomITensor(SType,i,j,k)
    @test norm(A)≈sqrt(scalar(dag(A)*A))
  end
  @testset "Test add ITensors" begin
    A = randomITensor(SType,i,j,k)
    B = randomITensor(SType,k,i,j)
    C = A+B
    for ii ∈ 1:dim(i), jj ∈ 1:dim(j), kk ∈ 1:dim(k)
      @test C[i(ii),j(jj),k(kk)]==A[j(jj),i(ii),k(kk)]+B[i(ii),k(kk),j(jj)]
    end
    @test array(permute(C,i,j,k))==array(permute(A,i,j,k))+array(permute(B,i,j,k))
  end

  @testset "Test factorizations of an ITensor" begin

    A = randomITensor(SType,i,j,k,l)

    @testset "Test SVD of an ITensor" begin
      U,S,V,u,v = svd(A,(j,l))
      @test store(S) isa Diag{Vector{Float64}}
      @test A≈U*S*V
      @test U*dag(prime(U,u))≈δ(SType,u,u') atol=1e-14
      @test V*dag(prime(V,v))≈δ(SType,v,v') atol=1e-14
    end

    @testset "Test SVD truncation" begin
        ii = Index(4)
        jj = Index(4)
        T = randomITensor(ComplexF64,ii,jj)
        U,S,V = svd(T,ii;maxdim=2)
        u,s,v = svd(matrix(T))
        @test norm(U*S*V-T)≈sqrt(s[3]^2+s[4]^2)
    end

    @testset "Test QR decomposition of an ITensor" begin
      Q,R,q = qr(A,(i,l))
      q = commonindex(Q,R)
      @test A≈Q*R
      @test Q*dag(prime(Q,q))≈δ(SType,q,q') atol=1e-14
    end

    @testset "Test polar decomposition of an ITensor" begin
      U,P,u = polar(A,(k,l))
      @test A≈U*P
      #Note: this is only satisfied when left dimensions
      #are greater than right dimensions
      UUᵀ =  U*dag(prime(U,u))

      # TODO: use a combiner to combine the u indices to make
      # this test simpler
      for ii ∈ 1:dim(u[1]), jj ∈ 1:dim(u[2]), iip ∈ 1:dim(u[1]), jjp ∈ 1:dim(u[2])
        val = UUᵀ[u[1](ii),u[2](jj),u[1]'(iip),u[2]'(jjp)]
        if ii==iip && jj==jjp
          @test val≈one(SType) atol=1e-14
        else
          @test val≈zero(SType) atol=1e-14
        end
      end
    end

    @testset "Test Hermitian eigendecomposition of an ITensor" begin
      is = IndexSet(i,j)
      T = randomITensor(is...,prime(is)...)
      T = T + swapprime(dag(T),0,1)
      U,D,u = eigenHermitian(T)
      @test T ≈ U*D*prime(dag(U))
      UUᴴ =  U*prime(dag(U),u)
      @test UUᴴ ≈ δ(u,u') atol=1e-14
    end

  end # End ITensor factorization testset

end # End Dense storage test

end # End Dense ITensor basic functionality


text:
function get_sim_data_indirect(mechanism::Mechanism,env::Environment,Δt::Real;
                               relax_comp=false)
    vs = VariableSelector()
    add_var!(vs, :qnext, num_positions(mechanism))
    add_var!(vs, :vnext, num_velocities(mechanism))
    for i = 1:length(env.contacts)
        β_dim = size(env.contacts[i].obstacle.basis,2)
        add_var!(vs, Symbol("β", i), β_dim)
        add_var!(vs, Symbol("λ", i), 1)
        add_var!(vs, Symbol("c_n", i), 1)
        if relax_comp
            add_var!(vs, Symbol("slack", i), 3)
        end
    end

    cs = ConstraintSelector()
    add_eq!(cs, :kin, num_positions(mechanism))
    add_eq!(cs, :dyn, num_velocities(mechanism))
    for i = 1:length(env.contacts)
        β_dim = size(env.contacts[i].obstacle.basis,2)
        add_ineq!(cs, Symbol("β_pos", i), β_dim)
        add_ineq!(cs, Symbol("λ_pos", i), 1)
        add_ineq!(cs, Symbol("c_n_pos", i), 1)
        add_ineq!(cs, Symbol("ϕ_pos", i), 1)
        add_ineq!(cs, Symbol("fric_pos", i), β_dim)
        add_ineq!(cs, Symbol("cone_pos", i), 1)
        if relax_comp
            add_ineq!(cs, Symbol("slack_pos", i), 3)
            add_ineq!(cs, Symbol("ϕ_c_n_comp", i), 1)
            add_ineq!(cs, Symbol("fric_β_comp", i), β_dim)
            add_ineq!(cs, Symbol("cone_λ_comp", i), 1)
        else
            add_eq!(cs, Symbol("ϕ_c_n_comp", i), 1)
            add_eq!(cs, Symbol("fric_β_comp", i), β_dim)
            add_eq!(cs, Symbol("cone_λ_comp", i), 1)
        end
    end

    state_cache = [StateCache(mechanism) for n = 1:2]
    envj_cache = [EnvironmentJacobianCache(env) for n = 1:2]

    generate_solver_fn = :generate_solver_fn_sim_indirect
    extract_sol = :extract_sol_sim_indirect

    sim_data = SimData(mechanism,env,
                       state_cache,envj_cache,
                       Δt,vs,cs,generate_solver_fn,extract_sol,
                       [],[],[],[],[],[],[],[],[],1,[],[],[],[])

    sim_data
end

function extract_sol_sim_indirect(sim_data::SimData, results::AbstractArray{T,2}) where T
    vs = sim_data.vs
    relax_comp = haskey(vs.vars, :slack1)
    env = sim_data.env
    N = size(results,2)

    qtraj = Array{Array{Float64,1},1}(undef, 0)
    vtraj = Array{Array{Float64,1},1}(undef, 0)
    utraj = Array{Array{Float64,1},1}(undef, 0)
    contact_traj = Array{Array{Array{Float64,1},1},1}(undef, 0)
    slack_traj = Array{Array{Float64,1},1}(undef, 0)
    x = MechanismState(sim_data.mechanism)
    for n = 1:N
        set_configuration!(x, vs(results[:,n], Symbol("qnext")))
        normalize_configuration!(x)
        qnext = configuration(x)
        push!(qtraj, qnext)
        push!(vtraj, vs(results[:,n], Symbol("vnext")))
        contact_sol = Array{Array{Float64,1},1}(undef, 0)
        for i = 1:length(env.contacts)
            if relax_comp
                push!(slack_traj, vs(results[:,n], Symbol("slack", i)))
            end
            push!(contact_sol, vcat(vs(results[:,n], Symbol("c_n", i)),
                               vs(results[:,n], Symbol("β", i)),
                               vs(results[:,n], Symbol("λ", i))))
        end
        push!(contact_traj, contact_sol)
    end

    # some other usefull vectors
    ttraj = [(i-1)*sim_data.Δt for i = 1:N]
    qv_mat = vcat(hcat(qtraj...),hcat(vtraj...))

    qtraj, vtraj, utraj, contact_traj, slack_traj, ttraj, qv_mat, results
end

function generate_solver_fn_sim_indirect(sim_data,q0,v0,u0)
    x0 = sim_data.state_cache[1][Float64]
    Δt = sim_data.Δt
    vs = sim_data.vs
    cs = sim_data.cs

    relax_comp = haskey(vs.vars, :slack1)
    num_contacts = length(sim_data.env.contacts)
    num_vel = num_velocities(sim_data.mechanism)
    world = root_body(sim_data.mechanism)
    world_frame = default_frame(world)

    set_configuration!(x0, q0)
    set_velocity!(x0, v0)
    setdirty!(x0)
    H = mass_matrix(x0)

    function eval_obj(x::AbstractArray{T}) where T
        f = 0.

        if relax_comp
            for i = 1:num_contacts
                slack = vs(x, Symbol("slack", i))
                f += sum(slack)
            end
        end

        f
    end

    function eval_cons(x::AbstractArray{T}) where T
        xn = sim_data.state_cache[2][T]
        envj = sim_data.envj_cache[2][T]

        contact_bias = Vector{T}(undef, num_vel)
        g = Vector{T}(undef, cs.num_eqs + cs.num_ineqs) # TODO preallocate

        qnext = vs(x, :qnext)
        vnext = vs(x, :vnext)

        set_configuration!(xn, qnext)
        set_velocity!(xn, vnext)
        setdirty!(xn)

        normalize_configuration!(xn)

        config_derivative = configuration_derivative(xn) # TODO preallocate
        dyn_bias = dynamics_bias(xn) # TODO preallocate
        if (num_contacts > 0)
            contact_jacobian!(envj, xn)
            contact_τ_indirect!(contact_bias, sim_data, envj, x)
        end

        g[cs(:kin)] .= qnext .- q0 .- Δt .* config_derivative
        g[cs(:dyn)] .= H * (vnext - v0) .- Δt .* (u0 .- dyn_bias .- contact_bias)
        for i = 1:num_contacts
            cj = envj.contact_jacobians[i]
            c = cj.contact

            # TODO preallocate
            contact_v = cj.contact_rot' * point_jacobian(xn, path(sim_data.mechanism, world, c.body), transform(xn, c.point, world_frame)).J * vnext

            β = vs(x, Symbol("β", i))
            λ = vs(x, Symbol("λ", i))
            c_n = vs(x, Symbol("c_n", i))

            Dtv = c.obstacle.basis' * contact_v

            g[cs(Symbol("β_pos", i))] .= -β
            g[cs(Symbol("λ_pos", i))] .= -λ
            g[cs(Symbol("c_n_pos", i))] .= -c_n
            g[cs(Symbol("ϕ_pos", i))] .= -envj.contact_jacobians[i].ϕ
            # λe + D'*v >= 0
            g[cs(Symbol("fric_pos", i))] .= -(λ .+ Dtv)
            # μ*c_n - sum(β) >= 0
            g[cs(Symbol("cone_pos", i))] .= -(c.obstacle.μ .* c_n .- sum(β))
            if !relax_comp
                # dist * c_n = 0
                g[cs(Symbol("ϕ_c_n_comp", i))] .= envj.contact_jacobians[i].ϕ .* c_n
                # (λe + Dtv)' * β = 0
                g[cs(Symbol("fric_β_comp", i))] .= (λ .+ Dtv) .* β
                # (μ * c_n - sum(β)) * λ = 0
                g[cs(Symbol("cone_λ_comp", i))] .= (c.obstacle.μ .* c_n .- sum(β)) .* λ
            else
                slack = vs(x, Symbol("slack", i))
                g[cs(Symbol("slack_pos", i))] .= -slack
                g[cs(Symbol("ϕ_c_n_comp", i))] .= envj.contact_jacobians[i].ϕ .* c_n .- slack[1]
                g[cs(Symbol("fric_β_comp", i))] .= (λ .+ Dtv) .* β .- slack[2]
                g[cs(Symbol("cone_λ_comp", i))] .= (c.obstacle.μ .* c_n .- sum(β)) .* λ .- slack[3]
            end
        end

        g
    end

    generate_autodiff_solver_fn(eval_obj,eval_cons,cs.eqs,cs.ineqs,vs.num_vars)
end

ext: jl


text:
function set_boundary(Qbase, cellxmax, cellymax, vecAx, vecAy, bdcon, Rd, g, nval, boundary_number, rank, size, comm)

    divn = zeros(Int, 4)
    temp = Array{Float64}(undef,2)

    # 2:imaginary cell
    icell = 2
    send_arrayx = zeros(cellxmax*icell*nval)
    send_arrayy = zeros(cellymax*icell*nval)
    rec_arrayx = zeros(cellxmax*icell*nval)
    rec_arrayy = zeros(cellymax*icell*nval)

    #println(boundary_number)

    # x-
    if boundary_number[1] == -41
        # parallel with communication

        ite = 1
        for l in 1:nval
            for j in 1:cellymax
                for i in icell+1:icell+icell  # 周期境界条件
                    send_arrayy[ite] = Qbase[i,j,l]
                    ite += 1
                end
            end
        end

        # 送る
        send_rank = rank + (size^0.5-1) * (size^0.5)
        MPI.Isend(send_arrayy, send_rank, rank+200, comm)

        # cellから受け取る
        rec = MPI.Irecv!(rec_arrayy, send_rank, send_rank+100, comm)

        # wait
        MPI.Waitall!([rec])

        # 代入
        ite = 1
        for l in 1:nval
            for j in 1:cellymax
                for i in icell+1:icell+icell
                    Qbase[i,j,l] = rec_arrayy[ite]
                    ite += 1
                end
            end
        end

    elseif boundary_number[1] < 0
        Qbase = boundary_condition_xm(Qbase, cellxmax, cellymax, vecAx, vecAy, bdcon, Rd, g, nval)
    else
        # parallel with communication

        ite = 1
        for l in 1:nval
            for j in 1:cellymax
                for i in 1:icell
                    send_arrayy[ite] = Qbase[i,j,l]
                    ite += 1
                end
            end
        end

        # south に送る
        MPI.Isend(send_arrayy, boundary_number[2], rank+200, comm)

        # north of south cell から受け取る
        rec = MPI.Irecv!(rec_arrayy, boundary_number[1], boundary_number[1]+100, comm)

        # wait
        MPI.Waitall!([rec])

        # 代入
        ite = 1
        for l in 1:nval
            for j in 1:cellymax
                for i in 1:icell
                    Qbase[i,j,l] = rec_arrayy[ite]
                    ite += 1
                end
            end
        end

        # use ?
        divn[1] = 1
    end

    # x+
    if  boundary_number[2] == -41
        # parallel with communication

        ite = 1
        for l in 1:nval
            for j in 1:cellymax
                for i in cellxmax-icell-(icell+1):cellxmax-icell  # 周期境界条件
                    send_arrayy[ite] = Qbase[i,j,l]
                    ite += 1
                end
            end
        end

        # 送る
        send_rank = rank - (size^0.5-1) * (size^0.5)
        MPI.Isend(send_arrayy, send_rank, rank+200, comm)

        # cellから受け取る
        rec = MPI.Irecv!(rec_arrayy, send_rank, send_rank+100, comm)

        # wait
        MPI.Waitall!([rec])

        # 代入
        ite = 1
        for l in 1:nval
            for j in 1:cellymax
                for i in cellxmax-icell-(icell+1):cellxmax-icell  # 周期境界条件
                    Qbase[i,j,l] = rec_arrayy[ite]
                    ite += 1
                end
            end
        end

    elseif boundary_number[2] < 0
        Qbase = boundary_condition_xp(Qbase, cellxmax, cellymax, vecAx, vecAy, bdcon, Rd, g, nval)
    else
        # parallel with communication

        ite = 1
        for l in 1:nval
            for j in 1:cellymax
                for i in cellxmax-icell+1:cellxmax
                    send_arrayy[ite] = Qbase[i,j,l]
                    ite += 1
                end
            end
        end

        # south に送る
        MPI.Isend(send_arrayy, boundary_number[1], rank+100, comm)

        # north of south cell から受け取る
        rec = MPI.Irecv!(rec_arrayy, boundary_number[2], boundary_number[2]+200, comm)

        # wait
        MPI.Waitall!([rec])

        # 代入
        ite = 1
        for l in 1:nval
            for j in 1:cellymax
                for i in cellxmax-icell+1:cellxmax
                    Qbase[i,j,l] = rec_arrayy[ite]
                    ite += 1
                end
            end
        end

        # use ?
        divn[2] = 1
    end

    # y-
    if  boundary_number[3] == -41
        # parallel with communication

        ite = 1
        for l in 1:nval
            for j in icell+1:icell+icell  # 周期境界条件
                for i in 1:cellxmax
                    send_arrayy[ite] = Qbase[i,j,l]
                    ite += 1
                end
            end
        end

        # 送る
        send_rank = rank + (size^0.5-1)
        MPI.Isend(send_arrayy, send_rank, rank+200, comm)

        # cellから受け取る
        rec = MPI.Irecv!(rec_arrayy, send_rank, send_rank+100, comm)

        # wait
        MPI.Waitall!([rec])

        # 代入
        ite = 1
        for l in 1:nval
            for j in icell+1:icell+icell  # 周期境界条件
                for i in 1:cellxmax
                    Qbase[i,j,l] = rec_arrayy[ite]
                    ite += 1
                end
            end
        end

    elseif boundary_number[3] < 0
        Qbase = boundary_condition_ym(Qbase, cellxmax, cellymax, vecAx, vecAy, bdcon, Rd, g, nval)
    else
        # parallel with communication

        ite = 1
        for l in 1:nval
            for j in 1:icell
                for i in 1:cellymax
                    send_arrayy[ite] = Qbase[i,j,l]
                    ite += 1
                end
            end
        end

        # south に送る
        MPI.Isend(send_arrayy, boundary_number[4], rank+400, comm)

        # north of south cell から受け取る
        rec = MPI.Irecv!(rec_arrayy, boundary_number[3], boundary_number[3]+300, comm)

        # wait
        MPI.Waitall!([rec])

        ite = 1
        for l in 1:nval
            for j in 1:icell
                for i in 1:cellxmax
                    Qbase[i,j,l] = send_arrayy[ite]
                    ite += 1
                end
            end
        end

        # use ?
        divn[3] = 1
    end

    # y+
    if boundary_number[4] == -41

        ite = 1
        for l in 1:nval
            for j in cellymax-icell-(icell+1):cellymax-icell  # 周期境界条件
                for i in 1:cellxmax
                    send_arrayy[ite] = Qbase[i,j,l]
                    ite += 1
                end
            end
        end

        # 送る
        send_rank = rank - (size^0.5-1)
        MPI.Isend(send_arrayy, send_rank, rank+200, comm)

        # cellから受け取る
        rec = MPI.Irecv!(rec_arrayy, send_rank, send_rank+100, comm)

        # wait
        MPI.Waitall!([rec])

        # 代入
        ite = 1
        for l in 1:nval
            for j in cellymax-icell-(icell+1):cellymax-icell  # 周期境界条件
                for i in 1:cellxmax
                    Qbase[i,j,l] = rec_arrayy[ite]
                    ite += 1
                end
            end
        end
    elseif boundary_number[4] < 0
        Qbase = boundary_condition_yp(Qbase, cellxmax, cellymax, vecAx, vecAy, bdcon, Rd, g, nval)
    else
        # parallel with communication

        ite = 1
        for l in 1:nval
            for j in cellymax-icell+1:cellymax
                for i in 1:cellymax
                    send_arrayy[ite] = Qbase[i,j,l]
                    ite += 1
                end
            end
        end

        # south に送る
        MPI.Isend(send_arrayy, boundary_number[3], rank+300, comm)

        # north of south cell から受け取る
        rec = MPI.Irecv!(rec_arrayy, boundary_number[4], boundary_number[4]+400, comm)

        # wait
        MPI.Waitall!([rec])

        ite = 1
        for l in 1:nval
            for j in cellymax-icell+1:cellymax
                for i in 1:cellxmax
                    Qbase[i,j,l] = send_arrayy[ite]
                    ite += 1
                end
            end
        end

        # use ?
        divn[4] = 1
    end

    return Qbase
end

# ------------------------------------
# boundary conditions x- (x[1] & x[2])
# ------------------------------------
function boundary_condition_xm(Qbase, cellxmax, cellymax, vecAx, vecAy, bdcon, Rd, g, nval)
    """
    bdcon[i][j]
    i : bd number(x-, x+, y-, y+)
    j=1-6 : "bd1_con":"2",
            "bd1_rho":"1.0",
            "bd1_u"  :"300.0",
            "bd1_v"  :"0.0",
            "bd1_p"  :"1.0",
            "bd1_T"  :"300.0",
    """

    # bd1 = x-
    if Int(bdcon[1][1]) == 11
        for j in 1:cellymax
            for l in 1:nval
                Qbase[1,j,l] = bdcon[1][l+1]
            end
        end
    elseif Int(bdcon[1][1]) == 21
        for j in 1:cellymax
            for l in 1:nval
                Qbase[1,j,l] = Qbase[2,j,l]
            end
        end
    elseif Int(bdcon[1][1]) == 31
        for j in 1:cellymax
            for l in 1:nval
                Qbase[1,j,l] = Qbase[2,j,l]
            end

            xvec = vecAx[2,j,1]
            yvec = vecAx[2,j,2]
            u = Qbase[2,j,2]
            v = Qbase[2,j,3]

            Qbase[1,j,2] = ((-xvec^2+yvec^2)*u-2*xvec*yvec*v)/(xvec^2+yvec^2)
            Qbase[1,j,3] = (-2*xvec*yvec*u+(xvec^2-yvec^2)*v)/(xvec^2+yvec^2)
        end
    elseif Int(bdcon[1][1]) == 32
        for j in 1:cellymax
            for l in 1:nval
                Qbase[1,j,l] = Qbase[2,j,l]
            end

            u = Qbase[2,j,2]
            v = Qbase[2,j,3]

            Qbase[1,j,2] = -u
            Qbase[1,j,3] = -v
        end
    elseif Int(bdcon[1][1]) == 41
        for j in 1:cellymax
            for l in 1:nval
                Qbase[1,j,l] = Qbase[cellxmax-1,j,l]
            end
        end
    else
        println("------------------------")
        println(" boundary error ")
        println(" please check boundary number ")
        println("------------------------")
        throw(UndefVarError(:x))
    end

    return Qbase
end

# ------------------------------------
# boundary conditions x+ (x[cellxmax-1] & x[cellxmax])
# ------------------------------------
function boundary_condition_xp(Qbase, cellxmax, cellymax, vecAx, vecAy, bdcon, Rd, g, nval)

    # bd2 = x+
    if Int(bdcon[2][1]) == 11
        for j in 1:cellymax
            for l in 1:nval
                Qbase[cellxmax,j,l] = bdcon[2][l+1]
            end
        end
    elseif Int(bdcon[2][1]) == 21
        for j in 1:cellymax
            for l in 1:nval
                Qbase[cellxmax,j,l] = Qbase[cellxmax-1,j,l]
            end
        end
    elseif Int(bdcon[2][1]) == 31
        for j in 1:cellymax
            for l in 1:nval
                Qbase[cellxmax,j,l] = Qbase[cellxmax-1,j,l]
            end

            xvec = vecAx[cellymax,j,1]
            yvec = vecAx[cellymax,j,2]
            u = Qbase[cellxmax-1,j,2]
            v = Qbase[cellxmax-1,j,3]

            Qbase[cellxmax,j,2] = ((-xvec^2+yvec^2)*u-2*xvec*yvec*v)/(xvec^2+yvec^2)
            Qbase[cellxmax,j,3] = (-2*xvec*yvec*u+(xvec^2-yvec^2)*v)/(xvec^2+yvec^2)
        end
    elseif Int(bdcon[2][1]) == 32
        for j in 1:cellymax
            for l in 1:nval
                Qbase[cellxmax,j,l] = Qbase[cellxmax-1,j,l]
            end
            u = Qbase[cellxmax-1,j,2]
            v = Qbase[cellxmax-1,j,3]

            Qbase[cellxmax,j,2] = -u
            Qbase[cellxmax,j,3] = -v
        end
    elseif Int(bdcon[2][1]) == 41
        for j in 1:cellymax
            for l in 1:nval
                Qbase[cellxmax,j,l] = Qbase[2,j,l]
            end
        end
    else
        println("------------------------")
        println(" boundary error ")
        println(" please check boundary number ")
        println("------------------------")
        throw(UndefVarError(:x))
    end

    return Qbase
end

# ------------------------------------
# boundary conditions y- (y[1] & y[2])
# ------------------------------------
function boundary_condition_ym(Qbase, cellxmax, cellymax, vecAx, vecAy, bdcon, Rd, g, nval)

    # bd3 = y-
    if Int(bdcon[3][1]) == 11
        for i in 1:cellxmax
            for l in 1:nval
                Qbase[i,1,l] = bdcon[3][l+1]
            end
        end
    elseif Int(bdcon[3][1]) == 21
        for i in 1:cellxmax
            for l in 1:nval
                Qbase[i,1,l] = Qbase[i,2,l]
            end
        end
    elseif Int(bdcon[3][1]) == 31
        for i in 1:cellxmax
            for l in 1:nval
                Qbase[i,1,l] = Qbase[i,2,l]
            end

            xvec = vecAy[i,2,1]
            yvec = vecAy[i,2,2]
            u = Qbase[i,2,2]
            v = Qbase[i,2,3]

            Qbase[i,1,2] = ((-xvec^2+yvec^2)*u-2*xvec*yvec*v)/(xvec^2+yvec^2)
            Qbase[i,1,3] = (-2*xvec*yvec*u+(xvec^2-yvec^2)*v)/(xvec^2+yvec^2)
        end
    elseif Int(bdcon[3][1]) == 32
        for i in 1:cellxmax
            for l in 1:nval
                Qbase[i,1,l] = Qbase[i,2,l]
            end

            u = Qbase[i,2,2]
            v = Qbase[i,2,3]

            Qbase[i,1,2] = -u
            Qbase[i,1,3] = -v
        end
    elseif Int(bdcon[3][1]) == 41
        for i in 1:cellxmax
            for l in 1:nval
                Qbase[i,1,l] = Qbase[i,cellymax-1,l]
            end
        end
    elseif Int(bdcon[3][1]) == 33
        for i in 1:cellxmax
            for l in 1:nval
                Qbase[i,1,l] = Qbase[i,2,l]
            end

            u = Qbase[i,2,2]
            v = Qbase[i,2,3]

            Qbase[i,1,2] = -u
            Qbase[i,1,3] = -v

            Tw = bdcon[3][nval+2]
            p  = Qbase[i,1,4]
            Qbase[i,1,1] = p/(Tw*Rd)
        end
    else
        println("------------------------")
        println(" boundary error ")
        println(" please check boundary number ")
        println("------------------------")
        throw(UndefVarError(:x))
    end

    return Qbase
end

# ------------------------------------
# boundary conditions y+ (y[cellymax-1] & y[cellymax])
# ------------------------------------
function boundary_condition_yp(Qbase, cellxmax, cellymax, vecAx, vecAy, bdcon, Rd, g, nval)

    # bd4 = y+
    if Int(bdcon[4][1]) == 11
        for i in 1:cellxmax
            for l in 1:nval
                Qbase[i,cellymax,l] = bdcon[4][l+1]
            end
        end
    elseif Int(bdcon[4][1]) == 21
        for i in 1:cellxmax
            for l in 1:nval
                Qbase[i,cellymax,l] = Qbase[i,cellymax-1,l]
            end
        end
    elseif Int(bdcon[4][1]) == 31
        for i in 1:cellxmax
            for l in 1:nval
                Qbase[i,cellymax,l] = Qbase[i,cellymax-1,l]
            end

            xvec = vecAy[i,cellymax,1]
            yvec = vecAy[i,cellymax,2]
            u = Qbase[i,cellymax-1,2]
            v = Qbase[i,cellymax-1,3]

            Qbase[i,cellymax,2] = ((-xvec^2+yvec^2)*u-2*xvec*yvec*v)/(xvec^2+yvec^2)
            Qbase[i,cellymax,3] = (-2*xvec*yvec*u+(xvec^2-yvec^2)*v)/(xvec^2+yvec^2)
        end
    elseif Int(bdcon[4][1]) == 32
        for i in 1:cellxmax
            for l in 1:nval
                Qbase[i,cellymax,l] = Qbase[i,cellymax-1,l]
            end

            u = Qbase[i,cellymax-1,2]
            v = Qbase[i,cellymax-1,3]

            Qbase[i,cellymax,2] = -u
            Qbase[i,cellymax,3] = -v
        end
    elseif Int(bdcon[4][1]) == 41
        for i in 1:cellxmax
            for l in 1:nval
                Qbase[i,cellymax,l] = Qbase[i,2,l]
            end
        end
    elseif Int(bdcon[4][1]) == 101
        temp1 = Int64(round(cellxmax/4)+1) # 1/4流出
        temp2 = Int64((temp1-1)*3+1)       # 1/4~3/4流入
        for i in 1:temp1
            for l in 1:nval
                Qbase[i,cellymax,l] = Qbase[i,cellymax-1,l]
            end
        end
        for i in temp1+1:temp2
            for l in 1:nval
                Qbase[i,cellymax,l] = bdcon[4][l+1]
            end
        end
        for i in temp2+1:cellxmax
            for l in 1:nval
                Qbase[i,cellymax,l] = Qbase[i,cellymax-1,l]
            end
        end
    elseif Int(bdcon[4][1]) == 102
        temp1 = Int64(round(cellxmax/4)+1) # 1/4流出
        temp2 = Int64((temp1-1)*3+1)       # 1/4~3/4流入
        for i in 1:temp1
            for l in 1:nval
                Qbase[i,cellymax,l] = Qbase[i,cellymax-1,l]
            end
        end
        for i in temp1+1:temp2
            for l in 1:nval
                Qbase[i,cellymax,l] = bdcon[4][l+1]
            end
            T = bdcon[4][nval+2]
            p = (Qbase[i,cellymax,1]*Rd) * T
            Qbase[i,cellymax,4] = p
        end
        for i in temp2+1:cellxmax
            for l in 1:nval
                Qbase[i,cellymax,l] = Qbase[i,cellymax-1,l]
            end
        end
    elseif Int(bdcon[4][1]) == 12
        for i in 1:cellxmax
            for l in 1:nval
                Qbase[i,cellymax,l] = bdcon[4][l+1]
            end
            T = bdcon[4][nval+2]
            p = (Qbase[i,cellymax,1]*Rd) * T
            Qbase[i,cellymax,4] = p
        end
    elseif Int(bdcon[4][1]) == 9
        for i in 1:cellxmax
            for l in 1:nval
                Qbase[i,cellymax,l] = bdcon[4][l+1]
            end

            Ayx = vecAy[i,cellymax,1] / (vecAy[i,cellymax,1]^2 + vecAy[i,cellymax,2]^2)^0.5
            Ayy = vecAy[i,cellymax,2] / (vecAy[i,cellymax,1]^2 + vecAy[i,cellymax,2]^2)^0.5

            #=
            Un = Qbase[i,cellymax-1,2]*Ayx + Qbase[i,cellymax-1,2]*Ayy
            if Un < 0

            end
            =#
        end
    else
        println("------------------------")
        println(" boundary error ")
        println(" please check boundary number ")
        println("------------------------")
        throw(UndefVarError(:x))
    end

    return Qbase
end

# ------------------------------------
# check boundary condition
# ------------------------------------
function check_bd(bdcon)
    bdnum = [11, 12, 101, 102, 21, 31, 32, 33, 41]
    for l in 1:4
        di = 0
        for nn in bdnum
            if Int(bdcon[l][1]) == nn
                di = 1
            end
        end

        if di == 0
            println("\n check boundary condition ! \n")
            throw(UndefVarError(:x))
        end
    end
end

# ------------------------------------
# set boundary condition number for mpi
# ------------------------------------
function set_boundary_number(bdcon, rank, size)
    #=
    bd_NSEW = [N, S, E, W]
    北南東西に対応する境界条件を入れる．
    number <  0 : 境界条件,bdcon
    number >= 0 : 対応するセルのrank
    =#
    #=
    evaluation system

    parallel num = 4

      rank      div(rank, size^0.5)    rem(rank, size^0.5)

      2 | 3            1 | 1                  0 | 1
     -------          -------                -------
      0 | 1            0 | 0                  0 | 1

    parallel num = 16

          rank              div(rank, size^0.5)       rem(rank, size^0.5)

    12 | 13 | 14 | 15          3  | 3 | 3  | 3          0  | 1 | 2  | 3
    ------------------        ------------------       ------------------
    8  | 9  | 10 | 11          2  | 2 | 2  | 2          0  | 1 | 2  | 3
    ------------------        ------------------       ------------------
    4  | 5  | 6  | 7           1  | 1 | 1  | 1          0  | 1 | 2  | 3
    ------------------        ------------------       ------------------
    0  | 1  | 2  | 3           0  | 0 | 0  | 0          0  | 1 | 2  | 3

    0 : upper or left edge
    sqrt(parallel)-1 : lower or right edge

    =#
    bd_NSEW = zeros(Int, 4)
    y = div(rank, size^0.5)
    x = rem(rank, size^0.5)

    if y == 0 && y == Int(size^0.5 - 1)
        bd_NSEW[1] = - Int(bdcon[2][1])
        bd_NSEW[2] = - Int(bdcon[1][1])
    elseif y == 0
        bd_NSEW[1] = rank + size
        bd_NSEW[2] = - Int(bdcon[1][1])
    elseif y == Int(size^0.5 - 1)
        bd_NSEW[1] = - Int(bdcon[2][1])
        bd_NSEW[2] = rank - size
    else
        bd_NSEW[1] = rank + size
        bd_NSEW[2] = rank - size
    end

    if x == 0 && x == Int(size^0.5 - 1)
        bd_NSEW[3] = - Int(bdcon[4][1])
        bd_NSEW[4] = - Int(bdcon[3][1])
    elseif x == 0
        bd_NSEW[3] = rank + 1
        bd_NSEW[4] = - Int(bdcon[3][1])
    elseif x == Int(size^0.5 - 1)
        bd_NSEW[3] = - Int(bdcon[4][1])
        bd_NSEW[4] = rank - 1
    else
        bd_NSEW[3] = rank + 1
        bd_NSEW[4] = rank - 1
    end
    return bd_NSEW
end


text:
module StructJuMP

using JuMP # To reexport, should be using (not import)
import MathProgBase
import ReverseDiffSparse
using Printf, Libdl

# These modules could be optional.
import StructJuMPSolverInterface
import MPI

export StructuredModel, getStructure, getparent, getchildren, getProcIdxSet,
       num_scenarios, @second_stage, getprobability, getMyRank, mpiWrapper
       #BendersBridge, DLP, loadAndSolveConicProblem
# Macro to exportall
macro exportall(pkg)
    Expr(:export, names(JuMP)...)
end
@exportall JuMP

# ---------------
# StructureData
# ---------------
if @isdefined(MPI)
    mutable struct MPIWrapper
        comm::MPI.Comm
        init::Function

        function MPIWrapper()
            instance = new(MPI.Comm(-1))
            finalizer(freeMPIWrapper, instance)

            instance.init = function(ucomm::MPI.Comm)
                if @isdefined(MPI) && MPI.Initialized() && ucomm.val == -1
                    instance.comm = MPI.COMM_WORLD
                elseif @isdefined(MPI) && !MPI.Initialized()
                    MPI.Init()
                    instance.comm = MPI.COMM_WORLD
                elseif @isdefined(MPI) && MPI.Initialized() && ucomm.val != -1
                    instance.comm = ucomm
                elseif @isdefined(MPI) && MPI.Finalized()
                    error("MPI is already finalized!")
                else
                    #doing nothing
                end
            end

            return instance
        end

    end
    function freeMPIWrapper(instance::MPIWrapper)
        if @isdefined(MPI) && MPI.Initialized() && !MPI.Finalized()
            #MPI.Finalize()
            #println("Finalized would have been called")
        end
    end

    const mpiWrapper = MPIWrapper();

    mutable struct StructureData
        probability::Dict{Int,Float64}
        children::Dict{Int,JuMP.Model}
        parent
        num_scen::Int
        othermap::Dict{JuMP.Variable,JuMP.Variable}
        mpiWrapper::MPIWrapper
    end
else
    mutable struct StructureData
        probability::Dict{Int,Float64}
        children::Dict{Int,JuMP.Model}
        parent
        num_scen::Int
        othermap::Dict{JuMP.Variable,JuMP.Variable}
    end
end
default_probability(m::JuMP.Model) = 1 / num_scenarios(m)
default_probability(::Nothing) = 1.0

# ---------------
# StructuredModel
# ---------------

function structprinthook(io::IO, m::Model)
    print(io, m, ignore_print_hook=true)
    print(io, "*** children ***\n")
    # TODO it would be nice to indent all the children
    # essentially wrap the IO object (subclass it) to add 4 spaces before each line
    # this would then recursively be wrapped as more stages are added
    for (id, mm) in getchildren(m)
      @printf(io, "Child ID %d:\n", id)
      print(io, mm)
      print(io, "\n")
    end
end

# Constructor with the number of scenarios
function StructuredModel(;solver=JuMP.UnsetSolver(), parent=nothing, same_children_as=nothing, id=0, comm=@isdefined(MPI) ? MPI.Comm(-1) : -1, num_scenarios::Int=0, prob::Float64=default_probability(parent))
    m = JuMP.Model(solver=solver)
    if parent === nothing
        id = 0
        if @isdefined(MPI)
            mpiWrapper.init(comm)
        end
        if @isdefined(StructJuMPSolverInterface)
            JuMP.setsolvehook(m,StructJuMPSolverInterface.sj_solve)
        end
    else
        @assert id != 0
        stoch = getStructure(parent)
        stoch.children[id] = m
        stoch.probability[id] = prob
    end

    if same_children_as !== nothing
      if !isa(same_children_as, JuMP.Model) || !haskey(same_children_as.ext, :Stochastic)
        error("The JuMP model given for the argument `same_children_as' is not valid. Please create it using the `StructuredModel' function.")
      end
      probability = same_children_as.ext[:Stochastic].probability
      children = same_children_as.ext[:Stochastic].children
    else
      probability = Dict{Int, Float64}()
      children = Dict{Int, JuMP.Model}()
    end
    if @isdefined(MPI)
        m.ext[:Stochastic] = StructureData(probability, children, parent, num_scenarios, Dict{JuMP.Variable,JuMP.Variable}(), mpiWrapper)
    else
        m.ext[:Stochastic] = StructureData(probability, children, parent, num_scenarios, Dict{JuMP.Variable,JuMP.Variable}())
    end

    # Printing children is important as well
    JuMP.setprinthook(m, structprinthook)

    m
end

# -------------
# Get functions
# -------------

getStructure(m::JuMP.Model)   = m.ext[:Stochastic]
getparent(m::JuMP.Model)      = getStructure(m).parent
getchildren(m::JuMP.Model)    = getStructure(m).children::Dict{Int,JuMP.Model}
getprobability(m::JuMP.Model) = getStructure(m).probability::Dict{Int, Float64}
num_scenarios(m::JuMP.Model)  = getStructure(m).num_scen::Int

function getMyRank()
    myrank = 0;
    mysize = 1;
    if @isdefined(MPI) && !MPI.Initialized()
        MPI.Init()
    end
    if @isdefined(MPI) && MPI.Initialized() && !MPI.Finalized()
        comm = MPI.COMM_WORLD
        mysize = MPI.Comm_size(comm)
        myrank = MPI.Comm_rank(comm)
    end
    return myrank,mysize
end

function getProcIdxSet(numScens::Integer)
    mysize = 1;
    myrank = 0;
    if @isdefined(MPI) == true && MPI.Initialized() == true
        comm = MPI.COMM_WORLD
        mysize = MPI.Comm_size(comm)
        myrank = MPI.Comm_rank(comm)
    end
    # Why don't we just take a round-and-robin?
    proc_idx_set = Int[];
    for s = myrank:mysize:(numScens-1)
        push!(proc_idx_set, s+1);
    end
    return proc_idx_set;
end

function getProcIdxSet(m::JuMP.Model)
    numScens = num_scenarios(m)
    return getProcIdxSet(numScens)
end

macro second_stage(m,ind,code)
    return quote
        proc_idx_set = getProcIdxSet($(esc(m)))
        for $(esc(ind)) in proc_idx_set
            $(esc(code))
        end
    end
end

include("nlp.jl")

#include("BendersBridge.jl")

end


text:
"""
    foldxd(step, xform::Transducer, array; [init, simd, basesize, threads_basesize, pool])

e**X**tended **d**istributed fold (reduce).  This is a distributed
`reduce` based on extended fold protocol defined in Transducers.jl.

Unlike [`foldxt`](@ref), early termination by [`reduced`](@ref) is not
supported yet.

Use [`dcollect`](@ref) or [`dcopy`](@ref) to collect results into a
container.

See also: [Parallel processing tutorial](@ref tutorial-parallel),
[`foldxl`](@ref), [`foldxt`](@ref).

!!! compat "Transducers.jl 0.4.3"

    New in version 0.4.3.

# Keyword Arguments
- `pool::AbstractWorkerPool`: Passed to `Distributed.remotecall`.
- `basesize::Integer = amount(array) ÷ nworkers()`: A size of chunk in
  `array` that is processed by each worker.  A smaller size may be
  required when computation time for processing each item can
  fluctuate a lot.
- `threads_basesize::Integer = basesize ÷ nthreads()`: A size of chunk
  in `array` that is processed by each task in each worker process.
  The default setting assumes that the number of threads used in all
  workers are the same.  For heterogeneous setup where each worker
  process has different number of threads, it may be required to use
  smaller `threads_basesize` _and_ `basesize` to get a good performance.
- For other keyword arguments, see [`foldl`](@ref).

# Examples
```jldoctest
julia> using Transducers

julia> foldxd(+, 1:3 |> Map(exp) |> Map(log))
6.0

julia> 1:5 |> Filter(isodd) |> foldxd(+)
9

julia> foldxd(TeeRF(min, max), [5, 2, 6, 8, 3])
(2, 8)
```
"""
foldxd
foldxd(rf; kw...) = itr -> foldxd(rf, itr; kw...)
foldxd(step, xform::Transducer, itr; init = DefaultInit, kwargs...) =
    unreduced(dtransduce(xform, Completing(step), init, itr; kwargs...))
foldxd(step::F, foldable; kwargs...) where {F} =
    foldxd(step, extract_transducer(foldable)...; kwargs...)

"""
    dtransduce(xform::Transducer, step, init, array; [simd, basesize, threads_basesize, pool])

See [`foldxd`](@ref) and [`transduce`](@ref).
"""
function dtransduce(
    xform::Transducer,
    step,
    init,
    coll0;
    simd::SIMDFlag = Val(false),
    basesize::Union{Integer,Nothing} = nothing,
    threads_basesize::Union{Integer,Nothing} = nothing,
    pool::Distributed.AbstractWorkerPool = Distributed.default_worker_pool(),
    _remote_reduce = _transduce_assoc_nocomplete,
)
    rf0 = _reducingfunction(xform, step; init = init)
    rf1, coll = retransform(rf0, coll0)
    rf = maybe_usesimd(rf1, simd)
    isempty(coll) && return check_empty(rf, complete(rf, start(rf, init)))
    load_me_everywhere()
    basesize = if basesize === nothing
        max(1, amount(coll) ÷ Distributed.nworkers())
    else
        Int(basesize)
    end
    threads_basesize = if threads_basesize === nothing
        max(1, basesize ÷ Threads.nthreads())
    else
        Int(threads_basesize)
    end
    @argcheck basesize > 0
    @argcheck threads_basesize > 0
    futures = map(firstindex(coll):basesize:lastindex(coll)) do start
        Distributed.remotecall(
            _remote_reduce,
            pool,
            rf,
            init,
            coll[start:min(end, start - 1 + basesize)],
            threads_basesize,
        )
    end
    # TODO: Cancel remote computation when there is a Reduced.
    result = Distributed.remotecall_fetch(pool, futures) do futures
        results = map(fetch, futures)
        return complete(rf, combine_all(rf, results))
        # TODO: call combine asynchronously
    end
    return check_empty(rf, result)
end
# Note: Calling combine on remote is useful not only for computation
# off-loading but also for making it actually work. If `rf` contains any
# stateful transducers that use any anonymous functions defined in `Main`,
# `unwrap` will detect that the private states as incompatible. So, we need to
# use `rf` that is deserialized by Distributed, when combining and completing
# the stateful reductions.

function load_me_everywhere()
    pkgid = Base.PkgId(@__MODULE__)
    @everywhere Base.require($pkgid)
end

"""
    dcopy(xf::Transducer, T, reducible; [basesize, threads_basesize]) :: Union{T, Empty{T}}
    dcopy(xf::Transducer, reducible::T; [basesize, threads_basesize]) :: Union{T, Empty{T}}
    dcopy([T,] itr; [basesize, threads_basesize]) :: Union{T, Empty{T}}

Distributed.jl-based parallel version of [`copy`](@ref).  Keyword
arguments are passed to [`foldxd`](@ref).  For examples, see
[`tcopy`](@ref).

See also: [Parallel processing tutorial](@ref tutorial-parallel)
(especially [Example: parallel `collect`](@ref tutorial-parallel-collect)).

!!! compat "Transducers.jl 0.4.5"

    New in version 0.4.5.

!!! compat "Transducers.jl 0.4.8"

    `dcopy` now accepts iterator comprehensions and eductions.
"""
dcopy(xf, T, reducible; kwargs...) =
    foldxd(append!!, Map(SingletonVector) ∘ xf, reducible; init = Empty(T), kwargs...)
dcopy(xf, reducible; kwargs...) = dcopy(xf, _materializer(reducible), reducible; kwargs...)

function dcopy(::Type{T}, itr; kwargs...) where {T}
    xf, foldable = extract_transducer(itr)
    return dcopy(xf, T, foldable; kwargs...)
end

function dcopy(itr; kwargs...)
    xf, foldable = extract_transducer(itr)
    return dcopy(xf, foldable; kwargs...)
end

"""
    dcollect(xf::Transducer, reducible; [basesize, threads_basesize]) :: Union{Vector, Empty{Vector}}
    dcollect(itr; [basesize, threads_basesize]) :: Union{Vector, Empty{Vector}}

Distributed.jl-based parallel version of [`collect`](@ref).
This is just a short-hand notation of `dcopy(xf, Vector, reducible)`.
Use [`dcopy`](@ref) to get a container other than a `Vector`.

See also: [Parallel processing tutorial](@ref tutorial-parallel)
(especially [Example: parallel `collect`](@ref tutorial-parallel-collect)).

!!! compat "Transducers.jl 0.4.5"

    New in version 0.4.5.

!!! compat "Transducers.jl 0.4.8"

    `dcollect` now accepts iterator comprehensions and eductions.
"""
dcollect(xf, reducible; kwargs...) = dcopy(xf, Vector, reducible; kwargs...)
dcollect(itr; kwargs...) = dcollect(extract_transducer(itr)...; kwargs...)


text:
using Photometry.Aperture: bounds

@testset "Apertures" begin
    ap_rect = RectangularAperture(50, 40, 10, 10, 0)
    @test center(ap_rect) == (50, 40)
    @test bounds(ap_rect) == (45, 55, 35, 45)
    @test size(ap_rect) == (11, 11)
    @test size(ap_rect, 1) == 11
    @test RectangularAperture([50, 40], 10, 10, 0) == ap_rect

    ap_ann = RectangularAnnulus(50, 40, 5, 10, 10, 0)
    @test bounds(ap_ann) == (45, 55, 35, 45)
    @test size(ap_ann) == (11, 11)
    @test RectangularAnnulus([50, 40], 5, 10, 10, 0) == ap_ann
end

@testset "Rectangle Aperture" begin
    ap0 = RectangularAperture(0, 0, 0, 0, 0)

    @test sprint(show, ap0) == "RectangularAperture(0, 0, w=0, h=0, θ=0°)"

    ap1 = RectangularAperture(0, 0, 1, 1, 0)

    @test sprint(show, ap1) == "RectangularAperture(0, 0, w=1, h=1, θ=0°)"
end

@testset "Rectangle Annulus" begin

    ap1 = RectangularAnnulus(0, 0, 1, 1, 1, 0)
    @test center(ap1) == (0, 0)
    @test sprint(show, ap1) == "RectangularAnnulus(0.0, 0.0, w_in=1.0, w_out=1.0, h_in=1.0, h_out=1.0, θ=0.0°)"
end


text:
X = [+1 -1 -1 -1 -1 +1;
         -1 +1 -1 -1 +1 -1;
         -1 -1 +1 +1 -1 -1;
         -1 -1 +1 +1 -1 -1;
         -1 -1 +1 +1 -1 -1;
         -1 +1 -1 -1 +1 -1;
         +1 -1 -1 -1 -1 +1.]

O = [+1 +1 +1 +1 +1 +1;
         +1 -1 -1 -1 -1 +1;
         +1 -1 -1 -1 -1 +1;
         +1 -1 -1 -1 -1 +1;
         +1 -1 -1 -1 -1 +1;
         +1 -1 -1 -1 -1 +1;
         +1 +1 +1 +1 +1 +1.]

F1 = [-1 -1 -1 -1 -1 -1;
          -1 +1 +1 +1 +1 -1;
          -1 +1 -1 -1 -1 -1;
          -1 +1 +1 +1 +1 -1;
          -1 +1 -1 -1 -1 -1;
          -1 +1 -1 -1 -1 -1;
          -1 -1 -1 -1 -1 -1.]

F2 = [+1 +1 +1 +1 +1 -1;
          +1 -1 -1 -1 -1 -1;
          +1 -1 -1 -1 -1 -1;
          +1 +1 +1 +1 +1 -1;
          +1 -1 -1 -1 -1 -1;
          +1 -1 -1 -1 -1 -1;
          +1 -1 -1 -1 -1 -1.]


text:
using SimpleNavierStokes, Test

@testset "Lid-Driven cavity runs" begin
    @test LidDrivenCavity(tfinal=1).tfinal == 1
    @test LidDrivenCavity(tfinal=2).tfinal == 2
end


text:
D=DBayesianNode(:D, ["d1","d2","d_water"])

a1=CBayesianNode(:a1)
a2=CBayesianNode(:a2)
a3=CBayesianNode(:a3)

bn=BayesianNetwork([D,a1,a2])

add_edge!(bn,D,a1)
add_edge!(bn,D,a2)

analysis=[a1,a2]

# Distributions

using Distributions

pdfNodeDict = Dict{Symbol, Array{Function,1}}()

for a in analysis
    a_d1_dist = Normal(rand(1:200), rand(1:15))
    a_d1_func = x -> pdf(a_d1_dist,x)
    a_d2_dist = Normal(rand(200:400), rand(1:15))
    a_d2_func = x -> pdf(a_d2_dist,x)

    samples = [ rand(dist) for i=1:100, dist in [a_d1_dist, a_d2_dist] ]

    a_d_water_pdf = x -> 1/(maximum(samples) - minimum(samples))

    pdfNodeDict[a.label] = [a_d1_func, a_d2_func, a_d_water_pdf]

    add_probability!(bn, P(a.label|:D), ProbabilityDensityDistribution(states(D), [a_d1_func, a_d2_func, a_d_water_pdf]))
end

l_states_d = length(states(D))
add_probability!(bn, P(:D), ProbabilityDistribution(states(D), fill(1/l_states_d, l_states_d)))
inputDict = Dict{Symbol,Any}({:a1 => 30, :a2 => 30})
testCPD = P(:D|inputDict)
diagProbs = query(bn,testCPD)

# Normalize results of query so we can compare
normalized_diagProbs = normalize_factor(convert(Array{Float64,1},collect(values(diagProbs))))

# For all a in analysis, run every function of a and give it the input defined in inputDict. The reducing steps performs an element-wise multiplication, and finally it's normalized.
calculatedDiagProbs = normalize_factor(reduce(.*, [ map(x -> x(inputDict[a.label]), pdfNodeDict[a.label]) for a in analysis ]))

# Test if the two normalized arrays are approximately equal
@test sum(map(x -> isapprox(normalized_diagProbs[x[1]], calculatedDiagProbs[x[1]]), enumerate(calculatedDiagProbs))) == l_states_d

#bn.query(P(:D|{:a1 => 30, :a2 => 250, a_d_water_pdf => 0}))


text:
lll(H::AbstractArray{Rational{Td},2},δ=.75,η=.51) where {Td<:Integer} =
    l2(H,Rational{Td},δ,η)

"""
    B = l2(H::AbstractArray{Td,2},TG::Type{Tg},δ=.75,η=.51) where {Td<:Number,Tg<:Number}

Do L2 lattice reduction of matrix `H` using optional parameters `δ` and `η`,
using a Gram matrix of type TG. The output is `B`, an LLL-reduced basis such
that `H*T = B`, where `T` is a unimodular matrix which is not returned.

Follows L2 from Damien Stehle, "Floating-Point LLL: Theoretical and
Practical Aspects", chapter 5 of "The LLL Algorithm, survey and
applications", Springer-Verlag, 2009.

The `l2` function is still in a raw (alpha) state and may change or be removed.
The `lll` function in this library also does LLL reduction, and functions on
a wider set of input types, such as complex data. When `l2` works, it is
generally faster than `lll` on small bases, say of dimensions less than 80.

# Examples
```jldoctest
julia> H= [1 2; 3 4];B,_ = LLLplus.l2(H); B
┌ Warning: l2 is in a raw (alpha) state and may change. See the help text.
└ @ LLLplus src/l2.jl:45
2×2 Matrix{Int64}:
 1  -1
 1   1

julia> H= [.5 2; 3 4]; B,_= LLLplus.l2(H); B
┌ Warning: l2 is in a raw (alpha) state and may change. See the help text.
└ @ LLLplus src/l2.jl:45
2×2 Matrix{Float64}:
 1.5  -1.0
 1.0   2.0

julia> N=30;H = randn(N,N); B,T = LLLplus.l2(H);
┌ Warning: l2 is in a raw (alpha) state and may change. See the help text.
└ @ LLLplus src/l2.jl:45

```
"""
function l2(H::AbstractArray{Td,2},TG::Type{Tg}=Td,δ=.75,η=.51) where
    {Td<:Any,Tg<:Any}

    @warn "l2 is in a raw (alpha) state and may change. See the help text." maxlog=1 _file="src/l2.jl"

    if !(0.25 < δ < 1.0)
        error("δ must be between 1/4 and 1.");
    end
    if !(0.5 < η < sqrt(δ))
        error("η must be between 1/2 and sqrt(δ).");
    end
    n,d = size(H)
    Ti= LLLplus.getIntType(Td)
    Tf = float(Td)

    # Follows Stehle
    ϵ = .001  # ϵ = eps(Td) may be too small
    C = ϵ
    if Tf <: Complex
        @error "`l2` does not handle complex data; try `lll`."
        return
    end
    Tfe = real(Tf)
    l = precision(Tfe)
    left = d^2*( ((1+η)^2+ϵ)/(δ-ϵ^2) )^d * 2^(-l+10+C*d)
    right = minimum([ϵ η-.5 1-δ])
    # if left>right
    #     @warn "Precision of l=$l bits for $Tf is not sufficient "*
    #           "to satisfy Stehle's requirement [survey book, ch 5, thm 2]."
    # end
    if l<1.6*d
        @warn "Precision of l=$l bits does not even satisfy Stele's "*
        "asymptotic requirement of l≥1.6d."
    end

    δ= Tfe(δ);
    η= Tfe(η)
    B = copy(H)
    G = Tg.(B)'*B
    ηb= (η+.5)/2
    δb= (δ+1)/2

    μ = Matrix{Tf}(I,n,d)
    r = zeros(Tf,n,d)
    s = Vector{Tf}(undef,d)
    X = Vector{Ti}(undef,d)

    r[1,1] = G[1,1]

    κ=2;
    while κ<=d
        lazysizereduction!(ηb,κ,B,G,r,μ,s,X,n,d,Tg)

        κp=κ
        while κ>=2 && δb*real(r[κ-1,κ-1]) > real(s[κ-1])
            κ-=1
        end
        for i=1:κ-1
            μ[κ,i] = μ[κp,i]
            r[κ,i] = r[κp,i]
        end
        r[κ,κ] = s[κ]

        if κ!=κp
            bκp = B[:,κp]
            copyto!(B,n*κ+1,B,n*(κ-1)+1,n*length(κ:κp-1))
            B[:,κ] .= bκp

            if κ<d
                for i=κ:κp
                    @inbounds for d1=1:d
                        G[d1,i] = B[1,d1]*B[1,i]
                        for d2=2:d
                            G[d1,i] += B[d2,d1]*B[d2,i]
                        end
                    end
                end
                # maybe make G a `Hermitian` view of an underlying matrix,
                # set values using G.data[...] above, and skip the loop
                # below? see also discourse: https://bit.ly/36QguCG

                for i=κ:κp
                    for dx=1:d
                        G[i,dx] = G[dx,i]
                    end
                end
            end
        end
        κ+=1
    end
    T = Ti.(round.(inv(H)*B))
    return B,T
end

function lazysizereduction!(ηb,κ,B,G,r,μ,s,X,n,d,Tg)

    @label startCholesky
    # update the cholesky factorization
    i = κ
    for j = 1:i-1
        r[i,j] = G[i,j]
        for k = 1:j-1
            r[i,j] -= μ[j,k]*r[i,k]
        end
        μ[i,j] = r[i,j]/r[j,j]
    end
    s[1] = G[i,i]
    for j = 2:i
        s[j] = s[j-1]-μ[i,j-1]*r[i,j-1]
    end
    r[i,i] = s[i]

    μκ = view(μ,κ,1:κ-1)
    if maximum(abs.(μκ))≤ηb
        return
    else
        @inbounds for i = κ-1:-1:1
            X[i] = round(μ[κ,i])
            for j=1:i-1
                μ[κ,j]-=X[i]*μ[i,j]
            end
        end
        @inbounds for i = κ-1:-1:1
            for nx=1:n
                B[nx,κ] -=X[i]*B[nx,i]
            end
        end
        @inbounds for d1=1:d
            G[d1,κ] = B[1,d1]*B[1,κ]
            for d2=2:d
                G[d1,κ] += B[d2,d1]*B[d2,κ]
            end
        end
         for dx=1:d
            G[κ,dx] = G[dx,κ]
        end
    end
    @goto startCholesky
end
#=
"""
    B = l2avx(H::AbstractArray{Td,2},TG::Type{Tg},δ=.75,η=.51) where
                         {Td<:Number,Tg<:Number}

A version of the `l2` function with a few calls to the `@avx` macro
from the `LoopVectorization.jl` package.  See the `l2` help text

# Examples
```julia
julia> using LLLplus
julia> using LoopVectorization
julia> H= [1 2; 3 4];B = l2avx(H)
┌ Warning: l2avx is in a raw (alpha) state and may change. See the help text.
└ @ LLLplus ~/shared/LLLplus/src/l2.jl:42
2×2 Matrix{Int64}:
 1  -1
 1   1
```
"""
function l2avx(H::AbstractArray{Td,2},TG::Type{Tg}=Td,δ=.75,η=.51) where
    {Td<:Number,Tg<:Number}

    @warn "l2avx is in a raw (alpha) state and may change. See the help text." maxlog=1

    if !(0.25 < δ < 1.0)
        error("δ must be between 1/4 and 1.");
    end
    if !(0.5 < η < sqrt(δ))
        error("η must be between 1/2 and sqrt(δ).");
    end
    n,d = size(H)
    Ti= LLLplus.getIntType(Td)

    Tf = float(Td)

    # Follows Stehle
    ϵ = .001  # ϵ = eps(Td) may be too small
    C = ϵ
    if Tf <: Complex
        @error "`l2avx` does not handle complex data; try `lll`."
        return
    end
    Tfe = real(Tf)
    l = precision(Tfe)
    left = d^2*( ((1+η)^2+ϵ)/(δ-ϵ^2) )^d * 2^(-l+10+C*d)
    right = minimum([ϵ η-.5 1-δ])
    # if left>right
    #     @warn "Precision of l=$l bits for $Tf is not sufficient "*
    #           "to satisfy Stehle's requirement [survey book, ch 5, thm 2]."
    # end
    if l<1.6*d
        @warn "Precision of l=$l bits does not even satisfy Stele's "*
              "asymptotic requirement of l≥1.6d."
    end

    δ= Tfe(δ);
    η= Tfe(η)
    B = copy(H)
    G = Tg.(B)'*B
    ηb= (η+.5)/2
    δb= (δ+1)/2

    μ = Matrix{Tf}(I,n,d)
    r = zeros(Tf,n,d)
    s = Vector{Tf}(undef,d)
    X = Vector{Ti}(undef,d)

    r[1,1] = G[1,1]

    κ=2;
    while κ<=d
        lazysizereductionAVX!(ηb,κ,B,G,r,μ,s,X,n,d,Tg)

        κp=κ
        while κ>=2 && δb*real(r[κ-1,κ-1]) > real(s[κ-1])
            κ-=1
        end
        for i=1:κ-1
            μ[κ,i] = μ[κp,i]
            r[κ,i] = r[κp,i]
        end
        r[κ,κ] = s[κ]

        if κ!=κp
            bκp = B[:,κp]
            copyto!(B,n*κ+1,B,n*(κ-1)+1,n*length(κ:κp-1))
            B[:,κ] .= bκp

            if κ<d
                @avx for i=κ:κp
                    for d1=1:d
                        G[d1,i] = B[1,d1]*B[1,i]
                        for d2=2:d
                            G[d1,i] += B[d2,d1]*B[d2,i]
                        end
                    end
                end
                # maybe make G a `Hermitian` view of an underlying matrix,
                # set values using G.data[...] above, and skip the loop
                # below? see also discourse: https://bit.ly/36QguCG

                @inbounds for i=κ:κp
                    for dx=1:d
                        G[i,dx] = G[dx,i]
                    end
                end
            end
        end
        κ+=1
    end
    return B
end

function lazysizereductionAVX!(ηb,κ,B,G,r,μ,s,X,n,d,Tg)

    @label startCholesky
    # update the cholesky factorization
    i = κ
    @inbounds for j = 1:i-1
        r[i,j] = G[i,j]
        for k = 1:j-1
            r[i,j] -= μ[j,k]*r[i,k]
        end
        μ[i,j] = r[i,j]/r[j,j]
    end
    s[1] = G[i,i]
    @inbounds for j = 2:i
        s[j] = s[j-1]-μ[i,j-1]*r[i,j-1]
    end
    r[i,i] = s[i]

    μκ = view(μ,κ,1:κ-1)
    if maximum(abs.(μκ))≤ηb
        return
    else
        @inbounds for i = κ-1:-1:1
            X[i] = round(μ[κ,i])
            for j=1:i-1
                μ[κ,j]-=X[i]*μ[i,j]
            end
        end
        @avx for nx=1:n
            for i = 1:κ-1
                B[nx,κ] -=X[i]*B[nx,i]
            end
        end
        @avx for d1=1:d
            G[d1,κ] = B[1,d1]*B[1,κ]
            for d2=2:d
                G[d1,κ] += B[d2,d1]*B[d2,κ]
            end
        end
        for dx=1:d
            G[κ,dx] = G[dx,κ]
        end
    end
    @goto startCholesky
end
=#

"""
    dataTypeForGram(nbits,d)

Given a bit depth `nbits`, dimension `d` of a integer square matrix, return
the integer datatype needed to represent the Gram matrix in full
precision. This is conservative for most matrices, yet in some corner cases,
the deeper bit depth in the datatype returned is needed.

"""
function dataTypeForGram(nbits,d)
    # nbits*2-1 bits required to represent output of each multiply
    # additional log2(d) bits required for sums. So for each element of the
    # Gram matrix we need `bitsRequired` bits.
    bitsRequired = Int64(ceil(nbits*2-1+log2(d)))
    return intTypeGivenBitsRequired(bitsRequired)
end

"""
    dataTypeForGram(B::AbstractArray{Td,2}) where {Td<:Number}

Given a matrix B of eltype Td, return a datatype for representing the Gram
matrix of B. When Td<:Integer, we calculate the Gram matrix as BigInt, then
use its maximum value to get an `nbits` that to calculate a Gram data type
from `dataTypeForGram(nbits,d)`. Yes, it may be overkill to calculate the
Gram to find out the bit precision needed to represent it.

"""
function dataTypeForGram(B::AbstractArray{Td,2}) where {Td<:Number}
    if Td<:AbstractFloat
        return Float64  # maybe return Td?
    elseif Td<:Integer
        G = big.(B)'*B
        bitsRequired = Int64(ceil(log2(maximum(abs.(G)))))
        return intTypeGivenBitsRequired(bitsRequired)
    #elseif Complex
        #?
    else
        return Td
    end
end

"""
    intTypeGivenBitsRequired(bitsRequired)

Given the number of `bitsRequired`, return the smallest Integer type that
can represent that many bits, as well as a sign bit.

"""
function intTypeGivenBitsRequired(bitsRequired)
#    intTypes = [Int8, Int16, Int32, Int64,Int128,Int256,Int512,Int1024]
    intTypes = [Int8, Int16, Int32, Int64,Int128]
    intBits = sizeof.(intTypes)*8 .-1 # 1 bit for sign. could add 1 for spare
    ix = findfirst(bitsRequired .< intBits)
    return !isnothing(ix) ? intTypes[ix] : BigInt
end


text:
using Base.Threads

function showResult(dtStr::String ,AlarmDynStatus::Array{Bool,1}, AlarmStaticStatus::Array{Bool,1},CumSumDataobj::spaceCumSumData )
    segid1= CumSumDataobj.startID
    segid2= CumSumDataobj.endID
    if true #segid1 > 1
        println( "$dtStr: Dyn, ")
        for a in AlarmDynStatus[1:end ]
            a == true ? print( "Y,") : print(",")
        end
        println()
        println( "$dtStr: Sta, ")
        for a in AlarmStaticStatus[1:end ]
            a == true ? print( "Y,") : print(",")
        end
        println()
    end
end
function AlarmCreditCalAna( CumSumDataobj::spaceCumSumData)
    temobj = CumSumDataobj
    rectsumObj =temobj.rectSum

    print("+$(temobj.dt[1:8]):$(temobj.dt[9:14]) , ")

    e =  temobj.AlrmProcEndtID-temobj.AlrmProcStartID +1

    a= temobj.AlrmProcStartID - temobj.startID  + 1
    b= temobj.AlrmProcEndtID  - temobj.startID  + 1
    viewa  = @view rectsumObj[ a: b]
    # doesWinding( CumSumDataobj)
    # println("$a   $b   $(size(viewa ))")
    for k = 1 : e
        # print( rectsumObj[ k ] ,"," )
        if viewa[ k ] >= threshold
                print( viewa[ k ] ,"," )
        else
            　print(",")
        end
    end
    # println("-------")
end
include("../../../Jl-Aux/glbCfg/internalPara.jl")
function postProcFeatureData(dataObject::spaceCumSumData ,data::Array{UInt16,1} )
    # println("   ------------------       <?>            0     -------------         ",length(data ))
    cumSumProc( dataObject,  data  )
    # println("   ------------------       <?>        1         -------------         ",length(data ))

    vsc = glbvsc
    vTcNum = glbvTcNum

    # @info vsc, vTcNum
    arng = getSpaceLoopRng( dataObject )
    # println( "---------------",arng )
    # println("   ------------------       <?>        2         -------------         ",arng)
   @threads for vpos = arng
        # println("   ------------------       <?>        3         -------------         ",arng)
        calRectSum( dataObject, vpos ,vsc ,vTcNum   )
        # println("   ------------------       <?>        3         -------------         ",arng)
        CalDynmicTh( dataObject,  vpos )
    end
    # println("   ------------------       <?>        4         -------------         ",arng)
    AlarmRTCalAna(dataObject)
end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
function solveValueFunction()

#This function obtains the value function for each time period and
#the policy function (i.e. optimal next-period asset choice) for each time
#period. From there we can work the optimal consumption level.

#The approach taken is by backwards recursion. The optimisation each period
#is carried out using 'fminbnd'. This is an in-built optimiser in Matlab.
#The optimisation routine it uses is known as the 'golden search method'

#global T r tol minCons
#global numPointsA numPointsY Agrid Ygrid incTransitionMrx
#global Agrid1 EV1

#% ------------------------------------------------------------------------
# GENERATE MATRICES TO STORE NUMERICAL APPROXIMATIONS AND INITIATE AS NAN

# Matrices to hold the policy, value and marginal utility functions
V        = Array(Float64,T+1, numPointsA, numPointsY);
policyA1 = Array(Float64,T,   numPointsA, numPointsY);
policyC  = Array(Float64,T,   numPointsA, numPointsY);

#Matrices to hold expected value and marginal utility functions
EV  = Array(Float64,T+1, numPointsA, numPointsY);

#% ------------------------------------------------------------------------
#Set the terminal value function and expected value function to 0

EV[T + 1, :,:]  = .0;          # continuation value at T-1
V[T + 1,:,:]    = .0;
#% ------------------------------------------------------------------------
# SOLVE RECURSIVELY THE CONSUMER'S PROBLEM, STARTING AT TIME T-1 AND MOVING
# BACKWARDS TO ZERO, ONE PERIOD AT A TIME

for ixt=T:-1:1                                # Loop from time T-1 to 1
    Agrid1 = vec(Agrid[ixt + 1, :]);          # The grid on assets tomorrow, vector

    for ixA = 1:1:numPointsA                  # points on asset grid

        # STEP 1. solve problem at grid points in assets and income
        # ---------------------------------------------------------
        A    = Agrid[ixt, ixA];            # assets today
        lbA1 = Agrid[ixt + 1, 1];          # lower bound: assets tomorrow

        for ixY = 1:1:numPointsY               # points on income grid

            # Value of income and information for optimisation
            Y    = Ygrid[ixt, ixY];            # income today
            ubA1 = (A + Y - minCons)*(1+r);    # upper bound: assets tomorrow
            EV1  = vec(EV[ixt + 1,:, ixY]);    # relevant section of EV matrix (in assets tomorrow), vector

            intfV  = InterpIrregular(Agrid1,EV1, BCnil, InterpLinear); # Define the interpolation

            # Compute solution
            if (ubA1 - lbA1 < minCons)                       # if liquidity constrained
                negV = objectivefunc(lbA1, A, Y,Agrid1,EV1,intfV);
                policyA1[ixt,ixA,ixY] = lbA1;
            else                                             # if interior solution

                function obj(A1)
                    return objectivefunc(A1, A, Y, Agrid1,EV1,intfV)
                end

                Res = optimize(obj,lbA1,ubA1);          #Find the solution...

                policyA1[ixt,ixA,ixY]=Res.minimum;
                negV=Res.f_minimum;

                if isnan(negV)
                    println("Sendo error en ixt=$ixt, ixA=$ixA, ixY=$ixY")
                end

            end # if (ubA1 - lbA1 < minCons)

            # Store solution and its value
            policyC[ixt, ixA, ixY] = A + Y - policyA1[ixt, ixA, ixY]/(1+r);
            V[ixt, ixA, ixY]       = -negV;
        end #ixY

        # STEP 2. integrate out income today conditional on income
        # yesterday to get EV and EdU
        # --------------------------------------------------------
        realisedV = squeeze(V[ixt, ixA, :],1);
        for ixY = 1:1:numPointsY
            d=incTransitionMrx[ixY,:]*realisedV'
            EV[ixt, ixA, ixY]  = d[1,1];
        end #ixY
    end #ixA

    println("Passed period $ixt of $T")

end #ixt

return policyA1, policyC, V, EV

end #function

max_stars_repo_name: floswald/ucl-econ-julia

max_issues_repo_name: floswald/ucl-econ-julia

max_forks_repo_name: floswald/ucl-econ-julia


text:
using Flux
println("NNPDE_tests")
using DiffEqFlux
println("Starting Soon!")
using ModelingToolkit
using DiffEqBase
using Test, NeuralPDE
using GalacticOptim
using Optim
using CUDA
using Quadrature
using QuasiMonteCarlo
import ModelingToolkit: Interval, infimum, supremum

using Random
Random.seed!(100)

cb = function (p,l)
    println("Current loss is: $l")
    return false
end
CUDA.allowscalar(false)
#const gpuones = cu(ones(1))

## ODE
println("ode")
@parameters θ
@variables u(..)
Dθ = Differential(θ)

# 1D ODE
eq = Dθ(u(θ)) ~ θ^3 + 2.f0*θ + (θ^2)*((1.f0+3*(θ^2))/(1.f0+θ+(θ^3))) - u(θ)*(θ + ((1.f0+3.f0*(θ^2))/(1.f0+θ+θ^3)))

# Initial and boundary conditions
bcs = [u(0.) ~ 1.0f0]

# Space and time domains
domains = [θ ∈ Interval(0.0,1.0)]
# Discretization
dt = 0.1f0
# Neural network
inner = 20
chain = Chain(Dense(1,inner,Flux.σ),
              Dense(inner,inner,Flux.σ),
              Dense(inner,inner,Flux.σ),
              Dense(inner,inner,Flux.σ),
              Dense(inner,inner,Flux.σ),
              Dense(inner,1)) |> gpu

initθ = DiffEqFlux.initial_params(chain) |> gpu

strategy = NeuralPDE.GridTraining(dt)
discretization = NeuralPDE.PhysicsInformedNN(chain,
                                             strategy;
                                             init_params = initθ
                                             )

@named pde_system = PDESystem(eq,bcs,domains,[θ],[u(θ)])
prob = NeuralPDE.discretize(pde_system,discretization)
symprob = NeuralPDE.symbolic_discretize(pde_system,discretization)
res = GalacticOptim.solve(prob, ADAM(1e-2); maxiters=2000)
phi = discretization.phi

analytic_sol_func(t) = exp(-(t^2)/2)/(1+t+t^3) + t^2
ts = [infimum(d.domain):dt/10:supremum(d.domain) for d in domains][1]
u_real  = [analytic_sol_func(t) for t in ts]
u_predict  = [first(Array(phi([t],res.minimizer))) for t in ts]

@test u_predict ≈ u_real atol = 0.2

# t_plot = collect(ts)
# plot(t_plot ,u_real)
# plot!(t_plot ,u_predict)

## 1D PDE Dirichlet boundary conditions
println("1D PDE Dirichlet boundary conditions")
@parameters t x
@variables u(..)
Dt = Differential(t)
Dxx = Differential(x)^2

eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
bcs = [u(0,x) ~ cos(x),
        u(t,0) ~ exp(-t),
        u(t,1) ~ exp(-t) * cos(1)]

domains = [t ∈ Interval(0.0,1.0),
          x ∈ Interval(0.0,1.0)]

@named pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t, x)])

inner = 30
chain = FastChain(FastDense(2,inner,Flux.σ),
                  FastDense(inner,inner,Flux.σ),
                  FastDense(inner,inner,Flux.σ),
                  FastDense(inner,inner,Flux.σ),
                  FastDense(inner,inner,Flux.σ),
                  FastDense(inner,inner,Flux.σ),
                  FastDense(inner,1))#,(u,p)->gpuones .* u)

strategy = NeuralPDE.StochasticTraining(500)
initθ = CuArray(Float64.(DiffEqFlux.initial_params(chain)))
discretization = NeuralPDE.PhysicsInformedNN(chain,
                                             strategy;
                                             init_params = initθ)
prob = NeuralPDE.discretize(pdesys,discretization)
symprob = NeuralPDE.symbolic_discretize(pdesys,discretization)

res = GalacticOptim.solve(prob, ADAM(0.01);maxiters=1000)
prob = remake(prob,u0=res.minimizer)
res = GalacticOptim.solve(prob,ADAM(0.001);maxiters=1000)
phi = discretization.phi

u_exact = (t,x) -> exp.(-t) * cos.(x)
ts,xs = [infimum(d.domain):0.01:supremum(d.domain) for d in domains]
u_predict = reshape([first(Array(phi([t,x],res.minimizer))) for t in ts for x in xs],(length(ts),length(xs)))
u_real = reshape([u_exact(t,x) for t in ts  for x in xs ], (length(ts),length(xs)))
diff_u = abs.(u_predict .- u_real)

@test u_predict ≈ u_real atol = 1.0

# p1 = plot(ts, xs, u_real, linetype=:contourf,title = "analytic");
# p2 = plot(ts, xs, u_predict, linetype=:contourf,title = "predict");
# p3 = plot(ts, xs, diff_u,linetype=:contourf,title = "error");
# plot(p1,p2,p3)

## 1D PDE Neumann boundary conditions and Float64 accuracy
println("1D PDE Neumann boundary conditions and Float64 accuracy")
@parameters t x
@variables u(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

# 1D PDE and boundary conditions
eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
bcs = [u(0,x) ~ cos(x),
        Dx(u(t,0)) ~ 0.0,
        Dx(u(t,1)) ~ -exp(-t) * sin(1.0)]

# Space and time domains
domains = [t ∈ Interval(0.0,1.0),
        x ∈ Interval(0.0,1.0)]

# PDE system
@named pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t, x)])

inner = 20
chain = FastChain(FastDense(2,inner,Flux.σ),
                  FastDense(inner,inner,Flux.σ),
                  FastDense(inner,inner,Flux.σ),
                  FastDense(inner,inner,Flux.σ),
                  FastDense(inner,1))

strategy = NeuralPDE.QuasiRandomTraining(500; #points
                                         sampling_alg = SobolSample(),
                                         resampling = false,
                                         minibatch = 30)
initθ = CuArray(Float64.(DiffEqFlux.initial_params(chain)))
discretization = NeuralPDE.PhysicsInformedNN(chain,
                                             strategy;
                                             init_params = initθ)
prob = NeuralPDE.discretize(pdesys,discretization)
symprob = NeuralPDE.symbolic_discretize(pdesys,discretization)

res = GalacticOptim.solve(prob, ADAM(0.1); maxiters=2000)
prob = remake(prob,u0=res.minimizer)
res = GalacticOptim.solve(prob,ADAM(0.01); maxiters=2000)
phi = discretization.phi

u_exact = (t,x) -> exp(-t) * cos(x)
ts,xs = [infimum(d.domain):0.01:supremum(d.domain) for d in domains]
u_predict = reshape([first(Array(phi([t,x],res.minimizer))) for t in ts for x in xs],(length(ts),length(xs)))
u_real = reshape([u_exact(t,x) for t in ts  for x in xs ], (length(ts),length(xs)))
diff_u = abs.(u_predict .- u_real)

@test u_predict ≈ u_real atol = 1.0

# p1 = plot(ts, xs, u_real, linetype=:contourf,title = "analytic");
# p2 = plot(ts, xs, u_predict, linetype=:contourf,title = "predict");
# p3 = plot(ts, xs, diff_u,linetype=:contourf,title = "error");
# plot(p1,p2,p3)

## 2D PDE
println("2D PDE")
@parameters t x y
@variables u(..)
Dxx = Differential(x)^2
Dyy = Differential(y)^2
Dt = Differential(t)
t_min= 0.0
t_max = 2.0
x_min = 0.0
x_max = 2.0
y_min = 0.0
y_max = 2.0

# 3D PDE
eq  = Dt(u(t,x,y)) ~ Dxx(u(t,x,y)) + Dyy(u(t,x,y))

analytic_sol_func(t,x,y) = exp(x+y)*cos(x+y+4t)
# Initial and boundary conditions
bcs = [u(t_min,x,y) ~ analytic_sol_func(t_min,x,y),
       u(t,x_min,y) ~ analytic_sol_func(t,x_min,y),
       u(t,x_max,y) ~ analytic_sol_func(t,x_max,y),
       u(t,x,y_min) ~ analytic_sol_func(t,x,y_min),
       u(t,x,y_max) ~ analytic_sol_func(t,x,y_max)]

# Space and time domains
domains = [t ∈ Interval(t_min,t_max),
           x ∈ Interval(x_min,x_max),
           y ∈ Interval(y_min,y_max)]

# Neural network
inner = 25
chain = FastChain(FastDense(3,inner,Flux.σ),
                  FastDense(inner,inner,Flux.σ),
                  FastDense(inner,inner,Flux.σ),
                  FastDense(inner,inner,Flux.σ),
                  FastDense(inner,1))#,(u,p)->gpuones .* u)

initθ = CuArray(Float64.(DiffEqFlux.initial_params(chain)))

strategy = NeuralPDE.GridTraining(0.05)
discretization = NeuralPDE.PhysicsInformedNN(chain,
                                             strategy;
                                             init_params = initθ)

@named pde_system = PDESystem(eq,bcs,domains,[t,x,y],[u(t,x,y)])
prob = NeuralPDE.discretize(pde_system,discretization)
symprob = NeuralPDE.symbolic_discretize(pde_system,discretization)

res = GalacticOptim.solve(prob,ADAM(0.01);maxiters=2500)
prob = remake(prob,u0=res.minimizer)
res = GalacticOptim.solve(prob,ADAM(0.001);maxiters=2500)
@show res.original

phi = discretization.phi
ts,xs,ys = [infimum(d.domain):0.1:supremum(d.domain) for d in domains]
u_real = [analytic_sol_func(t,x,y) for t in ts for x in xs for y in ys]
u_predict = [first(Array(phi([t, x, y], res.minimizer))) for t in ts for x in xs for y in ys]

@test u_predict ≈ u_real rtol = 0.2

# using Plots
# using Printf
#
# function plot_(res)
#     # Animate
#     anim = @animate for (i, t) in enumerate(0:0.05:t_max)
#         @info "Animating frame $i..."
#         u_real = reshape([analytic_sol_func(t,x,y) for x in xs for y in ys], (length(xs),length(ys)))
#         u_predict = reshape([Array(phi([t, x, y], res.minimizer))[1] for x in xs for y in ys], length(xs), length(ys))
#         u_error = abs.(u_predict .- u_real)
#         title = @sprintf("predict t = %.3f", t)
#         p1 = plot(xs, ys, u_predict,st=:surface, label="", title=title)
#         title = @sprintf("real")
#         p2 = plot(xs, ys, u_real,st=:surface, label="", title=title)
#         title = @sprintf("error")
#         p3 = plot(xs, ys, u_error, st=:contourf,label="", title=title)
#         plot(p1,p2,p3)
#     end
#     gif(anim,"3pde.gif", fps=10)
# end
#
# plot_(res)


text:
sqlitetype{T<:Integer}(::Type{T}) = "INT"
sqlitetype{T<:AbstractFloat}(::Type{T}) = "REAL"
sqlitetype{T<:AbstractString}(::Type{T}) = "TEXT"
sqlitetype(::Type{NullType}) = "NULL"
sqlitetype(x) = "BLOB"

function createtable!(db::DB, name::AbstractString, schema::Data.Schema; temp::Bool=false, ifnotexists::Bool=true)
    rows, cols = size(schema)
    temp = temp ? "TEMP" : ""
    ifnotexists = ifnotexists ? "IF NOT EXISTS" : ""
    header, types = Data.header(schema), Data.types(schema)
    columns = [string(esc_id(header[i]), ' ', sqlitetype(types[i])) for i = 1:cols]
    SQLite.execute!(db, "CREATE $temp TABLE $ifnotexists $(esc_id(name)) ($(join(columns, ',')))")
    return name
end

"""
independent SQLite.Sink constructor to create a new or wrap an existing SQLite table with name `name`.
must provide a `Data.Schema` through the `schema` argument
can optionally provide an existing SQLite table name or new name that a created SQLite table will be called through the `name` argument
`temp=true` will create a temporary SQLite table that will be destroyed automatically when the database is closed
`ifnotexists=false` will throw an error if `name` already exists in `db`
"""
function Sink(db::DB, name::AbstractString, schema::Data.Schema=Data.Schema(); temp::Bool=false, ifnotexists::Bool=true, append::Bool=false)
    cols = size(SQLite.query(db, "pragma table_info($name)"), 1)
    if cols == 0
        createtable!(db, name, schema)
        cols = size(SQLite.query(db, "pragma table_info($name)"), 1)
    else
        !append && execute!(db, "delete from $(esc_id(name))")
    end
    params = chop(repeat("?,", cols))
    stmt = SQLite.Stmt(db, "INSERT INTO $(esc_id(name)) VALUES ($params)")
    return Sink(db, name, stmt, "")
end

# DataStreams interface
Data.streamtypes{T<:SQLite.Sink}(::Type{T}) = [Data.Field]

function Sink{T}(sch::Data.Schema, ::Type{T}, append::Bool, ref::Vector{UInt8}, db::DB, name::AbstractString; kwargs...)
    sink = Sink(db, name, sch; append=append, kwargs...)
    execute!(sink.db, "PRAGMA synchronous = OFF;")
    sink.transaction = string("SQLITE",randstring(10))
    transaction(sink.db, sink.transaction)
    return sink
end
function Sink{T}(sink, sch::Data.Schema, ::Type{T}, append::Bool, ref::Vector{UInt8})
    execute!(sink.db, "PRAGMA synchronous = OFF;")
    sink.transaction = string("SQLITE", randstring(10))
    transaction(sink.db, sink.transaction)
    !append && execute!(sink.db, "delete from $(esc_id(sink.tablename))")
    return sink
end

function getbind!{T}(val::Nullable{T}, col, stmt)
    if isnull(val)
        SQLite.bind!(stmt, col, NULL)
    else
        SQLite.bind!(stmt, col, get(val))
    end
    return
end
if isdefined(:isna)
    function getbind!{T}(val::T, col, stmt)
        if isna(val)
            SQLite.bind!(stmt, col, NULL)
        else
            SQLite.bind!(stmt, col, val)
        end
        return
    end
else
    getbind!{T}(val::T, col, stmt) = SQLite.bind!(stmt, col, val)
end

function Data.streamto!{T}(sink::SQLite.Sink, ::Type{Data.Field}, val::T, row, col, sch)
    getbind!(val, col, sink.stmt)
    if col == size(sch, 2)
        SQLite.sqlite3_step(sink.stmt.handle)
        SQLite.sqlite3_reset(sink.stmt.handle)
    end
    return nothing
end

function Data.cleanup!(sink::SQLite.Sink)
    rollback(sink.db, sink.transaction)
    commit(sink.db, sink.transaction)
    execute!(sink.db, "PRAGMA synchronous = ON;")
    return nothing
end

function Data.close!(sink::SQLite.Sink)
    commit(sink.db, sink.transaction)
    execute!(sink.db, "PRAGMA synchronous = ON;")
    SQLite.execute!(sink.db, "ANALYZE $(esc_id(sink.tablename))")
    return nothing
end

"""
Load a Data.Source `source` into an SQLite table that will be named `tablename` (will be auto-generated if not specified).

`temp=true` will create a temporary SQLite table that will be destroyed automatically when the database is closed
`ifnotexists=false` will throw an error if `tablename` already exists in `db`
"""
function load{T}(db::SQLite.DB, name, ::Type{T}, args...; append::Bool=false, transforms::Dict=Dict{Int,Function}(), kwargs...)
    sink = Data.stream!(T(args...), SQLite.Sink, append, transforms, db, name; kwargs...)
    Data.close!(sink)
    return sink
end
function load{T}(db::SQLite.DB, name, source::T; append::Bool=false, transforms::Dict=Dict{Int,Function}(), kwargs...)
    sink = Data.stream!(source, SQLite.Sink, append, transforms, db, name; kwargs...)
    Data.close!(sink)
    return sink
end

load{T}(sink::Sink, ::Type{T}, args...; append::Bool=false, transforms::Dict=Dict{Int,Function}()) = (sink = Data.stream!(T(args...), sink, append, transforms); Data.close!(sink); return sink)
load(sink::Sink, source; append::Bool=false, transforms::Dict=Dict{Int,Function}()) = (sink = Data.stream!(source, sink, append, transforms); Data.close!(sink); return sink)

max_stars_repo_name: beffaCo/updateSqlite

max_issues_repo_name: beffaCo/updateSqlite

max_forks_repo_name: beffaCo/updateSqlite


text:
export scal,param,a,adag,f,fdag,OpSumAnalytic,ExpVal,Corr
export @boson_ops,@fermion_ops
export σx,σy,σz,σp,σm
export @Pr_str, @Pc_str, ∑

# define for σx, σy, σz
@enum SpatialIndex x=1 y=2 z=3
SpatialIndex(a::SpatialIndex) = a

const OpIndex = Union{Int,Symbol}
const OpIndices = NTuple{N,OpIndex} where N

abstract type Operator end
abstract type Scalar <: Operator; end
"Represent a scalar value (i.e., a number)"
struct scal{T<:Number} <: Scalar; v::T; end
"`param(g,(:i,:j),'n')`: represent a scalar named parameter ``g_{i,j}``. state can be purely real (`'r'`), not conjugated (`'n'`), or conjugated (`'c'`)"
struct param{T<:OpIndices} <: Scalar
    name::Symbol
    state::Char
    inds::T
    function param(name,state::Char,inds::OpIndex...)
        state in ('n','r','c') || throw(ArgumentError("state has to be one of n,r,c"))
        new{typeof(inds)}(name,state,inds)
    end
    param(name,inds::OpIndex...) = param(name,'n',inds...)
    param(name,state::Char,inds::Tuple) = param(name,state,inds...)
    param(name,inds::Tuple) = param(name,'n',inds...)
end
struct δ <: Scalar
    inds::Tuple{OpIndex,OpIndex}
    δ(inds) = δ(inds...)
    function δ(iA::OpIndex,iB::OpIndex)
        # sort indices
        if iA == iB
            scal(1)
        elseif iA isa Integer && iB isa Integer
            # e.g., δ_1,3 = 0 (integers are not symbolic!)
            scal(0)
        elseif sortsentinel(iB) < sortsentinel(iA)
            new((iB,iA))
        else
            new((iA,iB))
        end
    end
end
function δ(Ainds::OpIndices,Binds::OpIndices)::Operator
    if length(Ainds) == length(Binds)
        prod(collect(Operator,δ.(Ainds,Binds)))
    else
        scal(0)
    end
end

abstract type BaseOperator <: Operator; end
for (op,desc) in (
    (:BosonDestroy,"bosonic annihilation"),
    (:BosonCreate,"bosonic creation"),
    (:FermionDestroy,"fermionic annihilation"),
    (:FermionCreate,"fermionic creation"))
    @eval begin
        "`$($op)(name,inds)`: represent $($desc) operator ``name_{inds}``"
        struct $op{T<:OpIndices} <: BaseOperator
            name::Symbol
            inds::T
            $op(name::Symbol,inds::OpIndex...) = new{typeof(inds)}(name,inds)
            $op(name::Symbol,inds::Tuple) = $op(name,inds...)
        end
    end
end

"`@boson_ops name`: return functions for creating bosonic annihilation and creation operators with name `name` (i.e., wrappers of [`BosonDestroy`](@ref) and [`BosonCreate`](@ref))"
macro boson_ops(name)
    quote
        ann(args...) = BosonDestroy($(Meta.quot(name)),args...)
        cre(args...) = BosonCreate($(Meta.quot(name)),args...)
        ann,cre
    end
end

"`@fermion_ops name`: return functions for creating fermionic annihilation and creation operators with name `name` (i.e., wrappers of [`FermionDestroy`](@ref) and [`FermionCreate`](@ref))"
macro fermion_ops(name)
    quote
        ann(args...) = FermionDestroy($(Meta.quot(name)),args...)
        cre(args...) = FermionCreate($(Meta.quot(name)),args...)
        ann,cre
    end
end

a, adag = @boson_ops a
f, fdag = @fermion_ops f

for (op,desc,sym) in (
    (:σminus,"TLS annihilation","σ^-"),
    (:σplus,"TLS creation","σ^+"))
    @eval begin
        "`$($op)(inds)`: represent $($desc) operator ``$($sym)_{inds}``"
        struct $op{T<:OpIndices} <: BaseOperator
            inds::T
            $op(inds::OpIndex...) = new{typeof(inds)}(inds)
            $op(inds::Tuple) = $op(inds...)
        end
    end
end

"`σ(a,inds)`: represent Pauli matrix ``σ_{a,inds}`` for two-level system (TLS), where ``a ∈ \\{x,y,z\\}``."
struct σ{T<:OpIndices} <: BaseOperator
    a::SpatialIndex
    inds::T
    σ(a,inds::OpIndex...) = new{typeof(inds)}(SpatialIndex(a),inds)
    σ(a,inds::Tuple) = σ(a,inds...)
end

struct OpProd <: Operator; A::Operator; B::Operator; end
struct OpSum  <: Operator; A::Operator; B::Operator; end

"`OpSumAnalytic(i::Symbol,A::Operator)` or `∑(i,A)`: represent ``\\sum_{i} A``, with all possible values of ``i`` assumed to be included"
struct OpSumAnalytic <: Operator
    ind::Symbol
    A::Operator
    OpSumAnalytic(ind::Symbol,A::scal) = A.v == 0 ? A : new(ind,A)
    OpSumAnalytic(ind::Symbol,A::Operator) = begin
        if A isa OpProd && A.A isa scal
            A.A*OpSumAnalytic(ind,A.B)
        elseif A isa OpSumAnalytic && A.ind < ind
            # reorder nested sums by indices
            OpSumAnalytic(A.ind,OpSumAnalytic(ind,A.A))
        else
            for t in preftuple(A)
                # if there is a δ in the expression and it has the sum index, the sum disappears
                if t isa δ && ind in t.inds
                    irepl = ind==t.inds[1] ? t.inds[2] : t.inds[1]
                    return replace_index(A,ind,irepl)
                end
            end
            return new(ind,A)
        end
    end
    OpSumAnalytic(ind::Symbol,A::OpSum) = OpSumAnalytic(ind,A.A) + OpSumAnalytic(ind,A.B)
end

const ∑ = OpSumAnalytic

"`ExpVal(A::Operator)`: represent expectation value ``⟨A⟩``"
struct ExpVal <: Scalar
    A::Operator
    ExpVal(A::Scalar) = A
    ExpVal(A::OpSum) = ExpVal(A.A) + ExpVal(A.B)
    ExpVal(A::OpProd) = A.A isa Scalar ? A.A*ExpVal(A.B) : new(A)
    ExpVal(A::OpSumAnalytic) = OpSumAnalytic(A.ind,ExpVal(A.A))
    ExpVal(A::Operator) = new(A)
end

"""
    Corr(A::Operator)

Represents correlations of operator ``A``. ``A`` should be a product for this to
make sense, in which case ``⟨AB⟩_c = ⟨AB⟩ - ⟨A⟩⟨B⟩``, with corresponding
extensions for products of more operators.

See also: [`ascorr`](@ref)
"""
struct Corr <: Scalar
    A::Operator
    Corr(A::Scalar) = A
    Corr(A::OpSum) = Corr(A.A) + Corr(A.B)
    Corr(A::OpProd) = A.A isa Scalar ? A.A*Corr(A.B) : new(A)
    Corr(A::Operator) = new(A)
end

## functions for constructing Pauli operators depending on whether we use (σ+,σ-) or (σx,σy,σz) as the "basic" operators

const _using_σpm = Ref(false)
use_σpm(t::Bool=true) = (_using_σpm[] = t; nothing)
use_σxyz() = use_σpm(false)
using_σpm() = _using_σpm[]

"`σp(n)`: construct ``σ^+_n = \\frac12 σ_{x,n} + \\frac{i}{2} σ_{y,n}``"
σp(n...) = using_σpm() ? σplus(n...) : scal(1//2)*σ(x,n...) + scal(1im//2)*σ(y,n...)
"`σm(n)`: construct ``σ^-_n = \\frac12 σ_{x,n} - \\frac{i}{2} σ_{y,n}``"
σm(n...) = using_σpm() ? σminus(n...) : scal(1//2)*σ(x,n...) - scal(1im//2)*σ(y,n...)
"`σx(n)`: construct ``σ_{x,n}``"
σx(n...) = using_σpm() ? σminus(n...) + σplus(n...) : σ(x,n...)
"`σy(n)`: construct ``σ_{y,n}``"
σy(n...) = using_σpm() ? scal(1im) * (σminus(n...) - σplus(n...)) : σ(y,n...)
"`σz(n)`: construct ``σ_{z,n}``"
σz(n...) = using_σpm() ? scal(2)*σplus(n...)*σminus(n...) - scal(1) : σ(z,n...)

## functions for constructing `param`s with string macros,
# Pc"ω_i,j" = param(:ω,'n',:i,:j) (complex parameter)
# Pr"ω_i,j" = param(:ω,'r',:i,:j) (real parameter)

function parse_paramstr(s)
    s = split(s,"_")
    par = Symbol(s[1])
    @assert length(s) <= 2
    inds = length(s)==1 ? () : Meta.parse.(split(s[2],","))
    par,inds
end

macro Pc_str(s)
    par, inds = parse_paramstr(s)
    param(par,'n',inds...)
end
macro Pr_str(s)
    par, inds = parse_paramstr(s)
    param(par,'r',inds...)
end


text:
using HTTP, JSON

@testset "Netcode tests" begin
    begin
        d = Dict((("type", "gameStart"),))
        s = JSON.json(d)
        io = IOBuffer(s)
        @test LichessBot.eventsCallback(io) == 1
    end # begin
    begin
        io = IOBuffer("{\"type\":\"gameState\", \"moves\":\"a2a3\", \"btime\":5000, \"binc\":1000}")
        @test isnothing(gameCallback(io))
    end
end


text:
export iauSepp
"""
Angular separation between two p-vectors.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Given:
   a      double[3]    first p-vector (not necessarily unit length)
   b      double[3]    second p-vector (not necessarily unit length)

Returned (function value):
         double       angular separation (radians, always positive)

Notes:

   1. If either vector is null, a zero result is returned.

   2. The angular separation is most simply formulated in terms of
   scalar product.  However, this gives poor accuracy for angles
   near zero and pi.  The present algorithm uses both cross product
   and dot product, to deliver full accuracy whatever the size of
   the angle.

Called:
   iauPxp       vector product of two p-vectors
   iauPm        modulus of p-vector
   iauPdp       scalar product of two p-vectors

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# double iauSepp(double a[3], double b[3])

function iauSepp(a::AbstractVector{<:Real}, b::AbstractVector{<:Real})
   return ccall((:iauSepp, libsofa_c), Cdouble,
       (Ptr{Cdouble}, Ptr{Cdouble}),
       convert(Array{Float64, 1}, a),
       convert(Array{Float64, 1}, b))
end


text:
# This file was generated, do not modify it. # hide
avg = MyAverageTwo(regressor1=tree_booster,
                   regressor2=svm)

stack = MyTwoModelStack(regressor1=tree_booster,
                        regressor2=svm,
                        judge=forest) # judge=linear

all_models = [tree_booster, svm, forest, avg, stack];

for model in all_models
    print_performance(model, X, y)
end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

function first_8_words(
            compression_output::SVector{16, UInt32}
        )::SVector{8, UInt32}
    return SVector(
        compression_output[1], compression_output[2],
        compression_output[3], compression_output[4],
        compression_output[5], compression_output[6],
        compression_output[7], compression_output[8]
    )
end

struct Output
    input_chaining_value::SVector{8, UInt32}
    block_words::SVector{16, UInt32}
    counter::UInt64
    block_len::UInt32
    flags::UInt32
end

function chaining_value(self::Output)::SVector{8, UInt32}
    return first_8_words(compress(
        self.input_chaining_value,
        self.block_words,
        self.counter,
        self.block_len,
        self.flags
    ))
end

function root_output_bytes(
            self::Output,
            out_slice::AbstractVector{UInt8}
        )::Nothing
    local output_block_counter = UInt64(0)

    for out_block in partition(out_slice, 2 * OUT_LEN)
        local words = compress(
            self.input_chaining_value,
            self.block_words,
            output_block_counter,
            self.block_len,
            self.flags | ROOT
        )
        for (word, out_word) in zip(words, partition(out_block, 4))
            out_word[1] = UInt8(word & 0xff)
            if length(out_word) >= 2
                out_word[2] = UInt8((word >> 8) & 0xff)
            end
            if length(out_word) >= 3
                 out_word[3] = UInt8((word >> 16) & 0xff)
            end
            if length(out_word) >= 4
                out_word[4] = UInt8((word >> 24) & 0xff)
            end
        end
        output_block_counter += 1
    end
    nothing
end