text:
export meataxe, composition_factors, composition_series, submodules, maximal_submodules, minimal_submodules

####################################################################
#
#  Tools for MeatAxe
#
#####################################################################

#
# Given a matrix $M$ in echelon form and a vector, it returns
# the vector reduced with respect to $M$
#
function cleanvect(M::fq_nmod_mat, v::fq_nmod_mat)

  @assert rows(v)==1
  w=deepcopy(v)
  if iszero(v)
    return w
  end
  for i=1:rows(M)
    ind=1
    while M[i,ind]==0
      ind+=1
    end
    if iszero(w[1,ind])
      continue
    end
    mult=w[1,ind]//M[i,ind]
    w[1,ind]=parent(M[1,1])(0)
    for k=ind+1:cols(M)
      w[1,k]-= mult*M[i,k]
    end
  end
  return w

end

#
#  Given a matrix C containing the coordinates of vectors v_1,dots, v_k
#  in echelon form, the function computes a basis for the submodule they generate
#

function closure(C::fq_nmod_mat, G::Array{fq_nmod_mat,1})

  rref!(C)
  i=1
  while i <= rows(C)
    w=sub(C, i:i, 1:cols(C))
    for j=1:length(G)
      res=cleanvect(C,w*G[j])
      if !iszero(res)
        C=vcat(C,res)
        if rows(C)==cols(C)
          i=cols(C)+1
          break
        end
      end
    end
    i+=1
  end
  rref!(C)
  return C

end

#
#  Given a matrix C containing the coordinates of vectors v_1,dots, v_k,
#  the function computes a basis for the submodule they generate
#

function spinning(C::fq_nmod_mat,G::Array{fq_nmod_mat,1})

  B=deepcopy(C)
  X=rref(C)[2]
  i=1
  while i != rows(B)+1
    for j=1:length(G)
      el= sub(B, i:i, 1:cols(B)) * G[j]
      res= cleanvect(X,el)
      if !iszero(res)
        X=vcat(X,res)
        rref!(X)
        B=vcat(B,el)
        if rows(B)==cols(B)
          return B
        end
      end
    end
    i+=1
  end
  return B

end

#
#  Function to obtain the action of G on the quotient and on the submodule
#

function clean_and_quotient(M::fq_nmod_mat,N::fq_nmod_mat, pivotindex::Set{Int})

  coeff=zero_matrix(parent(M[1,1]),rows(N),rows(M))
  for i=1:rows(N)
    for j=1:rows(M)
      ind=1
      while iszero(M[j,ind])
        ind+=1
      end
      coeff[i,j]=N[i,ind]//M[j,ind]
      for s=1:cols(N)
        N[i,s]-=coeff[i,j]*M[j,s]
      end
    end
  end
  vec= zero_matrix(parent(M[1,1]),rows(N),cols(M)-length(pivotindex))
  for i=1:rows(N)
    pos=0
    for s=1:cols(M)
      if !(s in pivotindex)
        pos+=1
        vec[i,pos]=N[i,s]
      end
    end
  end
  return coeff, vec
end

#
#  Restriction of the action to the submodule generated by C and the quotient
#

function _split(C::fq_nmod_mat,G::Array{fq_nmod_mat,1})
# I am assuming that C is a Fp[G]-submodule

  equot=Array{fq_nmod_mat,1}(length(G))
  esub=Array{fq_nmod_mat,1}(length(G))
  pivotindex=Set{Int}()
  for i=1:rows(C)
    ind=1
    while iszero(C[i,ind])
      ind+=1
    end
    push!(pivotindex,ind)
  end
  for a=1:length(G)
    subm,vec=clean_and_quotient(C, C*G[a],pivotindex)
    esub[a]=subm
    s=zero_matrix(parent(C[1,1]),cols(G[1])-length(pivotindex),cols(G[1])-length(pivotindex))
    pos=0
    for i=1:rows(G[1])
      if !(i in pivotindex)
        m,vec=clean_and_quotient(C,sub(G[a],i:i,1:rows(G[1])),pivotindex)
        for j=1:cols(vec)
          s[i-pos,j]=vec[1,j]
        end
      else
        pos+=1
      end
    end
    equot[a]=s
  end
  return FqGModule(esub),FqGModule(equot),pivotindex

end

#
#  Restriction of the action to the submodule generated by C
#

function actsub(C::fq_nmod_mat,G::Array{fq_nmod_mat,1})

  esub=Array{fq_nmod_mat,1}(length(G))
  pivotindex=Set{Int}()
  for i=1:rows(C)
    ind=1
    while iszero(C[i,ind])
      ind+=1
    end
    push!(pivotindex,ind)
  end
  for a=1:length(G)
    subm,vec=clean_and_quotient(C, C*G[a],pivotindex)
    esub[a]=subm
  end
  return FqGModule(esub)
end

#
#  Restriction of the action to the quotient by the submodule generated by C
#

function actquo(C::fq_nmod_mat,G::Array{fq_nmod_mat,1})

  equot=Array{fq_nmod_mat,1}(length(G))
  pivotindex=Set{Int}()
  for i=1:rows(C)
    ind=1
    while iszero(C[i,ind])
      ind+=1
    end
    push!(pivotindex,ind)
  end
  for a=1:length(G)
    s=zero_matrix(parent(C[1,1]),cols(G[1])-length(pivotindex),cols(G[1])-length(pivotindex))
    pos=0
    for i=1:rows(G[1])
      if !(i in pivotindex)
        m,vec=clean_and_quotient(C,sub(G[a],i:i,1:rows(G[1])),pivotindex)
        for j=1:cols(vec)
          s[i-pos,j]=vec[1,j]
        end
      else
        pos+=1
      end
    end
    equot[a]=s
  end
  return FqGModule(equot), pivotindex

end

#
#  Function that determine if two G-modules are isomorphic, provided that the first is irreducible
#

function isisomorphic(M::FqGModule,N::FqGModule)

  @assert M.isirreducible
  @assert M.K==N.K
  @assert length(M.G)==length(N.G)
  if M.dim!=N.dim
    return false
  end

  if M.dim==1
    return M.G==N.G
  end

  K=M.K
  Kx,x=PolynomialRing(K, "x", cached=false)

  if length(M.G)==1
    f=charpoly(Kx,M.G[1])
    g=charpoly(Kx,N.G[1])
    if f==g
      return true
    else
      return false
    end
  end

  n=M.dim
  posfac=n

  f=Kx(1)
  G=deepcopy(M.G)
  H=deepcopy(N.G)

  rel=_relations(M,N)
  return iszero(rel[N.dim, N.dim])

  #=

  #
  #  Adding generators to obtain randomness
  #

  for i=1:max(length(M.G),9)
    l1=rand(1:length(G))
    l2=rand(1:length(G))
    while l1 !=l2
      l2=rand(1:length(G))
    end
    push!(G, G[l1]*G[l2])
    push!(H, H[l1]*H[l2])
  end

  #
  #  Now, we get peakwords
  #

  A=zero_matrix(K,n,n)
  B=zero_matrix(K,n,n)
  found=false

  while !found

    A=zero_matrix(K,n,n)
    B=zero_matrix(K,n,n)
    l1=rand(1:length(G))
    l2=rand(1:length(G))
    push!(G, G[l1]*G[l2])
    push!(H, H[l1]*H[l2])

    for i=1:length(G)
      s=rand(K)
      A+=s*G[i]
      B+=s*H[i]
    end

    cp=charpoly(Kx,A)
    cpB=charpoly(Kx,B)
    if cp!=cpB
      return false
    end
    sq=prod(collect(keys(factor_squarefree(cp).fac)))
    j=1
    while !isone(sq)
      g=gcd(x^(Int(order(K)^j))-x,sq)
      sq=divexact(sq,g)
      lf=factor(g)
      for t in keys(lf.fac)
        f=t
        S=_subst(t,A)
        a,kerA=nullspace(transpose(S))
        if a==1
          M.dim_spl_fld=1
          found=true
          break
        end
        kerA=transpose(kerA)
        posfac=gcd(posfac,a)
        if divisible(fmpz(posfac),a)
          v=sub(kerA, 1:1, 1:n)
          U=v
          T =spinning(v,G)
          G1=[T*mat*inv(T) for mat in M.G]
          i=2
          E=fq_nmod_mat[eye(T,a)]
          while rows(U)!= a
            w= sub(kerA, i:i, 1:n)
            z= cleanvect(U,w)
            if iszero(z)
              continue
            end
            O =spinning(w,G)
            G2=[O*mat*inv(O) for mat in M.G]
            if G1 == G2
              b=kerA*O
              x=transpose(solve(transpose(kerA),transpose(b)))
              push!(E,x)
              U=vcat(U,z)
              U=closure(U,E)
            else
              break
            end
            if rows(U)==a
              M.dim_spl_fld=a
              found=true
              break
            else
              i+=1
            end
          end
        end
        if found==true
          break
        end
      end
      j+=1
    end
  end
  #
  #  Get the standard basis
  #

  L=_subst(f,A)
  a,kerA=nullspace(transpose(L))

  I=_subst(f,B)
  b,kerB=nullspace(transpose(I))

  if a!=b
    return false
  end

  Q= spinning(transpose(sub(kerA, 1:n, 1:1)), M.G)
  W= spinning(transpose(sub(kerB, 1:n, 1:1)), N.G)

  #
  #  Check if the actions are conjugated
  #
  S=inv(W)*Q
  T=inv(S)
  for i=1:length(M.G)
    if S*M.G[i]* T != N.G[i]
      return false
    end
  end
  return true

  =#
end

function _solve_unique(A::fq_nmod_mat, B::fq_nmod_mat)
  X = zero_matrix(base_ring(A), cols(B), rows(A))

  #println("solving\n $A \n = $B * X")
  r, per, L, U = lufact(B) # P*M1 = L*U
  inv!(per)
  @assert B == per*L*U

  Ap = inv(per)*A
  Y = similar(A)

  #println("first solve\n $Ap = $L * Y")

  for i in 1:cols(Y)
    for j in 1:rows(Y)
      s = Ap[j, i]
      for k in 1:j-1
        s = s - Y[k, i]*L[j, k]
      end
      Y[j, i] = s
    end
  end

  @assert Ap == L*Y

  #println("solving \n $Y \n = $U * X")

  YY = sub(Y, 1:r, 1:cols(Y))
  UU = sub(U, 1:r, 1:r)
  X = inv(UU)*YY

  @assert Y == U * X

  @assert B*X == A
  return X
end

function _enum_el(K,v,dim)

  if dim == 0
    return [v]
  else
    list=[]
    push!(v,K(0))
    for x in K
      v[length(v)]=x
      push!(list,deepcopy(v))
    end
    list1=[]
    for x in list
      append!(list1,_enum_el(K,x, dim-1))
    end
    return list1
  end
end

function dual_space(M::FqGModule)

  G=fq_nmod_mat[transpose(g) for g in M.G]
  return FqGModule(G)

end

function _subst(f::Nemo.PolyElem{T}, a::fq_nmod_mat) where {T <: Nemo.RingElement}
   #S = parent(a)
   n = degree(f)
   if n < 0
      return similar(a)#S()
   elseif n == 0
      return coeff(f, 0)*eye(a)
   elseif n == 1
      return coeff(f, 0)*eye(a) + coeff(f, 1)*a
   end
   d1 = isqrt(n)
   d = div(n, d1)
   A = powers(a, d)
   s = coeff(f, d1*d)*A[1]
   for j = 1:min(n - d1*d, d - 1)
      c = coeff(f, d1*d + j)
      if !iszero(c)
         s += c*A[j + 1]
      end
   end
   for i = 1:d1
      s *= A[d + 1]
      s += coeff(f, (d1 - i)*d)*A[1]
      for j = 1:min(n - (d1 - i)*d, d - 1)
         c = coeff(f, (d1 - i)*d + j)
         if !iszero(c)
            s += c*A[j + 1]
         end
      end
   end
   return s
end

#################################################################
#
#  MeatAxe, Composition Factors and Composition Series
#
#################################################################

doc"""
***
    meataxe(M::FqGModule) -> Bool, MatElem

> Given module M, returns true if the module is irreducible (and the identity matrix) and false if the space is reducible, togheter with a basis of a submodule

"""

function meataxe(M::FqGModule)

  K=M.K
  Kx,x=PolynomialRing( K,"x", cached=false)
  n=M.dim
  H=M.G
  if M.dim==1
    M.isirreducible=true
    return true, eye(H[1],n)
  end

  if length(H)==1
    A=H[1]
    poly=charpoly(Kx,A)
    sq=factor_squarefree(poly)
    lf=factor(first(keys(sq.fac)))
    t=first(keys(lf.fac))
    if degree(t)==n
      M.isirreducible=true
      return true, eye(H[1],n)
    else
      N= _subst(t, A)
      kern=transpose(nullspace(transpose(N))[2])
      B=closure(sub(kern,1:1, 1:n),H)
      return false, B
    end
  end

  #
  #  Adding generators to obtain randomness
  #
  G=deepcopy(H)
  Gt=fq_nmod_mat[transpose(x) for x in M.G]

  for i=1:max(length(M.G),9)
    l1=rand(1:length(G))
    l2=rand(1:length(G))
    while l1 !=l2
      l2=rand(1:length(G))
    end
    push!(G, G[l1]*G[l2])
  end

  while true

  # At every step, we add a generator to the group.

    push!(G, G[rand(1:length(G))]*G[rand(1:length(G))])

  #
  # Choose a random combination of the actual generators of G
  #
    A=zero_matrix(K,n,n)
    for i=1:length(G)
      A+=rand(K)*G[i]
    end

  #
  # Compute the characteristic polynomial and, for irreducible factor f, try the Norton test
  #
    poly=charpoly(Kx,A)
    sqfpart=keys(factor_squarefree(poly).fac)
    for el in sqfpart
      sq=el
      i=1
      while !isone(sq)
        f=gcd(powmod(x, order(K)^i, sq)-x,sq)
        sq=divexact(sq,f)
        lf=factor(f)
        for t in keys(lf.fac)
          N = _subst(t, A)
          a,kern=nullspace(transpose(N))
          #
          #  Norton test
          #
          B=closure(transpose(sub(kern,1:n, 1:1)),M.G)
          if rows(B)!=n
            M.isirreducible=false
            return false, B
          end
          kernt=nullspace(N)[2]
          Bt=closure(transpose(sub(kernt,1:n,1:1)),Gt)
          if rows(Bt)!=n
            subst=transpose(nullspace(Bt)[2])
            @assert rows(subst)==rows(closure(subst,G))
            M.isirreducible=false
            return false, subst
          end
          if degree(t)==a
            #
            # f is a good factor, irreducibility!
            #
            M.isirreducible=true
            return true, eye(G[1],n)
          end
        end
        i+=1
      end
    end
  end
end

doc"""
***
    composition_series(M::FqGModule) -> Array{MatElem,1}

> Given a Fq[G]-module M, it returns a composition series for M, i.e. a sequence of submodules such that the quotient of two consecutive element is irreducible.

"""

function composition_series(M::FqGModule)

  if isdefined(M, :isirreducible) && M.isirreducible==true
    return [eye(M.G[1],M.dim)]
  end

  bool, C = meataxe(M)
  #
  #  If the module is irreducible, we return a basis of the space
  #
  if bool ==true
    return [eye(M.G[1],M.dim)]
  end
  #
  #  The module is reducible, so we call the algorithm on the quotient and on the subgroup
  #
  G=M.G
  K=M.K

  rref!(C)

  esub,equot,pivotindex=_split(C,G)
  sub_list = composition_series(esub)
  quot_list = composition_series(equot)
  #
  #  Now, we have to write the submodules of the quotient and of the submodule in terms of our basis
  #
  list=Array{fq_nmod_mat,1}(length(sub_list)+length(quot_list))
  for i=1:length(sub_list)
    list[i]=sub_list[i]*C
  end
  for z=1:length(quot_list)
    s=zero_matrix(K,rows(quot_list[z]), cols(C))
    for i=1:rows(quot_list[z])
      pos=0
      for j=1:cols(C)
        if j in pivotindex
          pos+=1
        else
          s[i,j]=quot_list[z][i,j-pos]
        end
      end
    end
    list[length(sub_list)+z]=vcat(C,s)
  end
  return list
end

doc"""
***
    composition_factors(M::FqGModule)

> Given a Fq[G]-module M, it returns, up to isomorphism, the composition factors of M with their multiplicity,
> i.e. the isomorphism classes of modules appearing in a composition series of M

"""

function composition_factors(M::FqGModule; dimension::Int=-1)

  if isdefined(M, :isirreducible) && M.isirreducible
    if dimension!= -1
      if M.dim==dimension
        return Tuple{FqGModule, Int}[(M,1)]
      else
        return Tuple{FqGModule, Int}[]
      end
    else
      return Tuple{FqGModule, Int}[(M,1)]
    end
  end

  K=M.K::FqNmodFiniteField

  bool, C = meataxe(M)
  #
  #  If the module is irreducible, we just return a basis of the space
  #
  if bool
    if dimension!= -1
      if M.dim==dimension
        return Tuple{FqGModule, Int}[(M,1)]
      else
        return Tuple{FqGModule, Int}[]
      end
    else
      return Tuple{FqGModule, Int}[(M,1)]
    end
  end
  G=M.G
  #
  #  The module is reducible, so we call the algorithm on the quotient and on the subgroup
  #

  rref!(C)

  sub,quot,pivotindex=_split(C,G)
  sub_list = composition_factors(sub)
  quot_list = composition_factors(quot)
  #
  #  Now, we check if the factors are isomorphic
  #

  for i=1:length(sub_list)
    for j=1:length(quot_list)
      if isisomorphic(sub_list[i][1], quot_list[j][1])
        sub_list[i]=(sub_list[i][1], sub_list[i][2]+quot_list[j][2])
        deleteat!(quot_list,j)
        break
      end
    end
  end
  return append!(sub_list, quot_list)
  #=
  for i=1:length(sub_list)
    for j=1:length(quot_list)
      if isisomorphic(sub_list[i][1], quot_list[j][1])
        sub_list[i][2]+=quot_list[j][2]
        deleteat!(quot_list,j)
        break
      end
    end
  end
  return append!(sub_list,quot_list)
  =#
end

function _relations(M::FqGModule, N::FqGModule)

  @assert M.isirreducible
  G=M.G
  H=N.G
  K=M.K
  n=M.dim

  sys=zero_matrix(K,2*N.dim,N.dim)
  matrices=fq_nmod_mat[]
  first=true
  B=zero_matrix(K,1,M.dim)
  B[1,1]=K(1)
  X=B
  push!(matrices, eye(B,N.dim))
  i=1
  while i<=rows(B)
    w=sub(B, i:i, 1:n)
    for j=1:length(G)
      v=w*G[j]
      res=cleanvect(X,v)
      if !iszero(res)
        X=rref(vcat(X,v))[2]
        B=vcat(B,v)
        push!(matrices, matrices[i]*H[j])
      else
        x=_solve_unique(transpose(v),transpose(B))
        A=sum([x[q,1]*matrices[q] for q=1:rows(x)])
        A=A-(matrices[i]*H[j])
        if first
          for s=1:N.dim
            for t=1:N.dim
              sys[s,t]=A[t,s]
            end
          end
          first=false
        else
          for s=1:N.dim
            for t=1:N.dim
              sys[N.dim+s,t]=A[t,s]
            end
          end
        end
        rref!(sys)
      end
    end
    if sys[N.dim,N.dim]!=0
      break
    end
    i=i+1
  end
  return sys
end

function _irrsubs(M::FqGModule, N::FqGModule)

  @assert M.isirreducible

  K=M.K
  rel=_relations(M,N)
  if rel[N.dim, N.dim]!=0
    return fq_nmod_mat[]
  end
  a,kern=nullspace(rel)
  kern=transpose(kern)
  if a==1
    return fq_nmod_mat[closure(kern, N.G)]
  end
  #
  #  Reduce the number of homomorphism to try by considering the action of G on the homomorphisms
  #
  vects=fq_nmod_mat[sub(kern, i:i, 1:N.dim) for i=1:a]
  i=1
  while i<length(vects)
    X=closure(vects[i],N.G)
    j=i+1
    while j<= length(vects)
      if iszero(cleanvect(X,vects[j]))
        deleteat!(vects,j)
      else
        j+=1
      end
    end
    i+=1
  end
  if length(vects)==1
    return fq_nmod_mat[closure(vects[1], N.G)]
  end

  #
  #  Try all the possibilities. (A recursive approach? I don't know if it is a smart idea...)
  #  Notice that we eliminate lots of candidates considering the action of the group on the homomorphisms space
  #
  candidate_comb=append!(_enum_el(K,[K(0)], length(vects)-1),_enum_el(K,[K(1)],length(vects)-1))
  deleteat!(candidate_comb,1)
  list=Array{fq_nmod_mat,1}(length(candidate_comb))
  for j=1:length(candidate_comb)
    list[j] = sum([candidate_comb[j][i]*vects[i] for i=1:length(vects)])
  end
  final_list=fq_nmod_mat[]
  push!(final_list, closure(list[1], N.G))
  for i = 2:length(list)
    reduce=true
    for j=1:length(final_list)
      w=cleanvect(final_list[j],list[i])
      if iszero(w)
        reduce=false
        break
      end
    end
    if reduce
      push!(final_list, closure(list[i],N.G))
    end
  end
  return final_list

end

doc"""
***
    minimal_submodules(M::FqGModule)

> Given a Fq[G]-module M, it returns all the minimal submodules of M

"""

function minimal_submodules(M::FqGModule, dim::Int=M.dim+1, lf=[])

  K=M.K
  n=M.dim

  if M.isirreducible==true
    return fq_nmod_mat[]
  end

  list=fq_nmod_mat[]
  if isempty(lf)
    lf=composition_factors(M)
  end
  if length(lf)==1 && lf[1][2]==1
    return fq_nmod_mat[]
  end
  if dim!=n+1
    lf=[x for x in lf if x[1].dim==dim]
  end
  if isempty(lf)
    return list
  end
  G=M.G
  for x in lf
    append!(list,Hecke._irrsubs(x[1],M))
  end
  return list
end

doc"""
***
    maximal_submodules(M::FqGModule)

> Given a $G$-module $M$, it returns all the maximal submodules of M

"""

function maximal_submodules(M::FqGModule, index::Int=M.dim, lf=[])

  M_dual=dual_space(M)
  minlist=minimal_submodules(M_dual, index+1, lf)
  maxlist=Array{fq_nmod_mat,1}(length(minlist))
  for j=1:length(minlist)
    maxlist[j]=transpose(nullspace(minlist[j])[2])
  end
  return maxlist

end

doc"""
***
    submodules(M::FqGModule)

> Given a $G$-module $M$, it returns all the submodules of M

"""

function submodules(M::FqGModule)

  K=M.K
  list=fq_nmod_mat[]
  lf=composition_factors(M)
  minlist=minimal_submodules(M, M.dim+1, lf)
  for x in minlist
    rref!(x)
    N, pivotindex =actquo(x,M.G)
    ls=submodules(N)
    for a in ls
      s=zero_matrix(K,rows(a), M.dim)
      for t=1:rows(a)
        pos=0
        for j=1:M.dim
          if j in pivotindex
            pos+=1
          else
            s[t,j]=a[t,j-pos]
          end
        end
      end
      push!(list,vcat(x,s))
    end
  end
  for x in list
    rref!(x)
  end
  i=2
  while i<length(list)
    j=i+1
    while j<=length(list)
      if rows(list[j])!=rows(list[i])
        j+=1
      elseif list[j]==list[i]
        deleteat!(list, j)
      else
        j+=1
      end
    end
    i+=1
  end
  return append!(list,minlist)

end

doc"""
***
    submodules(M::FqGModule, index::Int)

> Given a $G$-module $M$, it returns all the submodules of M of index q^index, where q is the order of the field

"""

function submodules(M::FqGModule, index::Int; comp_factors=Tuple{FqGModule, Int}[])

  K=M.K
  if index==M.dim
    return fq_nmod_mat[zero_matrix(K,1,M.dim)]
  end
  list=fq_nmod_mat[]
  if index>= M.dim/2
    if index== M.dim -1
      if isempty(comp_factors)
        lf=composition_factors(M, dimension=1)
      else
        lf=comp_factors
      end
      list=minimal_submodules(M,1,lf)
      return list
    end
    if isempty(comp_factors)
      lf=composition_factors(M)
    else
      lf=comp_factors
    end
    for i=1: M.dim-index-1
      minlist=minimal_submodules(M,i,lf)
      for x in minlist
        N, pivotindex= actquo(x, M.G)
        #
        #  Recover the composition factors of the quotient
        #
        Sub=actsub(x, M.G)
        lf1=[(x[1], x[2]) for x in lf]
        for j=1:length(lf1)
          if isisomorphic(lf1[j][1], Sub)
            if lf1[j][2]==1
              deleteat!(lf1,j)
            else
              lf1[j]=(lf1[j][1], lf1[j][2]-1)
            end
            break
          end
        end
        #
        #  Recursively ask for submodules and write their bases in terms of the given set of generators
        #
        ls=submodules(N,index, comp_factors=lf1)
        for a in ls
          s=zero_matrix(K,rows(a)+rows(x), M.dim)
          for t=1:rows(a)
            pos=0
            for j=1:M.dim
              if j in pivotindex
               pos+=1
             else
               s[t,j]=a[t,j-pos]
              end
            end
          end
          for t=rows(a)+1:rows(s)
            for j=1:cols(s)
              s[t,j]=x[t-rows(a),j]
            end
          end
          push!(list,s)
        end
      end
    end

  #
  #  Eliminating repeatitions
  #

    for x in list
      rref!(x)
    end
    i=1
    while i<=length(list)
      k=i+1
      while k<=length(list)
        if list[i]==list[k]
          deleteat!(list, k)
        else
          k+=1
        end
      end
      i+=1
    end
    append!(list,minimal_submodules(M,M.dim-index, lf))
  else
  #
  #  Duality
  #
    M_dual=dual_space(M)
    dlist=submodules(M_dual, M.dim-index)
    list=fq_nmod_mat[transpose(nullspace(x)[2]) for x in dlist]
  end
  return list

end


text:
# ---
# title: 46. Permutations
# id: problem46
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Backtracking
# link: <https://leetcode.com/problems/permutations/description/>
# hidden: true
# ---
#
# Given an array `nums` of distinct integers, return _all the possible
# permutations_. You can return the answer in **any order**.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,2,3]
#     Output: [[1,2,3],[1,3,2],[2,1,3],[2,3,1],[3,1,2],[3,2,1]]
#
#
# **Example 2:**
#
#
#
#     Input: nums = [0,1]
#     Output: [[0,1],[1,0]]
#
#
# **Example 3:**
#
#
#
#     Input: nums = [1]
#     Output: [[1]]
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 6`
#   * `-10 <= nums[i] <= 10`
#   * All the integers of `nums` are **unique**.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end

ext: jl


text:
using OrdinaryDiffEq, CUDA, Test
CUDA.allowscalar(false)

# https://github.com/SciML/OrdinaryDiffEq.jl/issues/1614
function f(du, u, p, t)
    @. du = u
end

problem = ODEProblem(f, CUDA.ones(1), (0.0, 1.0))
stiffalg = Rosenbrock23()
for alg in (AutoDP5, AutoTsit5, AutoVern6, AutoVern7, AutoVern8, AutoVern9)
    @test_nowarn solve(problem, alg(stiffalg), save_everystep=false, save_start=false)
end


text:
struct SmallSignalOutput
    reduced_jacobian::AbstractArray{Float64}
    eigenvalues::Vector{Complex{Float64}}
    eigenvectors::Matrix{Complex{Float64}}
    index::Dict
    stable::Bool
    operating_point::Vector{Float64}
    damping::Dict{String, Dict{Symbol, Float64}}
    participation_factors::Dict{String, Dict{Symbol, Array{Float64}}}
end

function _determine_stability(vals::Vector{Complex{Float64}})
    for real_eig in real(vals)
        real_eig > 0.0 && return false
    end
    return true
end

function _calculate_forwardiff_jacobian(
    sim::Simulation{ResidualModel},
    x_eval::Vector{Float64},
)
    var_count = get_variable_count(sim.inputs)
    dx0 = zeros(var_count) #Define a vector of zeros for the derivative
    sysf! = (out, x) -> system_implicit!(
        out,            #output of the function
        dx0,            #derivatives equal to zero
        x,              #states
        sim.inputs,     #Parameters
        0.0,            #time equals to zero.
    )
    out = zeros(var_count) #Define a vector of zeros for the output
    jacobian = ForwardDiff.jacobian(sysf!, out, x_eval)
    return jacobian
end

function _calculate_forwardiff_jacobian(
    sim::Simulation{MassMatrixModel},
    x_eval::Vector{Float64},
)
    var_count = get_variable_count(sim.inputs)
    sysf! = (dx, x) -> system_mass_matrix!(
        dx,            #derivatives equal to zero
        x,              #states
        sim.inputs,     #Parameters
        0.0,            #time equals to zero.
    )
    dx = zeros(var_count) #Define a vector of zeros for the output
    jacobian = ForwardDiff.jacobian(sysf!, dx, x_eval)
    return jacobian
end

"""
Finds the location of the differential states in the reduced Jacobian
"""
function _make_reduced_jacobian_index(global_index, diff_states)
    jac_index = Dict{String, Dict{Symbol, Int}}()
    for (device_name, device_index) in global_index
        jac_index[device_name] = Dict{Symbol, Int}()
        for (state, ix) in device_index
            state_is_differential = diff_states[ix]
            if state_is_differential
                jac_index[device_name][state] = sum(diff_states[1:ix])
            end
        end
    end
    return jac_index
end

function _reduce_jacobian(
    jacobian::Matrix{Float64},
    diff_states::Vector{Bool},
    mass_matrix::LinearAlgebra.Diagonal{Float64, Vector{Float64}},
    global_index::MAPPING_DICT,
)
    alg_states = .!diff_states
    fx = @view jacobian[diff_states, diff_states]
    gy_temp = jacobian[alg_states, alg_states]
    ixs = range(1, length = length(alg_states))
    # Locate indices where the algebraic variables are
    alg_ix = ixs[alg_states]
    # Update alg_changes to false if all elements on both row and columns are zeros
    alg_changes = trues(length(alg_ix))
    for (ix, row) in enumerate(eachrow(gy_temp))
        row_changes = !all(x -> x == 0.0, row)
        col_changes = !all(x -> x == 0.0, gy_temp[:, ix])
        alg_changes[ix] = col_changes | row_changes
    end
    if !all(alg_changes)
        # Update algebraic states
        alg_states[alg_ix] .= alg_changes
        for b in alg_ix[.!alg_changes]
            name, state = get_state_from_ix(global_index, b)
            @info "Algebraic state $(state) in device $(name) was removed from the reduced jacobian due to having only zeros in both rows and columns."
        end
    end
    gy = jacobian[alg_states, alg_states]
    fy = @view jacobian[diff_states, alg_states]
    gx = @view jacobian[alg_states, diff_states]
    M = @view mass_matrix[diff_states, diff_states]
    inv_diag_M = 1.0 ./ LinearAlgebra.diag(M)
    inv_g = inv(gy)
    reduced_jacobian = inv_diag_M .* (fx - fy * inv_g * gx)
    return reduced_jacobian
end

function _get_eigenvalues(reduced_jacobian::AbstractArray{Float64}, multimachine::Bool)
    eigen_vals, R_eigen_vect = LinearAlgebra.eigen(Matrix(reduced_jacobian))
    if multimachine
        @warn(
            "No Infinite Bus found. Confirm stability directly checking eigenvalues.\nIf all eigenvalues are on the left-half plane and only one eigenvalue is zero, the system is small signal stable."
        )
        @debug(eigen_vals)
    end
    return eigen_vals, R_eigen_vect
end

function _get_damping(
    eigen_vals::Vector{Complex{Float64}},
    jac_index::Dict{String, Dict{Symbol, Int}},
)
    damping_results = Dict{String, Dict{Symbol, Float64}}()
    for (device_name, device_index) in jac_index
        damping_results[device_name] = Dict{Symbol, Float64}()
        for (state, ix) in device_index
            isnothing(ix) && continue
            eigen_val = eigen_vals[ix]
            damping_results[device_name][state] =
                -1 * real(eigen_val) / sqrt(real(eigen_val)^2 + imag(eigen_val)^2)
        end
    end
    return damping_results
end

function _get_participation_factors(
    R_eigen_vect::Matrix{Complex{Float64}},
    jac_index::Dict{String, Dict{Symbol, Int}},
)
    L_eigen_vect = inv(R_eigen_vect)
    participation_factors = Dict{String, Dict{Symbol, Array{Float64}}}()
    for (device_name, device_index) in jac_index
        participation_factors[device_name] = Dict{Symbol, Array{Float64}}()
        for (state, ix) in device_index
            den = sum(abs.(L_eigen_vect[:, ix]) .* abs.(R_eigen_vect[ix, :]))
            participation_factors[device_name][state] =
                abs.(L_eigen_vect[:, ix]) .* abs.(R_eigen_vect[ix, :]) ./ den
        end
    end
    return participation_factors
end

function _small_signal_analysis(
    jacobian::Matrix{Float64},
    x_eval::Vector{Float64},
    inputs::SimulationInputs,
    multimachine::Bool,
)
    mass_matrix = get_mass_matrix(inputs)
    diff_states = get_DAE_vector(inputs)
    global_index = make_global_state_map(inputs)
    jac_index = _make_reduced_jacobian_index(global_index, diff_states)
    reduced_jacobian = _reduce_jacobian(jacobian, diff_states, mass_matrix, global_index)
    eigen_vals, R_eigen_vect = _get_eigenvalues(reduced_jacobian, multimachine)
    damping = _get_damping(eigen_vals, jac_index)
    stable = _determine_stability(eigen_vals)
    participation_factors = _get_participation_factors(R_eigen_vect, jac_index)
    return SmallSignalOutput(
        reduced_jacobian,
        eigen_vals,
        R_eigen_vect,
        jac_index,
        stable,
        x_eval,
        damping,
        participation_factors,
    )
end

function _small_signal_analysis(
    ::Type{T},
    inputs::SimulationInputs,
    x_eval::Vector{Float64},
    multimachine = true,
) where {T <: SimulationModel}
    jacwrapper = get_jacobian(T, inputs, x_eval, 0)
    return _small_signal_analysis(jacwrapper.Jv, jacwrapper.x, inputs, multimachine)
end

function small_signal_analysis(sim::Simulation{T}; kwargs...) where {T <: SimulationModel}
    inputs = get_simulation_inputs(sim)
    x_eval = get(kwargs, :operating_point, get_initial_conditions(sim))
    return _small_signal_analysis(T, inputs, x_eval, sim.multimachine)
end

function small_signal_analysis(::Type{T}, system::PSY.System) where {T <: SimulationModel}
    simulation_system = deepcopy(system)
    inputs = SimulationInputs(T, simulation_system, ReferenceBus)
    x0_init = get_flat_start(inputs)
    set_operating_point!(x0_init, inputs, system)
    return _small_signal_analysis(T, inputs, x0_init)
end


text:
using Test

if VERSION >= v"1.2.0-DEV.573"
    include("snoopi.jl")
end

if VERSION >= v"1.6.0-DEV.1190"  # https://github.com/JuliaLang/julia/pull/37749
    @testset "snoopi_deep" begin
        include("snoopi_deep.jl")
    end
end

if VERSION >= v"1.6.0-DEV.1192"  # https://github.com/JuliaLang/julia/pull/37136
    @testset "snoopl" begin
        include("snoopl.jl")
    end
end

using SnoopCompile

# issue #26
logfile = joinpath(tempdir(), "anon.log")
@snoopc logfile begin
    map(x->x^2, [1,2,3])
end
data = SnoopCompile.read(logfile)
pc = SnoopCompile.parcel(reverse!(data[2]))
@test length(pc[:Base]) <= 1

# issue #29
keep, pcstring, topmod, name = SnoopCompile.parse_call("Tuple{getfield(JLD, Symbol(\"##s27#8\")), Any, Any, Any, Any, Any}")
@test keep
@test pcstring == "Tuple{getfield(JLD, Symbol(\"##s27#8\")), Int, Int, Int, Int, Int}"
@test topmod == :JLD
@test name == "##s27#8"
logfile = joinpath(tempdir(), "isdefined.log")
@snoopc logfile begin
    @eval module IsDef
        @generated function tm(x)
            return :(typemax($x))
        end
    end
    IsDef.tm(1)
end
data = SnoopCompile.read(logfile)
pc = SnoopCompile.parcel(reverse!(data[2]))
@test any(startswith.(pc[:IsDef], "isdefined"))

#=
# Simple call
let str = "sum"
    keep, pcstring, topmod = SnoopCompile.parse_call("Foo.any($str)")
    @test keep
    @test pcstring == "Tuple{$str}"
    @test topmod == :Main
end

# Operator
let str = "Base.:*, Int, Int"
    keep, pcstring, topmod = SnoopCompile.parse_call("Foo.any($str)")
    @test keep
    @test pcstring == "Tuple{$str}"
    @test topmod == :Base
end

# Function as argument
let str = "typeof(Base.identity), Array{Bool, 1}"
    keep, pcstring, topmod = SnoopCompile.parse_call("Foo.any($str, Vararg{Any, N} where N)")
    @test keep
    @test pcstring == "Tuple{$str, Int}"
    @test topmod == :Base
end

# Anonymous function closure in a new module as argument
let func = (@eval Main module SnoopTestTemp
            func = () -> (y = 2; (x -> x > y))
        end).func
    str = "getfield(SnoopTestTemp, Symbol(\"$(typeof(func()))\")), Array{Float32, 1}"
    keep, pcstring, topmod = SnoopCompile.parse_call("Foo.any($str)")
    @test keep
    @test pcstring == "Tuple{$str}"
    @test topmod == :SnoopTestTemp
end

# Function as a type
let str = "typeof(Base.Sort.sort!), Array{Any, 1}, Base.Sort.MergeSortAlg, Base.Order.By{typeof(Base.string)}"
    keep, pcstring, topmod = SnoopCompile.parse_call("Foo.Bar.sort!($str)")
    @test keep
    @test pcstring == "Tuple{$str}"
    @test topmod == :Base
end
=#

include("colortypes.jl")

if isdefined(SnoopCompile, :invalidation_trees)
    include("snoopr.jl")
end


text:
t_using = (tstart = time(); using CairoMakie; time() - tstart)

t_plot = (tstart = time(); save("test.png", scatter(1:4)); time() - tstart)
rm("test.png")
println("($t_using, $t_plot)")


text:
using TenhouLogs
const Mj = TenhouLogs

function Mj.analyzer(::Val{Mj.matchend}, pst::Mj.PlayState)
    @show pst.scores
end

Mj.analyseDatabase("dataset\\scraw2019s4p.db";
        offset = 170000, limit = 0)


text:
using Documenter, Modia3D

makedocs(
  modules  = [Modia3D],
  sitename = "Modia3D",
  authors  = "Andrea Neumayr, Martin Otter (DLR-SR)",
  format   = Documenter.HTML(prettyurls = false),
  pages    = [
     "Home"   => "index.md",
     "Manual" => [
        "man/GettingStarted.md"
        "man/Examples.md"
        "man/Materials.md"
        "man/CollisionHandling.md"
        "man/Plans.md"
        ],
     "Library" => [
        "lib/Composition.md",
        "lib/Graphics.md",
        "lib/Solids.md",
        "lib/ForceElements.md",
        "lib/Basics.md"
     ],
     "Internal" => [
        "internal/ContactDetection.md"
     ]
  ]
)


text:
#################### Missing Values Sampler ####################

#################### Types and Constructors ####################

struct MISSTune <: SamplerTune
    logf::Function
    dims::Tuple
    valueinds::AbstractArray
    distrinds::AbstractArray

end

function MISSTune(f::Function)
    MISSTune(f, (), Int[], Int[])
end

function MISSTune(s::AbstractStochastic, gm::Function)
    MISSTune(s.distr, s.value, gm)
end

function MISSTune(d::Distribution, v, gm::Function)
    MISSTune(gm, (), findall(isnan.(v)), Int[])
end

function MISSTune(D::Array{UnivariateDistribution}, v::Array, gm::Function)
    inds = findall(isnan.(v))
    MISSTune(gm, dims(D), inds, inds)
end

function MISSTune(D::Array{MultivariateDistribution}, v::Array, gm::Function)
    isvalueinds = falses(size(v))
    isdistrinds = falses(size(D))
    for sub in CartesianIndices(size(D))
        n = length(D[sub])
        for i = 1:n
            if isnan(v[sub, i])
                isvalueinds[sub, i] = isdistrinds[sub] = true
            end
        end
    end
    MISSTune(gm, dims(D), findall(isvalueinds), findall(isdistrinds))
end

const MISSVariate = Sampler{MISSTune,T} where {T}

#################### Sampler Constructor ####################
"""
    MISS(params::ElementOrVector{Symbol})

Construct a `Sampler` object to sampling missing output values. The constructor
should only be used to sample stochastic nodes upon which no other stochastic
node depends. So-called ‘output nodes’ can be identified with the `keys()`
function. Moreover, when the `MISS` constructor is included in a vector of
`Sampler` objects to define a sampling scheme, it should be positioned at the
beginning of the vector. This ensures that missing output values are updated
before any other samplers are executed.

Returns a `Sampler{Dict{Symbol, MISSTune}}` type object.

* `params`: stochastic node(s) that contain missing values (`NaN`) to be updated with the sampler.
"""
function MISS(params::Symbol)
    lf(m, args...) = m
    params = asvec(params)

    Sampler(Float64[], params, MISSTune(lf), Symbol[], false)
end

#################### Sampling Functions ####################

function sample!(v::MISSVariate{T}, get_model::Function; gen::Int, model::Model, kwargs...) where {T}
    params = v.params
    if gen == 1
        for key in params
            v.tune = MISSTune(model[key], get_model)
        end
    end
    for key in params
        node = model[key]

        node[v.tune.valueinds] = rand(node, v.tune)
        update!(model, node.targets)
    end
    nothing
end

rand(s::AbstractStochastic, miss::MISSTune) = rand_sub(s.distr, miss)

function rand_sub(d::Distribution, miss::MISSTune)
    x = rand(d)
    Float64[x[i] for i in miss.valueinds]
end

function rand_sub(D::Array{UnivariateDistribution}, miss::MISSTune)
    Float64[rand(d) for d in D[miss.distrinds]]
end

function rand_sub(D::Array{MultivariateDistribution}, miss::MISSTune)
    X = Array{Float64}(undef, miss.dims)
    for i in miss.distrinds
        d = D[i]
        X[ind2sub(D, i)..., 1:length(d)] = rand(d)
    end
    X[miss.valueinds]
end

max_stars_repo_name: erathorn/MCPhylo

max_issues_repo_name: erathorn/MCPhylo

max_forks_repo_name: erathorn/MCPhylo


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule ModularForms_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("ModularForms")
JLLWrappers.@generate_main_file("ModularForms", UUID("af2f1493-5da4-5c9a-92e4-d8b087733894"))
end  # module ModularForms_jll


text:
const ANG_THD = 1e-2

mutable struct Connection
    source::Int
    dest::Int
    path::Int
    laneConnections::Array{Tuple{Int,Int}}
end

function Connection(source::Int,dest::Int)
    Connection(source,dest,0,Array{Tuple{Int,Int}}(undef, 0))
end

mutable struct Junction
    connections::Array{Connection}
end

function get_min_difference_angles(A::Float64,B::Float64)
    b=B
    as=[A A+2*pi A-2*pi];
    (abs_diff,ind_a)=findmin(abs.(as .- b))
    a=as[ind_a]
    return a,b,abs_diff
end

function get_min_difference_angles_pos(A::Float64,B::Float64)
    b = B
    as = [A A+2*pi A-2*pi];
    diffs = b .- as
    pos_inds = findall(x->x>=0,diffs)
    pos_as = as[pos_inds]
    pos_diffs = diffs[pos_inds]
    pos_diff,ind_a = findmin(abs.(pos_diffs))
    a = pos_as[ind_a]
    return a,b,pos_diff
end

function get_min_difference_angles_neg(A::Float64,B::Float64)
    b = B
    as = [A A+2*pi A-2*pi];
    diffs = b .- as
    neg_inds = findall(x->x<=0,diffs)
    neg_as = as[neg_inds]
    neg_diffs = diffs[neg_inds]
    neg_diff,ind_a = findmin(abs.(neg_diffs))
    neg_diff = - neg_diff
    a = neg_as[ind_a]
    return a,b,neg_diff
end

function get_intersection_point(A::VecSE2, B::VecSE2)
    x1 = A.x
    y1 = A.y
    x2 = B.x
    y2 = B.y
    θ1,θ2,_ = get_min_difference_angles(A.θ,B.θ)
    a1 = tan(θ1)
    a2 = tan(θ2)
    @assert a1 != a2
    @assert abs(a1) < 1e4 || abs(a2) < 1e4
    if abs(a1) < 1e4 && abs(a2) < 1e4
            x = -(-a1*x1+a2*x2+y1-y2)/(a1-a2)
            y = -(-a1*a2*x1+a1*a2*x2+a2*y1-a1*y2)/(a1-a2)
    else
        if abs(a1) > 1e4
            x = x1
            y = a2*x1-a2*x2+y2
        elseif abs(a2) > 1e4
            x = x2
            y = -a1*x1+a1*x2+y1
        end
    end
    #println("A ",A)
    #println("B ",B)
    #println("Intersection ",VecE2(x,y))
    return VecE2(x,y)
end

function get_center_of_circle(A::VecSE2, B::VecSE2)
    x1 = A.x
    y1 = A.y
    x2 = B.x
    y2 = B.y
    θ1,θ2,angle_diff = get_min_difference_angles(A.θ,B.θ)
    x = -((-x2*cos(θ2)*sin(θ1)+x1*cos(θ1)*sin(θ2)+y1*sin(θ1)*sin(θ2)-y2*sin(θ1)*sin(θ2))/(cos(θ2)*sin(θ1)-cos(θ1)*sin(θ2)))
    y = -((-x1*cos(θ1)*cos(θ2)+x2*cos(θ1)*cos(θ2)-y1*cos(θ2)*sin(θ1)+y2*cos(θ1)*sin(θ2))/(cos(θ2)*sin(θ1)-cos(θ1)*sin(θ2)))

    if abs(abs(angle_diff)-pi) < ANG_THD
        x = (x1+x2)/2
        y = (y1+y2)/2
    end
    return VecE2(x,y)
end

function connect_two_points_by_circle(A::VecSE2{T}, B::VecSE2{T};s_base::T=zero(T)) where T<:Real
    ncurvepts_per_turn = 25

    distance = norm(VecE2(A-B))

    C = get_center_of_circle(A, B)
    CA = VecE2(A.x-C.x,A.y-C.y)
    CB = VecE2(B.x-C.x,B.y-C.y)

    if (CA.x>ANG_THD && sin(A.θ)>ANG_THD) || (CA.x<-ANG_THD && sin(A.θ)<-ANG_THD) || (CA.y<-ANG_THD && cos(A.θ)>ANG_THD) || (CA.y>ANG_THD && cos(A.θ)<-ANG_THD)
        attangle_a,attangle_b,angle_diff = get_min_difference_angles_pos(A.θ,B.θ)
        posangle_a,posangle_b,angle_diff = get_min_difference_angles_pos(atan(CA.y,CA.x),atan(CB.y,CB.x))
    else
        attangle_a,attangle_b,angle_diff = get_min_difference_angles_neg(A.θ,B.θ)
        posangle_a,posangle_b,angle_diff = get_min_difference_angles_neg(atan(CA.y,CA.x),atan(CB.y,CB.x))
    end

    curvepts = Array{CurvePt{T}}(undef, ncurvepts_per_turn)

    #@assert abs(angle_diff) > ANG_THD
    if abs(abs(angle_diff)-pi) < ANG_THD
        radius = distance/2
    else
        radius = sqrt(distance^2/(2*(1-cos(angle_diff))))
    end
    r = radius
    for j in 1:ncurvepts_per_turn
        t = (j-1)/(ncurvepts_per_turn-1) # ∈ [0,1]
        s = r*abs(angle_diff)*t+s_base
        k = 1.0/r
        if attangle_b < attangle_a
            k = -k
        end
        kd = 0.0
        curvepts[j] = CurvePt(VecSE2(C + polar(r, lerp(posangle_a,posangle_b, t)), lerp(attangle_a,attangle_b,t)),  s,k,kd)
    end
    total_length = r*abs(angle_diff)+s_base
    return curvepts,total_length
end

function connect_two_points(A::VecSE2{T}, B::VecSE2{T}) where T<:Real
    #connect any two points with one points array
    _,_,abs_diff = get_min_difference_angles(B.θ,A.θ)

    curvepts = Array{CurvePt{T}}(undef, 0)
    push!(curvepts,CurvePt(A,  0.0,0.0,0.0))
    if abs_diff < ANG_THD
        push!(curvepts,CurvePt(B,  norm(VecE2(A-B)),0.0,0.0))
    else
        if abs(abs_diff - pi) < ANG_THD
            #println("difference is pi")
            BA = VecE2(B.x-A.x,B.y-A.y)
            pojection_length = BA.x*cos(A.θ)+BA.y*sin(A.θ)
            if pojection_length > ANG_THD
                P = VecSE2(VecE2(A + polar(pojection_length,A.θ)),A.θ)
                push!(curvepts,CurvePt(P, pojection_length,0.0,0.0))
                circlepts,total_length = connect_two_points_by_circle(P, B;s_base = pojection_length)
                curvepts = vcat(curvepts,circlepts)
                push!(curvepts,CurvePt(B,total_length,0.0,0.0))
            elseif pojection_length < -ANG_THD
                P = VecSE2(VecE2(B + polar(pojection_length,B.θ)),B.θ)
                circlepts,total_length = connect_two_points_by_circle(A, P;s_base = 0.0)
                curvepts = vcat(curvepts,circlepts)
                push!(curvepts,CurvePt(P, total_length,0.0,0.0))
                push!(curvepts,CurvePt(B,total_length+abs(pojection_length),0.0,0.0))
            else #no need of P
                circlepts,total_length = connect_two_points_by_circle(A, B;s_base = 0.0)
                curvepts = vcat(curvepts,circlepts)
            end
        else
            #println("difference is general")
            I = get_intersection_point(A, B)
            dis1 = norm(VecE2(A-I))
            dis2 = norm(VecE2(B-I))
            AI = VecE2(I.x-A.x,I.y-A.y)
            BI = VecE2(I.x-B.x,I.y-B.y)
            dis_diff = abs(dis1-dis2)
            projection_length = AI.x*cos(A.θ)+ AI.y*sin(A.θ)
            projection_lengthB = BI.x*cos(B.θ)+ BI.y*sin(B.θ)
            if projection_length>=0 && projection_lengthB>=0
                #println("0")
                P = VecSE2(I+polar(projection_lengthB,A.θ),A.θ)
                push!(curvepts,CurvePt(P, projection_length+projection_lengthB,0.0,0.0))
                circlepts,total_length = connect_two_points_by_circle(P, B;s_base = projection_length+projection_lengthB)
                curvepts = vcat(curvepts,circlepts)
                push!(curvepts,CurvePt(B,total_length,0.0,0.0))
            elseif (dis1-dis2 >= 0.0 && projection_length >= 0.0) || (dis1-dis2 <= 0.0 && projection_length <= 0.0)
                #println("1")
                P = VecSE2(VecE2(A + polar(dis_diff,A.θ)),A.θ)
                push!(curvepts,CurvePt(P, dis_diff,0.0,0.0))
                circlepts,total_length = connect_two_points_by_circle(P, B;s_base = dis_diff)
                curvepts = vcat(curvepts,circlepts)
                push!(curvepts,CurvePt(B,total_length,0.0,0.0))
            elseif (dis1-dis2 >= 0.0 && projection_length <= 0.0) || (dis1-dis2 <= 0.0 && projection_length >= 0.0)
                #println("2")
                P = VecSE2(VecE2(B + polar(-dis_diff,B.θ)),B.θ)
                circlepts,total_length = connect_two_points_by_circle(A, P;s_base = 0.0)
                curvepts = vcat(curvepts,circlepts)
                push!(curvepts,CurvePt(P, total_length,0.0,0.0))
                push!(curvepts,CurvePt(B,total_length+dis_diff,0.0,0.0))
            #else
            #    println("3")
            #    circlepts,total_length = connect_two_points_by_circle(A, B;s_base = 0.0)
            #    curvepts = vcat(curvepts,circlepts)
            end
        end
    end
    return curvepts
end

function connect_two_points_seperate(A::VecSE2{T}, B::VecSE2{T}) where T<:Real
    #connect any two points with an array of curpve points
    _,_,abs_diff = get_min_difference_angles(B.θ,A.θ)

    curvepts = Array{CurvePt{T}}[]
    if abs_diff < ANG_THD
        push!(curvepts,[CurvePt(A,  0.0,0.0,0.0)])
        push!(curvepts[1],CurvePt(B,  norm(VecE2(A-B)),0.0,0.0))
    else
        if abs(abs_diff - pi) < ANG_THD
            #println("difference is pi")
            BA = VecE2(B.x-A.x,B.y-A.y)
            pojection_length = BA.x*cos(A.θ)+BA.y*sin(A.θ)
            if pojection_length > ANG_THD
                push!(curvepts,[CurvePt(A,  0.0,0.0,0.0)])
                P = VecSE2(VecE2(A + polar(pojection_length,A.θ)),A.θ)
                push!(curvepts[1],CurvePt(P, pojection_length,0.0,0.0))
                circlepts,total_length = connect_two_points_by_circle(P, B;s_base = 0.0)
                push!(curvepts,circlepts)
                push!(curvepts[2],CurvePt(B,total_length,0.0,0.0))
            elseif pojection_length < -ANG_THD
                P = VecSE2(VecE2(B + polar(pojection_length,B.θ)),B.θ)
                circlepts,total_length = connect_two_points_by_circle(A, P;s_base = 0.0)
                push!(curvepts,circlepts)
                push!(curvepts,[CurvePt(P, 0.0,0.0,0.0)])
                push!(curvepts[2],CurvePt(B,abs(pojection_length),0.0,0.0))
            else #no need of P
                circlepts,total_length = connect_two_points_by_circle(A, B;s_base = 0.0)
                push!(curvepts,circlepts)
            end
        else
            #println("difference is general")
            I = get_intersection_point(A, B)
            dis1 = norm(VecE2(A-I))
            dis2 = norm(VecE2(B-I))
            AI = VecE2(I.x-A.x,I.y-A.y)
            BI = VecE2(I.x-B.x,I.y-B.y)
            dis_diff = abs(dis1-dis2)
            projection_length = AI.x*cos(A.θ)+ AI.y*sin(A.θ)
            projection_lengthB = BI.x*cos(B.θ)+ BI.y*sin(B.θ)
            if dis_diff < ANG_THD
                circlepts,total_length = connect_two_points_by_circle(A, B;s_base = 0.0)
                push!(curvepts,circlepts)
            else
                if projection_length>=0 && projection_lengthB>=0
                    #println("0")
                    push!(curvepts,[CurvePt(A,  0.0,0.0,0.0)])
                    P = VecSE2(I+polar(projection_lengthB,A.θ),A.θ)
                    push!(curvepts[1],CurvePt(P, projection_length+projection_lengthB,0.0,0.0))
                    circlepts,total_length = connect_two_points_by_circle(P, B;s_base = 0.0)
                    push!(curvepts,circlepts)
                elseif (dis1-dis2 >= 0.0 && projection_length >= 0.0) || (dis1-dis2 <= 0.0 && projection_length <= 0.0)
                    #println("1")
                    push!(curvepts,[CurvePt(A,  0.0,0.0,0.0)])
                    P = VecSE2(VecE2(A + polar(dis_diff,A.θ)),A.θ)
                    push!(curvepts[1],CurvePt(P, dis_diff,0.0,0.0))
                    circlepts,total_length = connect_two_points_by_circle(P, B;s_base = 0.0)
                    push!(curvepts,circlepts)
                elseif (dis1-dis2 >= 0.0 && projection_length <= 0.0) || (dis1-dis2 <= 0.0 && projection_length >= 0.0)
                    #println("2")
                    P = VecSE2(VecE2(B + polar(-dis_diff,B.θ)),B.θ)
                    circlepts,total_length = connect_two_points_by_circle(A, P;s_base = 0.0)
                    push!(curvepts,circlepts)
                    push!(curvepts,[CurvePt(P, 0.0,0.0,0.0)])
                    push!(curvepts[2],CurvePt(B,dis_diff,0.0,0.0))
                end
            end
        end
    end
    return curvepts
end

function connect_two_lane_general!(source::Lane, dest::Lane,seg_id::Int,lane_width::Float64;
    boundary_left::LaneBoundary=LaneBoundary(:solid, :white),
    boundary_right::LaneBoundary=LaneBoundary(:solid, :white),
    boundary_middle::LaneBoundary=LaneBoundary(:solid, :white)
    )
    cindS = curveindex_end(source.curve)
    cindD = CURVEINDEX_START
    A = source.curve[cindS].pos
    B = dest.curve[cindD].pos
    seg12 = RoadSegment(seg_id, Array{Lane}(undef, 1))
    curvepts = connect_two_points(A, B)
    tag12 = LaneTag(seg_id,1)

    seg12.lanes[1] = Lane(tag12, curvepts, width=lane_width,
                                        boundary_left=boundary_left, boundary_right=boundary_right,
                                        )
    connect!(source, seg12.lanes[1])
    connect!(seg12.lanes[1], dest)
    return seg12
end

function connect_two_lane_general!(source::Lane, dest::Lane,roadway::Roadway,lane_width::Float64;
    boundary_left::LaneBoundary=LaneBoundary(:solid, :white),
    boundary_right::LaneBoundary=LaneBoundary(:solid, :white),
    boundary_middle::LaneBoundary=LaneBoundary(:solid, :white)
    )
    cindS = curveindex_end(source.curve)
    cindD = CURVEINDEX_START
    A = source.curve[cindS].pos
    B = dest.curve[cindD].pos
    curveptss = connect_two_points_seperate(A::VecSE2, B::VecSE2)
    segs = []
    seg_id = length(roadway.segments)
    for curvepts in curveptss
        seg_id += 1
        seg12 = RoadSegment(seg_id, Array{Lane}(undef, 1))
        tag12 = LaneTag(seg_id,1)
        seg12.lanes[1] = Lane(tag12, curvepts, width=lane_width,
                                        boundary_left=boundary_left, boundary_right=boundary_right,
                                        )
        push!(segs,seg12)
    end

    before = source
    next = nothing
    for i = 1:length(segs)
        connect!(before,segs[i].lanes[1])
        if i == length(segs)
            next = dest
        else
            next = segs[i+1].lanes[1]
        end
        connect!(segs[i].lanes[1],next)
        before = segs[i].lanes[1]
    end
    for i = 1:length(segs)
        push!(roadway.segments, segs[i])
    end
end

function connect_two_seg!(source::RoadSegment{T},dest::RoadSegment{T},roadway::Roadway{T};
    connections::Array{Tuple{Int,Int}}=Array{Tuple{Int,Int}}(undef, 0),
    boundary_left::LaneBoundary=LaneBoundary(:solid, :white),
    boundary_right::LaneBoundary=LaneBoundary(:solid, :white),
    boundary_middle::LaneBoundary=LaneBoundary(:solid, :white)
    ) where T
    if isempty(connections)
        for i=1:length(source.lanes)
            push!(connections,(i,i))
        end
    end
    curveptsss = []
    for i = 1:length(connections)
        cindS = curveindex_end(source.lanes[connections[i][1]].curve)
        cindD = CURVEINDEX_START
        A = source.lanes[connections[i][1]].curve[cindS].pos
        B = dest.lanes[connections[i][2]].curve[cindD].pos
        curveptss = connect_two_points_seperate(A, B)
        push!(curveptsss,curveptss)
    end
    segs = []
    seg_id = length(roadway.segments)

    for i = 1:length(curveptsss[1])
        seg_id += 1
        seg12 = RoadSegment{T}(seg_id, Vector{Lane{T}}(undef, length(connections)))
        for j =  1:length(connections)
            tag12 = LaneTag(seg_id,j)
            seg12.lanes[j] = Lane(tag12, curveptsss[j][i], width=source.lanes[connections[j][1]].width,
                                        boundary_left=boundary_left, boundary_right=boundary_right,
                                        )
        end
        push!(segs,seg12)
    end

    for i = 1:length(connections)
        before = source.lanes[connections[i][1]]
        next = segs[1].lanes[i]
        connect!(before,next)
        before = next
        for j = 1:length(segs)
            if j == length(segs)
                next = dest.lanes[connections[i][2]]
            else
                next = segs[j+1].lanes[i]
            end
            connect!(before,next)
            before = next
        end
    end
    for i = 1:length(segs)
        push!(roadway.segments, segs[i])
    end

end

function connect_two_seg_general!(source::RoadSegment{T},dest::RoadSegment{T},roadway::Roadway{T};
    direction::Tuple{Int,Int}=(1,1),
    connections::Array{Tuple{Int,Int}}=Array{Tuple{Int,Int}}(undef, 0),
    boundary_left::LaneBoundary=LaneBoundary(:solid, :white),
    boundary_right::LaneBoundary=LaneBoundary(:solid, :white),
    boundary_middle::LaneBoundary=LaneBoundary(:solid, :white)
    ) where T<:Real
    if isempty(connections)
        for i=1:length(source.lanes)
            push!(connections,(i,i))
        end
    end
    curveptsss = []
    for i = 1:length(connections)
        cindS = curveindex_end(source.lanes[connections[i][1]].curve)
        if direction[1]==-1
            cindS = CURVEINDEX_START
        end
        print()
        cindD = CURVEINDEX_START
        if direction[2]==-1
            cindD = curveindex_end(dest.lanes[connections[i][2]].curve)
        end
        A = source.lanes[connections[i][1]].curve[cindS].pos
        if direction[1]==-1
            A = VecSE2(A.x,A.y,mod2pi(A.θ+pi))
        end
        print(A)
        B = dest.lanes[connections[i][2]].curve[cindD].pos
        if direction[2]==-1
            B = VecSE2(B.x,B.y,mod2pi(B.θ+pi))
        end
        print(B)
        curveptss = connect_two_points_seperate(A, B)
        push!(curveptsss,curveptss)
    end
    segs = []
    seg_id = length(roadway.segments)

    for i = 1:length(curveptsss[1])
        seg_id += 1
        seg12 = RoadSegment{T}(seg_id, Array{Lane{T}}(undef, length(connections)))
        for j =  1:length(connections)
            tag12 = LaneTag(seg_id,j)
            seg12.lanes[j] = Lane(tag12, curveptsss[j][i], width=source.lanes[connections[j][1]].width,
                                        boundary_left=boundary_left, boundary_right=boundary_right,
                                        )
        end
        push!(segs,seg12)
    end

    for i = 1:length(connections)
        before = source.lanes[connections[i][1]]
        next = segs[1].lanes[i]
        if direction[1] == 1
            connect_general!(before,next,(1,1))
        elseif direction[1] == -1
            connect_general!(before,next,(-1,1))
        end
        before = next
        for j = 1:length(segs)
            if j == length(segs)
                next = dest.lanes[connections[i][2]]
                if direction[2] == -1
                    connect_general!(before,next,(1,-1))
                elseif direction[1] == 1
                    connect_general!(before,next,(1,1))
                end
            else
                next = segs[j+1].lanes[i]
                connect_general!(before,next,(1,1))
            end
            before = next
        end
    end

    for i = 1:length(segs)
        push!(roadway.segments, segs[i])
    end

end

function connect_general!(source::Lane, dest::Lane,direction::Tuple{Int,Int}=(1,1))
    # place these at the front
    if direction[1]==1 && direction[2] ==1
        cindS = curveindex_end(source.curve)
        cindD = CURVEINDEX_START

        unshift!(source.exits,   LaneConnection(true,  cindS, RoadIndex(cindD, dest.tag)))
        unshift!(dest.entrances, LaneConnection(false, cindD, RoadIndex(cindS, source.tag)))
        (source, dest)
    elseif direction[1]==1 && direction[2] ==-1
        cindS = curveindex_end(source.curve)
        cindD = curveindex_end(dest.curve)

        unshift!(source.exits,   LaneConnection(true,  cindS, RoadIndex(cindD, dest.tag)))
        unshift!(dest.exits, LaneConnection(true, cindD, RoadIndex(cindS, source.tag)))
        (source, dest)
    elseif direction[1]==-1 && direction[2] ==1
        cindS = CURVEINDEX_START
        cindD = CURVEINDEX_START

        unshift!(source.entrances,   LaneConnection(false,  cindS, RoadIndex(cindD, dest.tag)))
        unshift!(dest.entrances, LaneConnection(false, cindD, RoadIndex(cindS, source.tag)))
        (source, dest)
    elseif direction[1]==-1 && direction[2] ==-1
        cindS = CURVEINDEX_START
        cindD = curveindex_end(dest.curve)

        unshift!(source.entrances,   LaneConnection(false,  cindS, RoadIndex(cindD, dest.tag)))
        unshift!(dest.exits, LaneConnection(true, cindD, RoadIndex(cindS, source.tag)))
        (source, dest)
    end
end

function add_connection!(connection::Connection,roadway::Roadway{T};
    boundary_left::LaneBoundary=LaneBoundary(:solid, :white),
    boundary_right::LaneBoundary=LaneBoundary(:solid, :white),
    boundary_middle::LaneBoundary=LaneBoundary(:solid, :white)
    ) where T
    source = roadway.segments[connection.source]
    dest = roadway.segments[connection.dest]
    if isempty(connection.laneConnections)
        for i=1:length(source.lanes)
            push!(connection.laneConnections,(i,i))
        end
    end
    curveptsss = []
    for i = 1:length(connection.laneConnections)
        cindS = curveindex_end(source.lanes[connection.laneConnections[i][1]].curve)
        cindD = CURVEINDEX_START
        A = source.lanes[connection.laneConnections[i][1]].curve[cindS].pos
        B = dest.lanes[connection.laneConnections[i][2]].curve[cindD].pos
        curveptss = connect_two_points_seperate(A, B)
        push!(curveptsss,curveptss)
    end
    segs = []
    seg_id = length(roadway.segments)

    for i = 1:length(curveptsss[1])
        seg_id += 1
        seg12 = RoadSegment{T}(seg_id, Vector{Lane{T}}(undef, length(connection.laneConnections)))
        if i==1
            connection.path = seg_id
        end
        for j =  1:length(connection.laneConnections)
            tag12 = LaneTag(seg_id,j)
            seg12.lanes[j] = Lane(tag12, curveptsss[j][i], width=source.lanes[connection.laneConnections[j][1]].width,
                                        boundary_left=boundary_left, boundary_right=boundary_right,
                                        )
        end
        push!(segs,seg12)
    end

    for i = 1:length(connection.laneConnections)
        before = source.lanes[connection.laneConnections[i][1]]
        next = nothing
        for j = 1:length(segs)
            connect!(before,segs[j].lanes[i])
            if j == length(segs)
                next = dest.lanes[connection.laneConnections[i][2]]
            else
                next = segs[j+1].lanes[i]
            end
            connect!(segs[j].lanes[i],next)
            before = segs[j].lanes[1]
        end
    end
    for i = 1:length(segs)
        push!(roadway.segments, segs[i])
    end

end

function add_junction!(junction::Junction,roadway::Roadway)
    for connection in junction.connections
        add_connection!(connection,roadway)
    end
end

function gen_connected_lanes(;nlanes::Int=1,roadlength::Float64=5.0,
    lane_width::Float64=DEFAULT_LANE_WIDTH, # [m]
    lane_widths::Vector{Float64} = fill(lane_width, nlanes),
    origin1::VecSE2{T} = VecSE2(0.0,0.0,0.0),
    origin2::VecSE2{T} = VecSE2(origin1.x+roadlength,origin1.y-roadlength,1.25*pi),
    boundary_leftmost::LaneBoundary=LaneBoundary(:solid, :white),
    boundary_rightmost::LaneBoundary=LaneBoundary(:solid, :white),
    boundary_middle::LaneBoundary=LaneBoundary(:broken, :white)
    ) where T<:Real

    segments = RoadSegment{T}[]
    seg1 = gen_straight_segment(1, nlanes, roadlength,
                                                origin=origin1, lane_widths=lane_widths,
                                                boundary_leftmost=boundary_leftmost,
                                                boundary_rightmost=boundary_rightmost,
                                                boundary_middle=boundary_middle)
    seg2 = gen_straight_segment(2, nlanes, roadlength,
                                                origin=origin2, lane_widths=lane_widths,
                                                boundary_leftmost=boundary_leftmost,
                                                boundary_rightmost=boundary_rightmost,
                                                boundary_middle=boundary_middle)

    retval = Roadway{T}()
    push!(retval.segments, seg1)
    push!(retval.segments, seg2)

    boundary_left =  boundary_middle
    boundary_right = boundary_middle

    connect_two_seg!(retval.segments[1],retval.segments[2],retval)
    retval
end

function gen_intersection(;nlanes::Int=2,roadlength::Float64=5.0,
    lane_width::Float64=DEFAULT_LANE_WIDTH, # [m]
    lane_widths::Vector{Float64} = fill(lane_width, nlanes),
    crosswidth::Float64=nlanes*lane_width,
    origin::VecSE2 = VecSE2(0.0,0.0,0.0),
    #connections::Array{Tuple{Int,Int}}=[(1,6),(3,8),(5,2),(7,4),(1,2),(3,4),(5,6),(7,8),(1,4),(3,6),(5,8),(7,2)],
    junction::Junction=Junction([Connection(1,6),Connection(3,8),Connection(5,2),Connection(7,4),Connection(1,2),Connection(3,4),Connection(5,6),Connection(7,8),Connection(1,4),Connection(3,6),Connection(5,8),Connection(7,2)]),
    boundary_leftmost::LaneBoundary=LaneBoundary(:solid, :white),
    boundary_rightmost::LaneBoundary=LaneBoundary(:solid, :white),
    boundary_middle::LaneBoundary=LaneBoundary(:broken, :white)
    )

    segments = RoadSegment[]
    seg1 = gen_straight_segment(1, nlanes, roadlength,
                                                origin=origin, lane_widths=lane_widths,
                                                boundary_leftmost=boundary_leftmost,
                                                boundary_rightmost=boundary_rightmost,
                                                boundary_middle=boundary_middle)
    seg2 = gen_straight_segment(2, nlanes, roadlength,
                      origin=VecSE2(origin.x+roadlength+0.5*lane_width,origin.y-0.5*lane_width,-pi/2),
                                                lane_widths=lane_widths,
                                                boundary_leftmost=boundary_leftmost,
                                                boundary_rightmost=boundary_rightmost,
                                                boundary_middle=boundary_middle)
    seg3 = gen_straight_segment(3, nlanes, roadlength,
        origin=VecSE2(origin.x+roadlength+(2*nlanes-0.5)*lane_width,origin.y-roadlength-0.5*lane_width,pi/2),
                                                lane_widths=lane_widths,
                                                boundary_leftmost=boundary_leftmost,
                                                boundary_rightmost=boundary_rightmost,
                                                boundary_middle=boundary_middle)
    seg4 = gen_straight_segment(4, nlanes, roadlength,
        origin=VecSE2(origin.x+roadlength+2*nlanes*lane_width,origin.y,0.0),
                                                lane_widths=lane_widths,
                                                boundary_leftmost=boundary_leftmost,
                                                boundary_rightmost=boundary_rightmost,
                                                boundary_middle=boundary_middle)

    seg5 = gen_straight_segment(5, nlanes, roadlength,
        origin=VecSE2(origin.x+roadlength+2*nlanes*lane_width+roadlength,origin.y+(2*nlanes-1.0)lane_width,pi),
                                                lane_widths=lane_widths,
                                                boundary_leftmost=boundary_leftmost,
                                                boundary_rightmost=boundary_rightmost,
                                                boundary_middle=boundary_middle)

    seg6 = gen_straight_segment(6, nlanes, roadlength,
        origin=VecSE2(origin.x+roadlength+(2*nlanes-0.5)*lane_width,origin.y+(2*nlanes-0.5)*lane_width,pi/2),
                                                lane_widths=lane_widths,
                                                boundary_leftmost=boundary_leftmost,
                                                boundary_rightmost=boundary_rightmost,
                                                boundary_middle=boundary_middle)
    seg7 = gen_straight_segment(7, nlanes, roadlength,
    origin=VecSE2(origin.x+roadlength+0.5*lane_width,origin.y+(2*nlanes-0.5)*lane_width+roadlength,-pi/2),
                                                lane_widths=lane_widths,
                                                boundary_leftmost=boundary_leftmost,
                                                boundary_rightmost=boundary_rightmost,
                                                boundary_middle=boundary_middle)
    seg8 = gen_straight_segment(8, nlanes, roadlength,
        origin=VecSE2(origin.x+roadlength,origin.y+(2*nlanes-1.0)*lane_width,pi),
                                                lane_widths=lane_widths,
                                                boundary_leftmost=boundary_leftmost,
                                                boundary_rightmost=boundary_rightmost,
                                                boundary_middle=boundary_middle)

    retval = Roadway()
    push!(retval.segments, seg1)
    push!(retval.segments, seg2)
    push!(retval.segments, seg3)
    push!(retval.segments, seg4)
    push!(retval.segments, seg5)
    push!(retval.segments, seg6)
    push!(retval.segments, seg7)
    push!(retval.segments, seg8)

    add_junction!(junction,retval)
    return retval,junction
end

function gen_loop_roadway(nlanes::Int;
    length::Float64=100.0,
    width::Float64=10.0,
    radius::Float64=25.0,
    ncurvepts_per_turn::Int=25, # includes start and end
    lane_width::Float64=DEFAULT_LANE_WIDTH, # [m]
    lane_widths::Vector{Float64} = fill(lane_width, nlanes),
    boundary_leftmost::LaneBoundary=LaneBoundary(:solid, :white),
    boundary_rightmost::LaneBoundary=LaneBoundary(:solid, :white),
    boundary_middle::LaneBoundary=LaneBoundary(:broken, :white),
    )

    ncurvepts_per_turn ≥ 2 || error("must have at least 2 pts per turn")

    A = VecSE2(0.0,0.0,0.0)
    B = VecSE2(length+radius, radius,pi/2)
    C = VecSE2(length, width + 2*radius,-pi)
    D = VecSE2(-radius, width + radius, -pi/2)

    seg1 = gen_straight_segment(1, nlanes, length,
                                                origin=A, lane_widths=lane_widths,
                                                boundary_leftmost=boundary_leftmost,
                                                boundary_rightmost=boundary_rightmost,
                                                boundary_middle=boundary_middle)
    seg2 = gen_straight_segment(2, nlanes, width,
                                                origin=B, lane_widths=lane_widths,
                                                boundary_leftmost=boundary_leftmost,
                                                boundary_rightmost=boundary_rightmost,
                                                boundary_middle=boundary_middle)
    seg3 = gen_straight_segment(3, nlanes, length,
                                                origin=C, lane_widths=lane_widths,
                                                boundary_leftmost=boundary_leftmost,
                                                boundary_rightmost=boundary_rightmost,
                                                boundary_middle=boundary_middle)
    seg4 = gen_straight_segment(4, nlanes, width,
                                                origin=D, lane_widths=lane_widths,
                                                boundary_leftmost=boundary_leftmost,
                                                boundary_rightmost=boundary_rightmost,
                                                boundary_middle=boundary_middle)
    retval = Roadway()
    push!(retval.segments, seg1)
    push!(retval.segments, seg2)
    push!(retval.segments, seg3)
    push!(retval.segments, seg4)

    for i=1:3
         connect_two_seg!(retval.segments[i],retval.segments[i+1],retval)
    end
    connect_two_seg!(retval.segments[4],retval.segments[1],retval)
    retval
end

function add_line!(origin::VecSE2,nlanes::Int,laneLength::Float64,roadway::Roadway;
    lane_width::Float64=DEFAULT_LANE_WIDTH, # [m]
    lane_widths::Vector{Float64} = fill(lane_width, nlanes),
    boundary_leftmost::LaneBoundary=LaneBoundary(:solid, :white),
    boundary_rightmost::LaneBoundary=LaneBoundary(:solid, :white),
    boundary_middle::LaneBoundary=LaneBoundary(:broken, :white)
    )
    id = length(roadway.segments) + 1
    seg = gen_straight_segment(id, nlanes, laneLength,
                                                origin=origin, lane_widths=lane_widths,
                                                boundary_leftmost=boundary_leftmost,
                                                boundary_rightmost=boundary_rightmost,
                                                boundary_middle=boundary_middle)
    push!(roadway.segments,seg)
end

max_stars_repo_name: sisl/UrbanDrivingSimulation

max_issues_repo_name: sisl/UrbanDrivingSimulation

max_forks_repo_name: sisl/UrbanDrivingSimulation


text:
@ghdef mutable struct CheckRun
    id::Union{Int, Nothing}
    head_sha::Union{String, Nothing}
    external_id::Union{String, Nothing}
    status::Union{String, Nothing}
    conclusion::Union{String, Nothing}
    started_at::Union{DateTime, Nothing}
    completed_at::Union{DateTime, Nothing}
    app::Union{App, Nothing}
    pull_requests::Union{Vector{PullRequest}, Nothing}
end
namefield(cr::CheckRun) = cr.id

@api_default function create_check_run(api::GitHubAPI, repo::Repo; headers = Dict(), kwargs...)
    headers["Accept"] = "application/vnd.github.antiope-preview+json"
    result = gh_post_json(api, "/repos/$(name(repo))/check-runs"; headers=headers, kwargs...)
    return CheckRun(result)
end

@api_default function update_check_run(api::GitHubAPI, repo::Repo, id::Int; headers = Dict(), kwargs...)
    headers["Accept"] = "application/vnd.github.antiope-preview+json"
    result = gh_patch_json(api, "/repos/$(name(repo))/check-runs/$(id)"; headers=headers, kwargs...)
    return CheckRun(result)
end


text:
"""
    parse_config(config_file)

Generate a YAML.jl config object by first passing `config_file` to Polyorder to do any configuration checks for its consistency. Here ·config_file` is path string.
"""
function parse_config(config_file)
    config = Polyorder.Config(config_file)
    data = StdString()
    # @show Polyorder.version(config)
    Polyorder.get_config_data(config, CxxRef(data))
    # @show data
    config = YAML.load(String(data))

    return config
end

function init_model(config_file)
    config = parse_config(config_file)

    modeltype = config["Model"]["model"]
    if uppercase(modeltype) == "AB" || uppercase(modeltype) == "AB"
        model = Polyorder.Model_AB(config_file)
    elseif uppercase(modeltype) == "AB_A" || uppercase(modeltype) == "AB+A"
        model = Polyorder.Model_AB_A(config_file)
    elseif uppercase(modeltype) == "AB3_A" || uppercase(modeltype) == "AB3+A"
        model = Polyorder.Model_AB3_A(config_file)
    else
        error("Model $modeltype is not supported by Polyorder!")
    end

    return model, config
end

function reset_model!(model, config)
    config_str = StdString(YAML.write(config))
    Polyorder.reset(model, CxxRef(config_str))
end

function run!(model, config)
    reset_model!(model, config)
    max_iter = config["Algorithm_SCFT"]["max_iter"]
    display_interval = config["IO"]["display_interval"]
    if display_interval < 1
        display_interval = max_iter + 1
    end
    thresh_residual = config["Algorithm_SCFT"]["thresh_residual"]
    for i in 1:max_iter
        Polyorder.update(model)
        Hi = Polyorder.Hi(model)
        H = Polyorder.H(model)
        err = Polyorder.residual_error(model)
        if i % display_interval == 0
            @info "SCFT iteration: $i"
            Polyorder.display(model)
        end
        if err < thresh_residual
            return H + Hi
        end
    end
    return Inf
end

function scft(x, model, config)
    dim = config["Grid"]["dimension"]
    a = config["UnitCell"]["a"]
    config["UnitCell"]["a"] = x
    if dim == 2
        b = config["UnitCell"]["b"]
        config["UnitCell"]["b"] = (x / a) * b
    end
    if dim == 3
        c = config["UnitCell"]["c"]
        config["UnitCell"]["c"] = (x / a) * c
    end

    @info "SCFT optimize cell: a = $a"

    F = run!(model, config)

    @assert F !== Inf "SCFT cannot achieve accuray within max_iter steps!"

    return F
end

"""
brent_init(guess, interval)

Return an interval `(a, b)` which feeds to Optim.solve with Optim.Brent() and the initial value for Brent method will be `guess` and `interval = b - a`.
"""
function brent_init(guess, interval)
    r = (2 - MathConstants.golden)
    a = guess - r * interval
    b = a + interval
    if a < 1.0
        a = 1.0
        b = (guess - a) / r + a
    end

    return a, b
end

function optimize_cell!(model, config)
    guess = config["Algorithm_Cell_Optimization"]["cell_guess"][1]
    interval = config["Algorithm_Cell_Optimization"]["cell_interval"][1]
    abs_tol = config["Algorithm_Cell_Optimization"]["tol_cell"]
    store_trace = config["Algorithm_Cell_Optimization"]["store_trace"]
    show_trace = config["Algorithm_Cell_Optimization"]["show_trace"]
    lower, upper = brent_init(guess, interval)
    result = Optim.optimize(x -> scft(x, model, config), lower, upper, Optim.Brent(); abs_tol=abs_tol, store_trace=store_trace, show_trace=show_trace)

    return result
end


text:
struct Wolff{TF<:AbstractFloat, TI<:Integer} <:AbstractSampler
    β::TF
    J::TF
    L::TI
    p::TF
end

"""
    `Wolff(β, J, L)` - Set up the Wolff sampler

### Fields
* `β` - Inverse temperature
* `J` - Coupling constant>0
* `L` - Lattice size, L×L
"""
function Wolff(β::TF, J::TF, L::TI) where{TF<:AbstractFloat, TI<:Integer}
    return Wolff{TF, TI}(β, J, L, 1 - exp(-2 * β * J))
end

mutable struct WolffState{TI<:Integer, Tx<:AbstractArray{TI}} <:AbstractSamplerState
    x::Tx
    D::TI
end

function InitState!(x, sampler::Wolff)

    return WolffState(x, zero(eltype(x)));
end

function InitState(x₀, sampler::Wolff)

    return WolffState(deepcopy(x₀),zero(eltype(x₀)));
end

function UpdateState!(state::WolffState, sampler::Wolff)

    i = rand(1:sampler.L);
    j = rand(1:sampler.L);
    state.D = state.x[i,j];
    state.x[i,j] *=-1;

    k,l = mod1(i+1,sampler.L), j;
    if (state.x[k,l]==state.D)
        if (rand() < sampler.p)
            state.x[k,l] *=-1;
            FlipNeighbors!(state.x, k, l, sampler.p, state.D, sampler.L)
        end
    end

    k,l = mod1(i-1,sampler.L), j;
    if (state.x[k,l]==state.D)
        if (rand() < sampler.p)
            state.x[k,l] *=-1;
            FlipNeighbors!(state.x, k, l, sampler.p, state.D, sampler.L)
        end
    end

    k,l = i, mod1(j+1,sampler.L);
    if (state.x[k,l]==state.D)
        if (rand() < sampler.p)
            state.x[k,l] *=-1;
            FlipNeighbors!(state.x, k, l, sampler.p, state.D, sampler.L)
        end
    end

    k,l = i, mod1(j-1,sampler.L);
    if (state.x[k,l]==state.D)
        if (rand() < sampler.p)
            state.x[k,l] *=-1;
            FlipNeighbors!(state.x, k, l, sampler.p, state.D, sampler.L)
        end
    end
    state
end

"""
`FlipNeighbors!` - Recursively flip neighbors in the Wolff algorithm for 2D
square lattice Ising

### Fields
* `x` - state of the lattice
* `i` - i-coordinate of current spin
* `j` - j-coordinate of current spin
* `p` - Flipping probability
* `D` - Flipping direction
* `L` - Lattice size, N×N
"""
function FlipNeighbors!(x, i, j, p, D, L)

    k,l = mod1(i+1,L), j;
    if (x[k,l]==D)
        if (rand() < p)
            x[k,l] *=-1;
            FlipNeighbors!(x, k, l, p, D, L)
        end
    end

    k,l = mod1(i-1,L), j;
    if (x[k,l]==D)
        if (rand() < p)
            x[k,l] *=-1;
            FlipNeighbors!(x, k, l, p, D, L)
        end
    end

    k,l = i, mod1(j+1,L);
    if (x[k,l]==D)
        if (rand() < p)
            x[k,l] *=-1;
            FlipNeighbors!(x, k, l, p, D, L)
        end
    end

    k,l = i, mod1(j-1,L);
    if (x[k,l]==D)
        if (rand() < p)
            x[k,l] *=-1;
            FlipNeighbors!(x, k, l, p, D, L)
        end
    end
    x
end


text:
#=
P08 (**) Eliminate consecutive duplicates of list elements.
If a list contains repeated elements they should be replaced with a single copy of the element. The order of the elements should not be changed.

Example:
* (compress '(a a a a b c c a a d e e e e))
(A B C A D E)
=#

function _compress(coll, compressed)
    if isempty(coll)
        compressed
    else
        compressed = compressed[end] == coll[1] ? compressed : vcat(compressed, coll[1])
        _compress(coll[2:end], compressed)
    end
end

function compress(coll)
    _compress(coll[2:end], [coll[1]])
end

@assert compress([:a, :a, :a, :a, :b, :c, :c, :a, :a, :d, :e, :e, :e, :e]) == [:a, :b, :c, :a, :d, :e]

println("Tests passed: JL-08.jl")


text:
include("organize_binvisits.jl")

import ..TransferOperatorRectangularBinning
"""
    estimate_transferoperator_from_binvisits(bv::BinVisits;
                allocate_frac = 1,
                boundary_condition = :none)

Estimate transfer operator from information about which bin
gets visited by which points of the orbit.

## Transition probabilities
To determine the transition probabilities, we pre-allocate two index vectors
(`I` and `J`) and a probability vector (`P`). After the transition probabilities
have been determined, the vectors are combined into a sparse matrix (the
transfer matrix), which typically have very few nonzero entries (often only a
few percent).

## Memory allocation
To not allocate too much memory, `allocate_frac` controls what fraction of the
total number of possible transitions (``n_{states}^2``) we pre-allocate for.

## Example usage

```
# Random orbit.
orbit = rand(1000, 3)
x, y, z = orbit[:, 1], orbit[:, 2], orbit[:, 3]

# Embedding E = {(y(t+1), y(t), x(t))}
E = cembed([x, y], [2, 2, 1], [1, 0, 0]);

# Which bin sizes to use along each dimension?
ϵ = [.4, .2, .4]

# Identify which bins of the partition resulting
# from using ϵ each point of the embedding visits.
visited_bins = assign_bin_labels(E, ϵ)

# Which are the visited bins, which points
# visits which bin, repetitions, etc...
binvisits = get_binvisits(visited_bins)

# Use that information to estimate transfer operator
TO = transferoperator(binvisits)

# Verify that TO is Markov (NB. last row might be zero!,
# then this test might fail)
all(sum(TO, 2) .≈ 1)
```
"""
function estimate_transferoperator_from_binvisits(bv::BinVisits;
                allocate_frac = 1.0,
                boundary_condition = :none)

    valid_boundary_conditions = [:none, :exclude, :circular, :invariantize]
    if !(boundary_condition ∈ valid_boundary_conditions)
        error("Boundary condition $boundary_condition not valid.")
    end

    if !(0 < allocate_frac <= 1)
        error("allocate_frac needs to be on the interval (0, 1]")
    end

    first_visited_by = bv.first_visited_by
    visitors = bv.visitors
    visits_whichbin = bv.visits_whichbin
    n_visited_bins = length(first_visited_by)

    # Initialise transfer (Perron-Frobenius) operator as a sparse matrix
    # (keep track of indices and values in separate columns for now)
    n_possible_nonzero_entries = n_visited_bins^2
    N = ceil(Int, n_possible_nonzero_entries / allocate_frac)

    I = zeros(Int32, N)
    J = zeros(Int32, N)
    P = zeros(Float64, N)

    # Preallocate target index for the case where there is only
    # one point of the orbit visiting a bin.
    target_bin_j::Int = 0
    n_visitsᵢ::Int = 0

    # Keep track of how many transitions we've considered.
    transition_counter = 0

    if boundary_condition == :circular
        #warn("Using circular boundary condition")
        append!(visits_whichbin, [1])
    elseif boundary_condition == :random
        #warn("Using random circular boundary condition")
		append!(visits_whichbin, [rand(1:length(visits_whichbin))])
	end

    # Loop over the visited bins bᵢ
    for i in 1:n_visited_bins
        # How many times is this bin visited?
        n_visitsᵢ = length(visitors[i])

        # If both conditions below are true, then there is just one
        # point visiting the i-th bin. If there is only one visiting point and
        # it happens to be the last, we skip it, because we don't know its
        # image.
        if n_visitsᵢ == 1 && !(i == visits_whichbin[end])
            transition_counter += 1
            # To which bin does the single point visiting bᵢ jump if we
            # shift it one time step ahead along its orbit?
            target_bin_j = visits_whichbin[visitors[i][1] + 1][1]

            # We now know that exactly one point (the i-th) does the
            # transition from i to the target j.
            I[transition_counter] = i
            J[transition_counter] = target_bin_j
            P[transition_counter] = 1
        end
        # If more than one point of the orbit visits the i-th bin, we
        # identify the visiting points and track which bins bⱼ they end up
        # in after the forward linear map of the points.
        if n_visitsᵢ > 1
            timeindices_visiting_pts = visitors[i]

            # TODO: Introduce circular boundary condition. Simply excluding
            # might lead to a cascade of loosing points.

            # If bᵢ is the bin visited by the last point in the orbit, then
            # the last entry of `visiting_pts` will be the time index of the
            # last point of the orbit. In the next time step, that point will
            # have nowhere to go along its orbit (precisely because it is the
            # last data point). Thus, we exclude it.
            if i == visits_whichbin[end]
                #warn("Removing last point")
                n_visitsᵢ = length(timeindices_visiting_pts) - 1
                timeindices_visiting_pts = timeindices_visiting_pts[1:(end - 1)]
            end

            # To which boxes do each of the visitors to bᵢ jump in the next
            # time step?
            target_bins = visits_whichbin[timeindices_visiting_pts .+ 1]
            unique_target_bins = unique(target_bins)

            # Count how many points jump from the i-th bin to each of
            # the unique target bins, and use that to calculate the transition
            # probability from bᵢ to bⱼ.
            for j in 1:length(unique_target_bins)
                n_transitions_i_to_j = sum(target_bins .== unique_target_bins[j])

                transition_counter += 1
                I[transition_counter] = i
                J[transition_counter] = unique_target_bins[j]
                P[transition_counter] = n_transitions_i_to_j / n_visitsᵢ
            end
        end
    end

    # Combine indices and values into a sparse representation of the Perron-
    # Frobenius operator (transfer operator). Filter out the nonzero entries.
    # The number of rows in the matrix is given by the number of unique points
    # we have in the embedding.
    TO = Array(sparse(I[1:transition_counter],
                J[1:transition_counter],
                P[1:transition_counter],
                n_visited_bins, n_visited_bins))

    # There may be boxes which are visited by points of the orbit, but not by
    # any image points.
    # zc = zerocols(TO)
    # l = length(zc)
    # if l > 0
    #     col = zc[1]
    #     warn("There were $l all-zero columns. Column $col is the culprit -> normalizing ...")
    #     TO = TO ./ sum(TO, 2)
    # end
    # If the first bin is not visited, make every point starting in that bin
    # jump to the same bin with probability one.
    # zc = zerocols(TO)
    # if length(zc) > 0
    #     warn("There were $l all-zero columns: column $col")
    #     @show zc
    #     i = zc[1]
    #     TO[i, 1] = 1.0
    # end

    TransferOperatorRectangularBinning(TO)
end


text:
using Images, ColorSchemes, ImageFiltering, ImageTransformations
using Pipe
using CUDA

# Function for counting the steps to divergence.
function recurse_steps_to_satisfy(
    f::Function, condition::Function, max_steps::Integer)
    z = Complex(0.0, 0.0)
    for i = 1:max_steps
        z = f(z)
        if condition(z)
            return i
        end
    end
    return max_steps + 1
end

function mandelbrot_function(c::Complex)
    return (z) -> z * z + c
end

function mandelbrot_condition(z::Complex)
    return abs2(z) >= 4
end

function mandelbrot_steps(c::Complex, max_steps::Integer)
    return recurse_steps_to_satisfy(
        mandelbrot_function(c), mandelbrot_condition, max_steps)
end

function get_cmap(colorscheme::ColorScheme, max_steps::Integer)
    return [colorscheme[num_steps/(max_steps+1)]
            for num_steps = 1:max_steps+1]
end

function get_complex_domain(pixel_size, center_loc, width, height, resample_factor)
    horizontal_extent = width * pixel_size
    vertical_extent = height * pixel_size
    real_domain = (
        center_loc[1]-horizontal_extent/2:
        pixel_size/resample_factor:
        center_loc[1]+horizontal_extent/2
        )
    imaginary_domain = (
        center_loc[2]-vertical_extent/2:
        pixel_size/resample_factor:
        center_loc[2]+vertical_extent/2
        )

    # The imaginary domain is reversed so that the image
    # comes out with the correct orientation.
    imaginary_domain = reverse(imaginary_domain)

    complex_domain = [Complex(x, y) for y in imaginary_domain, x in real_domain]

    return complex_domain
end

# Color space.
colorscheme = ColorSchemes.inferno
max_steps = 2500

color_map = get_cmap(colorscheme, max_steps)

# Lay out the complex domain to evaluate.
pixel_size = 0.0000001
width = 1920
height = 1080
center_loc = (-0.75393, 0.05)
resample_factor = 2

complex_domain = get_complex_domain(pixel_size, center_loc, width, height, resample_factor)

# Antialiasing filter.
blur_kernel = Kernel.gaussian(resample_factor / 2)

@time begin
    image = @pipe complex_domain |>
        CuArray |> # Casting as CuArray sends the computation to the GPU.
        mandelbrot_steps.(_, max_steps) |>
        Array |>
        color_map[_] |>
        imfilter(_, blur_kernel) |>
        imresize(_, (height, width))
end

save("image.png", image)

max_stars_repo_name: nicolassmith/fractal-art

max_issues_repo_name: nicolassmith/fractal-art

max_forks_repo_name: nicolassmith/fractal-art


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct OperationProperties <: SwaggerModel
    serviceSpecification::Any # spec type: Union{ Nothing, ServiceSpecification } # spec name: serviceSpecification

    function OperationProperties(;serviceSpecification=nothing)
        o = new()
        validate_property(OperationProperties, Symbol("serviceSpecification"), serviceSpecification)
        setfield!(o, Symbol("serviceSpecification"), serviceSpecification)
        o
    end
end # type OperationProperties

const _property_map_OperationProperties = Dict{Symbol,Symbol}(Symbol("serviceSpecification")=>Symbol("serviceSpecification"))
const _property_types_OperationProperties = Dict{Symbol,String}(Symbol("serviceSpecification")=>"ServiceSpecification")
Base.propertynames(::Type{ OperationProperties }) = collect(keys(_property_map_OperationProperties))
Swagger.property_type(::Type{ OperationProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_OperationProperties[name]))}
Swagger.field_name(::Type{ OperationProperties }, property_name::Symbol) =  _property_map_OperationProperties[property_name]

function check_required(o::OperationProperties)
    true
end

function validate_property(::Type{ OperationProperties }, name::Symbol, val)
end


text:
module Firebase

# Includes
include("firestore.jl")
include("firestore-token.jl")
include("firebaseauth.jl")
include("convert.jl")
include("cloudstorage.jl")
include("realtime.jl")

# Imports
using HTTP
using GoogleCloud
using JSON
using Glob

# Utils
export get
export get_collection, get_singlepage,get_document, get_request
export convert_to_dict
export convert_field
export convert_timestamp
export convert_map
export convert_array

# Firebase Authentication related
export set_webapikey
export firebase_signup
export firebase_signin
export firebase_signinanon
export firebase_signinoauth
export firebase_fetchproviders
export firebase_passwordreset
export verify_passwordresetcode
export confirm_passwordreset
export firebase_changeemail
export firebase_changepassword
export firebase_updateprofile
export firebase_getuserinfo
export firebase_linkuserinfo
export firebase_unlinkprovide
export firebase_verifyemail
export firebase_confirmemailverify
export firebase_deleteuser

# Cloud Firestore
export init
export projectid
export get_token
export auth_header
export firestore_batchget
export firestore_batchwrite
export firestore_beginTransaction
export firestore_createdoc
export firestore_commit
export firestore_deletedoc
export firestore_getcollectionids
export firestore_getdoclist
export firestore_partitionquery
export firestore_patch
export firestore_rollback
export firestore_runquery

# Realtime Database
export realdb_init
export realdb_get
export realdb_post
export realdb_delete
export realdb_put
export realdb_patch

# Cloud Storage
export cloudstore_init
export cloudstore_get
export cloudstore_sendfile

end # module


text:
using DataFrames, NamedArrays
import DocOpt
using JLD

doc = """CSV Process Utility

Usage:
  process_csv.jl check <csvfile> [--exclude=<excludefile>]
  process_csv.jl statuscounts <csvfile> <csvfileout> [--exclude=<excludefile>] [--bestof=<list>]...
  process_csv.jl geomeans <csvfile> [--exclude=<excludefile>] [--bestof=<list>]...
  process_csv.jl perfprofile <csvfile> <output_jld> <solver>... [--exclude=<excludefile>] [--bestof=<list>]...
  process_csv.jl -h | --help

Options:
  -h --help               Show this screen.
  --exclude=<excludefile> A list of runs to exclude.
  --bestof=<list>         Create an artificial solver by taking the best of multiple solvers. Format: "A B C" creates a new solver A from best of B and C.
"""

arguments = DocOpt.docopt(doc)

# process the csv produced by process_awsoutput.jl
results = readtable(arguments["<csvfile>"])

# file listing manually excluded results
# format by line: solver instance ignored comment
# example: PAJ_CPLEX_MOSEK clay0304h.cbf # CPLEX gave wrong answer
if isa(arguments["--exclude"], String)
    excluded = [split(strip(l))[1:2] for l in split(readstring(arguments["--exclude"]),'\n') if strip(l) != ""]
else
    excluded = []
end
function isexcluded(solver,inst)
    for k in 1:length(excluded)
        if excluded[k][1] == solver && excluded[k][2] == inst
            return true
        end
    end
    return false
end

for i in 1:size(results,1)
    if isexcluded(results[i,:solver],results[i,:instance])
        results[i,:status] = "Excluded"
    end
    if isna(results[i,:status])
        results[i,:status] = "Blank"
    end
end

for solverlist in arguments["--bestof"]
    newsolver = split(solverlist)[1]
    best_of_solvers = split(solverlist)[2:end]
    newrows = []
    for rowlist in groupby(results,:instance)
        besttime = Inf
        beststatus = "Blank"
        inst = rowlist[1,:instance]
        for i in 1:size(rowlist,1)
            if rowlist[i,:solver] in best_of_solvers
                if rowlist[i,:status] == "Optimal"
                    besttime = min(besttime,rowlist[i,:solvertime])
                    beststatus = "Optimal"
                elseif rowlist[i,:status] in ["UserLimit","KilledTime","KilledMemory"]
                    if beststatus != "Optimal"
                        besttime = min(besttime,rowlist[i,:solvertime])
                        beststatus = "UserLimit"
                    end
                end
            end
        end
        newrow = Dict{Symbol,Any}(n => NA for n in names(results))
        newrow[:instance] = inst
        newrow[:solver] = newsolver
        newrow[:status] = beststatus
        newrow[:solvertime] = besttime
        push!(newrows,newrow)
    end
    for r in newrows
        push!(results, r)
    end
end

optimal_runs = results[results[:status] .== "Optimal",:]

function printdf(df)
    # lifted from DataFrames.show()
    rowindices1 = 1:size(df,1)
    rowindices2 = 1:0
    maxwidths = DataFrames.getmaxwidths(df, rowindices1, rowindices2, :Rows)
    DataFrames.showrows(STDOUT, df, 1:size(df,1), 1:0, maxwidths, false, :Row, true)
    println()
end

if arguments["check"]

    # problematic violations
    viol = optimal_runs[(optimal_runs[:max_linear_violation] .> 1e-3) | (optimal_runs[:max_soc_violation] .> 1e-3) | (optimal_runs[:max_socrot_violation] .> 1e-3) | (optimal_runs[:max_exp_violation] .> 1e-3) | (optimal_runs[:max_int_violation] .> 1e-3), :]

    if size(viol,1) == 0
        println("No problematic violations")
    else
        println("Problematic violations")
        println(viol)
    end

    # check for disagreements in objective value
    for optval_by_instance in groupby(optimal_runs, :instance)
        first_optval = optval_by_instance[1,:objval_reported]
        sense = optval_by_instance[1,:sense]
        inst = optval_by_instance[1,:instance]
        for i in 2:size(optval_by_instance,1)
            solver = optval_by_instance[i,:solver]
            optval = optval_by_instance[i,:objval_reported]
            if abs(optval-first_optval)/abs(first_optval+1e-5) > 1e-4
                println("Objective disagreement on instance $inst (sense = $sense)")
                printdf(optval_by_instance[:,[:solver,:objval_reported,:objbound]])
                break
            end
        end
    end

    # check for duplicated runs
    for g1 in groupby(results, :instance)
        for g2 in groupby(g1,:solver)
            if size(g2,1) > 1
                s = g2[1,:solver]
                inst = g2[1,:instance]
                println("Multiple results for solver $s instance $inst")
            end
        end
    end

elseif arguments["statuscounts"]

    # status counts by solver

    statuses = sort(collect(unique(results[:status])))
    all_solvers = sort(collect(unique(results[:solver])))
    status_table = NamedArray(zeros(Int,length(all_solvers),length(statuses)+1))
    setnames!(status_table, all_solvers, 1)
    setnames!(status_table, [statuses;"Total"], 2)
    for i in 1:size(results,1)
        status_table[results[i,:solver],results[i,:status]] += 1
        status_table[results[i,:solver],"Total"] += 1
    end

    df2 = DataFrame(convert(Array,status_table))
    names!(df2, convert(Vector{Symbol},names(status_table,2)))
    df2[:solver] = names(status_table,1)
    writetable(arguments["<csvfileout>"],df2)

elseif arguments["geomeans"]

    function shifted_geomean(table, field, shift, groupby)
        solvers = unique(table[groupby])
        sum_by = Dict(g => 0.0 for g in solvers)
        counter_by = Dict(g => 0 for g in solvers)
        # prod(x_i + s)^(1/n) - s
        # compute prod(x_i+s)^(1/n) as
        # exp(sum(log(x_i+s))/n)
        for i in 1:size(table,1)
            if isna(table[i,field])
                continue
            end
            sum_by[table[i,groupby]] += log(table[i,field] + shift)
            counter_by[table[i,groupby]] += 1
        end
        for g in solvers
            r = exp(sum_by[g]/counter_by[g]) - shift
            println("$g: $r")
        end
    end

    time_shift = 10.0
    subproblem_shift = 1.0

    # tidy data (needed for KilledX cases)
    for i in 1:size(results,1)
        if isna(results[i,:solvertime])
            results[i,:solvertime] = results[i,:timelimit]
        end
    end

    println("Shifted geomean of solve time on all instances")
    shifted_geomean(results, :solvertime, time_shift, :solver)
    println()

    println("Shifted geomean of solve time on optimal instances")
    shifted_geomean(optimal_runs, :solvertime, time_shift, :solver)
    println()

    println("Shifted geomean of conic subproblem count on all instances")
    shifted_geomean(results, :conic_subproblem_count, subproblem_shift, :solver)
    println()

    println("Shifted geomean of conic subproblem count on optimal instances")
    shifted_geomean(optimal_runs, :conic_subproblem_count, subproblem_shift, :solver)

elseif arguments["perfprofile"]
    solvers = arguments["<solver>"]
    # need a table of [time] X [solver] where each row is a solver
    time_rows = []
    itercount_rows = []
    instance_names = String[]
    for by_instance in groupby(optimal_runs, :instance)
        time_row = fill(Inf,length(solvers))
        itercount_row = fill(Inf,length(solvers))
        for i in 1:size(by_instance,1)
            push!(instance_names,by_instance[1,:instance])
            if by_instance[i,:status] == "Optimal"
                whichsolver = indexin([by_instance[i,:solver]],solvers)[1]
                if whichsolver != 0
                    time_row[whichsolver] = by_instance[i,:solvertime]
                    if !isna(by_instance[i,:conic_subproblem_count])
                        itercount_row[whichsolver] = by_instance[i,:conic_subproblem_count]
                    end
                end # else not plotting this solver
            end # else not solved, Inf by default
        end
        push!(time_rows,time_row')
        push!(itercount_rows,itercount_row')
    end
    time_mat = vcat(time_rows...)
    itercount_mat = vcat(itercount_rows...)
    @assert endswith(arguments["<output_jld>"],".jld")
    save(arguments["<output_jld>"], "time_table", time_mat, "itercount_table", itercount_mat, "solvers", solvers, "instance_names", instance_names)
    #p = performance_profile(mat, solvers, title="TITLE")
    #Plots.savefig("perf.png")
end

max_stars_repo_name: mlubin/PajaritoSupplement

max_issues_repo_name: mlubin/PajaritoSupplement

max_forks_repo_name: mlubin/PajaritoSupplement