text:
#=
Many of the base level routines here are adapted from the SC Toolbox by Toby Driscoll
and are inspired from the work of L. Trefethen, "Numerical Computation of the
Schwarz-Christoffel Transformation", STAN-CS-79-710, 1979.
=#

#= functions for the Schwarz-Christoffel exterior map from unit disk =#
function param(w::Vector{ComplexF64},β::Vector{Float64},
                 ζ0::Vector{ComplexF64},
                 qdat::Tuple{Array{Float64,2},Array{Float64,2}})
#=
Solve for the parameters of the exterior Schwarz-Christoffel mapping: the
prevertices `ζ` and the constant factor `c`. The routine requires the
vertices `w` in clockwise order, the exterior turning angles `β`,
a list of guesses for the prevertices `ζ0` (which can be blank), and
the Gauss-Jacobi quadrature data in `qdat`
=#

  n = length(w)
  if n == 2
    ζ = ComplexF64[-1,1]
  else
    len = abs.(diff(circshift(w,1)))
    nmlen = abs.(len[3:n-1]/len[2])
    if isempty(ζ0)
      y0 = zeros(n-1)
    else
      ζ0 = ζ0/ζ0[n]
      θ = angle.(ζ0)
      θ[θ.<=0] = θ[θ.<=0] .+ 2π
      dt = diff([0;θ[1:n-1];2π])
      @. y0 = log(dt[1:n-1]/dt[2:n])
    end

    depfun! = Depfun(β,nmlen,qdat)

    F0 = similar(y0)
    df = OnceDifferentiable(depfun!, y0, F0)
    sol = nlsolve(df,y0,show_trace = :false)

    ζ = zeros(ComplexF64,n)
    θ = zeros(n-1)
    y_to_ζ!(ζ,θ,sol.zero)

  end

  mid = ζ[1]*exp(0.5*im*angle(ζ[2]/ζ[1]))
  dequad = DQuad(β,qdat)
  c = (w[2]-w[1])/(dequad([ζ[1]],[mid],[1],ζ)[1]-dequad([ζ[2]],[mid],[2],ζ)[1])

  return ζ, c

end

function evaluate_exterior(ζ::Vector{ComplexF64},w::Vector{ComplexF64},
                  β::Vector{Float64},prev::Vector{ComplexF64},
                  c::ComplexF64,qdat::Tuple{Array{Float64,2},Array{Float64,2}})
  #=
  Evaluates the exterior Schwarz-Christoffel mapping at `ζ`, which is
  presumed to be inside the unit circle. The vector `w` are the vertices (in
  clockwise order), `β` are the exterior turning angles (also in clockwise order),
  `prev` are the prevertices on the unit circle, and `c` is the constant factor
  of the transformation. The tuple `qdat` contains the Gauss-Jacobi nodes and
  weights.
  =#

  if isempty(ζ)
    nothing
  end

  n = length(w)
  neval = length(ζ)
  tol = 10.0^(-size(qdat[1],1))

  # set up the integrator
  dequad = DQuad(β,qdat)

  # initialize the mapped evaluation points
  z = zeros(ComplexF64,neval)

  # find the closest prevertices to each evaluation point and their
  #  corresponding distances
  dz = abs.(repeat(ζ, 1, n) - repeat(transpose(prev), neval))
  #dz = abs.(hcat([ζ for i=1:n]...)-vcat([transpose(prev) for i=1:neval]...))
  (dist,ind) = findmin(dz, dims = 2)
  #sing = floor.(Int,(ind[:]-1)/neval)+1
  sing = [I[2] for I in ind]

  # find any prevertices in the evaluation list and set them equal to
  #  the corresponding vertices. The origin is also a singular point
  vertex = (dist[:] .< tol)
  z[vertex] = w[sing[vertex]]
  zerop = abs.(ζ) .< tol
  z[zerop] .= Inf
  vertex = vertex .| zerop

  # the starting (closest) singularities for each evaluation point
  prevs = prev[sing]

  # set the initial values of the non-vertices
  z[.!vertex] = w[sing[.!vertex]]

  # distance to singularity at the origin
  absζ = abs.(ζ)

  # unfinished cases
  unf = .!vertex

  # set the integration starting points
  ζold = copy(prevs)
  ζnew = copy(ζold)
  dist = ones(neval)
  while any(unf)
    # distance to integrate still
    dist[unf] = min.(1,2*absζ[unf]./abs.(ζ[unf]-ζold[unf]))

    # new integration end point
    ζnew[unf] = ζold[unf] + dist[unf].*(ζ[unf]-ζold[unf])

    # integrate
    z[unf] = z[unf] + c*dequad(ζold[unf],ζnew[unf],sing[unf],prev)

    # set new starting integration points for those that can be integrated
    #  further
    unf = dist .< 1
    ζold[unf] = ζnew[unf]

    # only the first step can have a singularity
    sing .= 0

  end

  return z

end

function evalderiv_exterior(ζ::Vector{ComplexF64},β::Vector{Float64},
                  prev::Vector{ComplexF64},c::ComplexF64)

#=
Evaluates the first and second derivative of the exterior Schwarz-Christoffel
mapping at `ζ`, which is
presumed to be inside the unit circle. The vector `β` are the exterior turning
angles (also in clockwise order), `prev` are the prevertices on the unit circle, and
`c` is the constant factor of the transformation.
=#
    n = length(prev)
    neval = length(ζ)
    β = [β;-2]
    terms = [hcat([1 .- ζ/prev[i] for i = 1:n]...) ζ]
    dz = c*exp.(log.(terms)*β)
    terms2 = [[-β[i]/prev[i] for i = 1:n];β[n+1]]
    ddz = dz.*((1.0./terms)*terms2)
    return dz, ddz

end

function evalinv_exterior(z::Vector{ComplexF64},w::Vector{ComplexF64},
                  β::Vector{Float64},prev::Vector{ComplexF64},
                  c::ComplexF64,qdat::Tuple{Array{Float64,2},Array{Float64,2}})

#=
Evaluates the inverse of the exterior Schwarz-Christoffel mapping, using a combination
of integration and Newton iteration, using techniques from Trefethen (1979).
=#

   n = length(w)
   ζ = zeros(ComplexF64,size(z))
   lenz = length(z)
   ζ0 = []
   maxiter = 10
   tol = 10.0^(-size(qdat[1],1))

   # Find z values close to vertices and set ζ to the corresponding
   # prevertices
   done = zeros(Bool,size(z))
   for j = 1:n
     idx = findall(abs.(z .- w[j]) .< 3*eps())
     ζ[idx] .= prev[j]
     done[idx] .= true
   end
   lenz -= sum(done)
   if lenz==0
     return ζ
   end

   # Now, for remaining z values, first try to integrate
   #  dζ/dt = (z - z(ζ₀))/z'(ζ) from t = 0 to t = 1,
   # with the initial condition ζ(0) = ζ₀.
   if isempty(ζ0)
     ζ0,z0 = initial_guess(z,w,β,prev,c,qdat)
   else
     z0 = evaluate_exterior(ζ0,w,β,prev,c,qdat)
     if length(ζ0)==1 && lenz > 1
       ζ0 = repeat(transpose(ζ0), lenz)
       z0 = repeat(transpose(z0), lenz)
     end
     z0 = z0[.~done]
     ζ0 = ζ0[.~done]
   end
   odetol = max(tol,1e-3)
   scale = z[.~done] - z0

   ζ0 = [real(ζ0);imag(ζ0)]

   f(ζ,p,t) = invfunc(ζ,scale,β,prev,c)
   tspan = (0.0,1.0)
   prob = ODEProblem(f,ζ0,tspan)
   sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
   lenu = length(ζ0)
   ζ[.~done] = sol.u[end][1:lenz]+im*sol.u[end][lenz+1:lenu];
   out = abs.(ζ) .> 1
   ζ[out] = sign.(ζ[out])

   # Now use Newton iterations to improve the solution
   ζn = ζ
   k = 0
   while ~all(done) && k < maxiter
     F = z[.~done] - evaluate_exterior(ζn[.~done],w,β,prev,c,qdat)
     M = length(F)
     dF, ddz = evalderiv_exterior(ζn[.~done],β,prev,c)
     ζn[.~done] = ζn[.~done] + F./dF

     done[.~done] = abs.(F).< tol
     k += 1
   end
   F = z[.~done] - evaluate_exterior(ζn[.~done],w,β,prev,c,qdat)
   if any(abs.(F).> tol)
     error("Check solution")
   end
   ζ = ζn

end

function invfunc(u,scale,β::Vector{Float64},prev::Vector{ComplexF64},c::ComplexF64)
    lenu = length(u)
    lenzp = Int(lenu/2)
    ζ = u[1:lenzp]+im*u[lenzp+1:lenu]

    dz, ddz = evalderiv_exterior(ζ,β,prev,c)
    f = scale./dz
    zdot = [real(f);imag(f)]
end

function initial_guess(z::Vector{ComplexF64},w::Vector{ComplexF64},
                  β::Vector{Float64},prev::Vector{ComplexF64},
                  c::ComplexF64,qdat::Tuple{Array{Float64,2},Array{Float64,2}})

  n = length(w)
  tol = 1000.0*10.0^(-size(qdat[1])[1])
  shape = copy(z)
  ζ0 = copy(z)
  z0 = copy(z)
  argz = angle.(prev);
  argz[argz.<=0] .+= 2π

  argw = cumsum([angle(w[2]-w[1]); -π*β[2:n]])

  infty = isinf.(w)
  fwd = circshift(1:n,-1)
  anchor = zero(w)
  anchor[.~infty] = w[.~infty]
  anchor[infty] = w[fwd[infty]]
  direcn = exp.(im*argw)
  direcn[infty] = -direcn[infty]
  len = abs.(w[fwd] - w)

  factor = 0.5
  done = zeros(Bool,length(z))
  M = length(z)
  iter = Int(0)

  A = zeros(Float64,2,2)

  ζbase = NaN*ones(ComplexF64,n)
  zbase = NaN*ones(ComplexF64,n)
  idx = []
  while M > 0
    for j = 1:n
      if j<n
        ζbase[j] = exp(im*(factor*argz[j] + (1-factor)*argz[j+1]))
      else
        ζbase[j] = exp(im*(factor*argz[n] + (1-factor)*(2π+argz[1])))
      end
    end
    zbase = evaluate_exterior(ζbase,w,β,prev,c,qdat)
    proj = real.( (zbase-anchor) .* conj.(direcn) )
    zbase = anchor + proj.*direcn
    if isempty(idx)
      dist,idxtemp = findmin(abs.( repeat(transpose(z[.~done]), n) - repeat(zbase, 1, M)), dims = 1)
      for k = 1:M
          push!(idx,idxtemp[k][1])
      end
    else
      idx[.~done] = idx[.~done].%n .+ 1
    end
    ζ0[.~done] = ζbase[idx[.~done]]
    z0[.~done] = zbase[idx[.~done]]

    for j = 1:n
      active = (idx.==j) .& (.~done)
      if any(active)

        done[active] = ones(Bool,sum(active))
        for k in [1:j-1;j+1:n]'
          A[:,1] = [real(direcn[k]);imag(direcn[k])]
          for p in findall(active)
            dif = z0[p]-z[p]
              A[:,2] = [real(dif);imag(dif)]
              if cond(A) < eps()
                zx = real( (z[p]-anchor[k]) / direcn[k] )
                z0x = real( (z0[p]-anchor[k]) / direcn[k] )
                if (zx*z0x < 0) || ((zx-len[k])*(z0x-len[k]) < 0)
                  done[p] = false
                end
              else
                dif = z0[p]-anchor[k]
                  s = A\[real(dif);imag(dif)]
                  if s[1]>=0 && s[1]<=len[k]
                    if abs(s[2]-1) < tol
                      ζ0[p] = ζbase[k]
                      z0[p] = zbase[k]
                    elseif abs(s[2]) < tol
                      if real( conj(z[p]-z0[p])*im*direcn[k] ) > 0
                        done[p] = false
                      end
                    elseif s[2] > 0 && s[2] < 1
                      done[p] = false
                    end
                  end
                end # cond(A)
              end # for p
            end # for k
            M = sum(.~done)
            if M == 0
              break
            end
          end # if any active
          if iter > 2*n
            error("Can''t seem to choose starting points.  Supply them manually.")
          else
            iter += 1
          end
          factor = rand(1)[1]
        end # for j

  end  # while M

  return ζ0, z0

end

####

function getcoefflist(power,div,mom)

  # mom are the moments, mom[1] is M₁, etc

  if power%div!=0
    error("Indivisible power")
  end
  pow = Int(power/div)

  # Find the set of multi-indices I for which
  # sum(k(t)*t) = power/div. Each row of I corresponds
  # to a different multi-index in the set
  I = [1]
  for j = 2:pow
    I = madvance(I)
  end

  # Find the coefficient for 1/ζ^(pow-1) in the fhat expansion
  coeff = 0
  for j = 1:size(I,1)
    sumI = sum(I[j,:])
    fact = 1
    for l = 1:pow
      il = I[j,l]
      fact *= mom[l]^il/factorial(il)/l^il
    end
    coeff += fact*(-1)^sumI
  end
  return -coeff/(pow-1)

end

function madvance(P)
  # Given a set P of multi-indices
  # that satisfy a certain moment condition
  # |p|_1 = m (for all p in P), find the set Pplus1
  # that satisfies |p'|_1 = m+1 (for all p' in Pplus1)
  nP = size(P,1)
  m = size(P,2)
  Ppad = [ones(nP) P zeros(nP)]
  rows = 0
  Pplus1 = Int64[]
  for k = 1:nP
    # Loop through each index in the kth
    #  member of P, and find non-zero indices for shifting
    for j = m:-1:0
      if Ppad[k,j+1] == 1
        rows += 1
        Pplus1 = [Pplus1;zeros(Int64,1,m+2)]
        Pplus1[rows,1:m+2] = Ppad[k,1:m+2]
        Pplus1[rows,j+2] = Pplus1[rows,j+2]+1
        Pplus1[rows,j+1] = Pplus1[rows,j+1]-1
      end
    end
  end
  Pplus1 = Pplus1[:,2:end]
  Pplus1 = sortslices(Pplus1, dims=1)
  dP = sum(abs.([transpose(Pplus1[1,:]);diff(Pplus1,1)]),2)
  return Pplus1[dP[:].!=0,:]

end

#######

struct DabsQuad{T,N,NQ}
  β :: Vector{T}
  qdat :: Tuple{Array{T,2},Array{T,2}}
end

function DabsQuad(β::Vector{T},tol::T) where T
  n = length(β)
  nqpts = max(ceil(Int,-log10(tol)),2)
  qdat = qdata(β,nqpts)
  DabsQuad{T,n,nqpts}(β,qdat)
end

function DabsQuad(β::Vector{T},
                  qdat::Tuple{Array{T,2},Array{T,2}}) where T
  n = length(β)
  nqpts = size(qdat[1],1)
  DabsQuad{T,n,nqpts}(β,qdat)
end

function (I::DabsQuad{T,N,NQ})(ζ1::Vector{ComplexF64},ζ2::Vector{ComplexF64},
                        sing1::Vector{Int64},ζ::Vector{ComplexF64}) where {T,N,NQ}
#=
This integrates |z'(λ)| from `ζ1` to `ζ2` on the unit circle, using
Gauss-Jacobi quadrature, where
`sing1` contains the index of the prevertex in `ζ` that `ζ1` corresponds
to, or 0 if `ζ1` is not a prevertex. Note that `ζ2` cannot be a prevertex.
=#

   (qnode,qwght) = I.qdat

   argz = angle.(ζ)

   argζ1 = angle.(ζ1)
   argζ2 = angle.(ζ2)
   ang21 = angle.(ζ2./ζ1)

   bigargz = transpose(repeat(argz, 1, NQ))

   discont = (argζ2-argζ1).*ang21 .< 0
   argζ2[discont] += 2π*sign.(ang21[discont])

   if isempty(sing1)
     sing1 = zeros(size(ζ1))
   end
   result = zeros(Float64,size(ζ1))

   nontriv = findall(ζ1.!=ζ2)
    #tic()
   for k in nontriv
     ζ1k, ζ2k, arg1k, arg2k, sing1k =
          ζ1[k], ζ2[k], argζ1[k], argζ2[k], sing1[k]
     ζs = vcat(ζ[1:sing1k-1],ζ[sing1k+1:end])
     dist = min(1,2*minimum(abs.(ζs .- ζ1k))/abs(ζ2k-ζ1k))
     argr = arg1k + dist*(arg2k-arg1k)
     ind = ((sing1k+N) % (N+1)) + 1
     ζnd = 0.5*((argr-arg1k)*qnode[:,ind] .+ argr .+ arg1k)
     wt = 0.5*abs(argr-arg1k)*qwght[:,ind]
     θ = (repeat(ζnd, 1, N) .- bigargz .+ 2π).%(2π)
     θ[θ.>π] = 2π .- θ[θ.>π]
     terms = 2sin.(0.5θ)
     if !any(terms==0.0)
        if sing1k > 0
            terms[:,sing1k] ./= abs.(ζnd .- arg1k)
            wt .*= (0.5*abs.(argr-arg1k)).^I.β[sing1k]
        end
        result[k] = transpose(exp.(log.(terms)*I.β))*wt
        while dist < 1.0
            argl = argr
            ζl = exp(im*argl)
            dist = min(1,2*minimum(abs.(ζ .- ζl)/abs(ζl-ζ2k)))
            argr = argl + dist*(arg2k-argl)
            ζnd = 0.5*((argr-argl)*qnode[:,N+1] .+ argr .+ argl)
            wt = 0.5*abs(argr-argl)*qwght[:,N+1]
            #θ = hcat([(ζnd - argz[i] + 2π).%(2π) for i = 1:N]...)
            θ = (repeat(ζnd, 1, N) .- bigargz .+ 2π).%(2π)
            θ[θ.>π] = 2π .- θ[θ.>π]
            terms = 2sin.(0.5θ)
            result[k] += transpose(exp.(log.(terms)*I.β))*wt
        end
     end
   end
    #toc()
    #nothing
   return result
end

function (I::DabsQuad{T,N,NQ})(ζ1::Vector{ComplexF64},ζ2::Vector{ComplexF64},
                        sing1::Vector{Int64},ζ::Vector{ComplexF64},pow::Int) where {T,N,NQ}

#=
This integrates λ⁻ᵏz'(λ) (where k is `pow`) from `ζ1` to `ζ2` on the unit circle,
using Gauss-Jacobi quadrature, where
`sing1` contains the index of the prevertex in `ζ` that `ζ1` corresponds
to, or 0 if `ζ1` is not a prevertex. Note that `ζ2` cannot be a prevertex.
=#

   (qnode,qwght) = I.qdat

   argz = angle.(ζ)

   argζ1 = angle.(ζ1)
   argζ2 = angle.(ζ2)
   ang21 = angle.(ζ2./ζ1)

   bigargz = transpose(repeat(argz, 1, NQ))

   discont = (argζ2-argζ1).*ang21 .< 0
   argζ2[discont] += 2π*sign.(ang21[discont])

   if isempty(sing1)
     sing1 = zeros(size(ζ1))
   end
   result = zeros(ComplexF64,size(ζ1))

   nontriv = findall(ζ1.!=ζ2)
    #tic()
   for k in nontriv
     ζ1k, ζ2k, arg1k, arg2k, sing1k =
          ζ1[k], ζ2[k], argζ1[k], argζ2[k], sing1[k]
     ζs = vcat(ζ[1:sing1k-1],ζ[sing1k+1:end])
     dist = min(1,2*minimum(abs.(ζs .- ζ1k))/abs(ζ2k-ζ1k))
     argr = arg1k + dist*(arg2k-arg1k)
     ind = ((sing1k+N) % (N+1)) + 1
     ζnd = 0.5*((argr-arg1k)*qnode[:,ind] .+ argr .+ arg1k)
     wt = 0.5*abs(argr-arg1k)*qwght[:,ind]
     θ = (repeat(ζnd, 1, N) .- bigargz .+ 2π).%(2π)
     #θ[θ.>π] = 2π-θ[θ.>π]
     #terms = 2sin.(0.5θ)
     terms = 1 .- exp.(im*θ)
     if !any(terms==0.0)
        #terms = hcat(terms,exp.(im*ζnd))
        if sing1k > 0
            terms[:,sing1k] ./= abs.(ζnd .- arg1k)
            wt .*= (0.5*abs.(argr-arg1k)).^I.β[sing1k]
        end
        result[k] = transpose(exp.(log.(terms)*I.β+im*(pow-1)*ζnd))*wt
        while dist < 1.0
            argl = argr
            ζl = exp(im*argl)
            dist = min(1,2*minimum(abs.(ζ .- ζl)/abs(ζl-ζ2k)))
            argr = argl + dist*(arg2k-argl)
            ζnd = 0.5*((argr-argl)*qnode[:,N+1] .+ argr .+ argl)
            wt = 0.5*abs(argr-argl)*qwght[:,N+1]
            #θ = hcat([(ζnd - argz[i] + 2π).%(2π) for i = 1:N]...)
            θ = (repeat(ζnd, 1, N) .- bigargz .+ 2π).%(2π)
            #θ[θ.>π] = 2π-θ[θ.>π]
            #terms = 2sin.(0.5θ)
            terms = 1 .- exp.(im*θ)
            #terms = hcat(terms,exp.(im*ζnd))
            result[k] += transpose(exp.(log.(terms)*I.β+im*(pow-1)*ζnd))*wt
        end
     end
   end
    #toc()
    #nothing
   return result
end

struct DQuad{T,N,NQ}
  β :: Vector{T}
  qdat :: Tuple{Array{T,2},Array{T,2}}
end

function DQuad(β::Vector{T},tol::T) where T
  n = length(β)
  nqpts = max(ceil(Int,-log10(tol)),2)
  qdat = qdata(β,nqpts)
  DQuad{T,n,nqpts}(β,qdat)
end

function DQuad(β::Vector{T},
                  qdat::Tuple{Array{T,2},Array{T,2}}) where T
  n = length(β)
  nqpts = size(qdat[1],1)
  DQuad{T,n,nqpts}(β,qdat)
end

function (I::DQuad{T,N,NQ})(ζ1::Vector{ComplexF64},ζ2::Vector{ComplexF64},
          sing1::Vector{Int64},ζ::Vector{ComplexF64};pow::Int=0) where {T,N,NQ}
#=
This integrates z'(λ) from `ζ1` to `ζ2` on a straight path in the circle
plane, where `sing1` contains the index of the prevertex in `ζ` that `ζ1` corresponds
to, or 0 if `ζ1` is not a prevertex. Note that `ζ2` cannot be a prevertex.
If the optional argument `pow` is included, then it integrates λ⁻ᵏz'(λ), where
k = `pow`.
=#

   (qnode,qwght) = I.qdat

   β = [I.β;-2+pow]

   bigζ = transpose(repeat(ζ, 1, NQ))

   if isempty(sing1)
     sing1 = zeros(Int,size(ζ1))
   end
   result = zeros(ComplexF64,size(ζ1))

   nontriv = findall(ζ1.!=ζ2)
    #tic()
   for k in nontriv
     ζ1k, ζ2k, sing1k = ζ1[k], ζ2[k], sing1[k]
     ζs = vcat(ζ[1:sing1k-1],ζ[sing1k+1:end])
     dist = min(1,2*minimum(abs.(ζs .- ζ1k))/abs(ζ2k-ζ1k))
     ζr = ζ1k + dist*(ζ2k-ζ1k)
     # Choose which type of Gauss-Jacobi weights to use based on whether
     # ζ1k is a prevertex.
     ind = sing1k + (N+1)*(sing1k==0)

     ζnd = 0.5*((ζr-ζ1k)*qnode[:,ind] .+ ζr .+ ζ1k)
     wt = 0.5*(ζr-ζ1k)*qwght[:,ind]
     terms = 1 .- repeat(ζnd, 1, N)./bigζ
     if !any(terms==0.0)
       terms = hcat(terms,ζnd)
        # If ζ1k is a prevertex, adjust the integrand so that it is properly
        # set up for the Gauss-Jacobi integration
        if sing1k > 0
            terms[:,sing1k] ./= abs.(terms[:,sing1k])
            wt .*= (0.5*abs.(ζr-ζ1k)).^β[sing1k]
        end
        result[k] = transpose(exp.(log.(terms)*β))*wt

        while dist < 1.0
            ζl = ζr
            dist = min(1,2*minimum(abs.(ζ .- ζl)/abs(ζl-ζ2k)))
            ζr = ζl + dist*(ζ2k-ζl)
            ζnd = 0.5*((ζr-ζl)*qnode[:,N+1] .+ ζr .+ ζl)
            wt = 0.5*(ζr-ζl)*qwght[:,N+1]
            terms = 1 .- repeat(ζnd, 1, N)./bigζ
            terms = hcat(terms,ζnd)
            result[k] += transpose(exp.(log.(terms)*β))*wt
        end
      end
    end
   return result
end

struct Depfun{T,N,NQ}
    ζ :: Vector{Complex{T}}
    β :: Vector{T}
    nmlen :: Vector{T}
    qdat :: Tuple{Array{T,2},Array{T,2}}
    θ    :: Vector{T}
    mid  :: Vector{Complex{T}}
    ints :: Vector{T}
    dabsquad :: DabsQuad{T,N,NQ}
end

function Depfun(β::Vector{T},nmlen::Vector{T},qdat:: Tuple{Array{T,2},Array{T,2}}) where T
    # should compute nmlen in here
    n = length(β)
    nqpts = size(qdat[1],1)
    ζ = zeros(ComplexF64,n)
    θ = zeros(n-1)
    mid = zeros(ComplexF64,n-2)
    ints = zeros(ComplexF64,n-2)
    dabsquad = DabsQuad(β,qdat)
    Depfun{T,n,nqpts}(ζ,β,nmlen,qdat,θ,mid,ints,dabsquad)
end

function (R::Depfun{T,N,NQ})(F,y) where {T,N,NQ}

  y_to_ζ!(R.ζ,R.θ,y)

  @. R.mid = exp(im*0.5*(R.θ[1:N-2]+R.θ[2:N-1]))

  #tic()
  R.ints .= R.dabsquad(R.ζ[1:N-2],R.mid,collect(1:N-2),R.ζ)
  R.ints .+= R.dabsquad(R.ζ[2:N-1],R.mid,collect(2:N-1),R.ζ)

  #toc()

  if N > 3
    @. F[1:N-3] = abs(R.ints[2:N-2])/abs(R.ints[1]) - R.nmlen
  end

  res = -sum(R.β./R.ζ)/R.ints[1]
  F[N-2] = real(res)
  F[N-1] = imag(res)
end

function y_to_ζ!(ζ::Vector{Complex{T}},
                    θ::Vector{T},y::Vector{T}) where T

  cs = cumsum(cumprod([1;exp.(-y)]))
  n = length(cs)
  @. θ = 2π*cs[1:n-1]/cs[n]
  ζ[n] = 1.0
  @. ζ[1:n-1] = exp(im*θ)

end


text:
import Base: parse
"""
    Features(data::String)

Convert a string to morphological feature object.
"""
struct Features <: AbstractFeature
    data::String
end

"""
    Prefix(data::Symbol, pos::AbstractPartOfSpeech)

Create a new `Prefix` object with data as the symbol of the morphological feature with pos as
the corresponding Part of Speech.
"""
struct Prefix <: AbstractFeature
    data::Symbol
    pos::AbstractPartOfSpeech
end

"""
    Suffix(data::Symbol, pos::AbstractPartOfSpeech)

Create a new `Suffix` object with data as the symbol of the morphological feature with pos as
the corresponding Part of Speech.
"""
struct Suffix <: AbstractFeature
    data::Symbol
    pos::AbstractPartOfSpeech
    feats::Array{AbstractFeature,1}
end

"""
    Stem(data::Symbol, pos::AbstractPartOfSpeech)

Create a new `Stem` object with data as the symbol of the morphological feature with `pos` as
the corresponding Part of Speech.
"""
struct Stem <: AbstractFeature
    data::Symbol
    pos::AbstractPartOfSpeech
    feats::Array{AbstractFeature,1}
end

"""
    Lemma(data::String)

Convert a string to a `Lemma` object.
"""
struct Lemma <: AbstractFeature
    data::String
end

"""
    Root(data::String)

Convert a string to a `Root` object.
"""
struct Root <: AbstractFeature
    data::String
end

"""
    Special(data::String)

Convert a string to a `Special` object.
"""
struct Special <: AbstractFeature
    data::String
end

function derivednoun(strs::Array{String,1})
    idx_pcpl = findfirst(x -> occursin(r"PCPL", x), strs)
    if !isa(idx_pcpl, Nothing)
        pos_lab = Symbol(string(strs[idx_pcpl-1], " ", strs[idx_pcpl]))
        return PARTOFSPEECH[pos_lab]
    else
        return nothing
    end
end

function accusative(s::String, strs::Array{String,1}, other_feats::Array{AbstractFeature,1})
    isstate = in("DEF", strs) || in("INDEF", strs)
    return isstate ?
        push!(other_feats, PARTOFSPEECH[Symbol(s * "C")]) :
        push!(other_feats, PARTOFSPEECH[Symbol(s)])
end

function defaultverb(pos::AbstractPartOfSpeech, other_feats::Array{AbstractFeature,1})
    feats = vcat(pos, other_feats)
    out = findfirst(x -> x isa Verb, feats)
    if !isa(out, Nothing)
        mood = findfirst(x -> x isa AbstractMood, feats)
        if mood isa Nothing
            push!(other_feats, Indicative())
        end

        voice = findfirst(x -> x isa AbstractVoice, feats)
        if voice isa Nothing
            push!(other_feats, Active())
        end

        vform = findfirst(x -> x isa AbstractForm, feats) # change this to AbstractForm
        if vform isa Nothing
            push!(other_feats, VerbFormI())
        end
        return other_feats
    else
        return other_feats
    end
end

function stem(feat::Features)
    occursin("STEM", feat.data) ||
        throw(DomainError(feat, "Expected STEM feature, got " * string(split(feat.data, "|")[1]) * "."))
    strs = split(feat.data, "|")[2:end]

    pos = string(split(strs[1], ":")[2])
    other_feats = AbstractFeature[]

    idx_lemma = findfirst(x -> occursin(r"^LEM", x), strs)
    idx_root = findfirst(x -> occursin(r"^ROOT", x), strs)
    idx_sp = findfirst(x -> occursin(r"^SP", x), strs)
    idx_mood = findfirst(x -> occursin(r"^MOOD", x), strs)
    if !isa(idx_lemma, Nothing)
        lemma = Lemma(string(split(strs[idx_lemma], ":")[2]))
        push!(other_feats, lemma)
    end
    if !isa(idx_root, Nothing)
        root = Root(string(split(strs[idx_root], ":")[2]))
        push!(other_feats, root)
    end
    if !isa(idx_sp, Nothing)
        sp = Special(string(split(strs[idx_sp], ":")[2]))
        push!(other_feats, sp)
    end
    if !isa(idx_mood, Nothing)
        mood = PARTOFSPEECH[Symbol(string(split(strs[idx_mood], ":")[2]))]
        push!(other_feats, mood)
    end

    idx_others = findall(x -> !occursin(r"POS|LEM|ROOT|SP|MOOD", x), strs)
    drvnoun = derivednoun(string.(strs[idx_others]))

    if !isa(drvnoun, Nothing)
        push!(other_feats, drvnoun)
    end

    idx_others = findall(x -> !occursin(r"POS|LEM|ROOT|SP|MOOD|PCPL", x), strs)
    remaining = string.(strs[idx_others])
    for s in remaining
        try
            if s === "ACC"
                other_feats = accusative(s, remaining, other_feats)
            else
                push!(other_feats, PARTOFSPEECH[Symbol(s)])
            end
        catch
            for c in s
                if c === 'P'
                    push!(other_feats, PARTOFSPEECH[Symbol(c * "L")])
                else
                    push!(other_feats, PARTOFSPEECH[Symbol(c)])
                end
            end
        end
    end
    other_feats = defaultverb(PARTOFSPEECH[Symbol(pos)], other_feats)
    return Stem(Symbol(pos), PARTOFSPEECH[Symbol(pos)], other_feats)
end

function prefix(feat::Features)
    occursin("PREFIX", feat.data) ||
        throw(DomainError(feat, "Expected PREFIX feature, got " * string(split(feat.data, "|")[1]) * "."))
    strs = string(split(feat.data, "|")[2])
    return Prefix(Symbol(strs), PARTOFSPEECH[Symbol(strs)])
end

function suffix(feat::Features)
    occursin("SUFFIX", feat.data) ||
        throw(DomainError(feat, "Expected SUFFIX feature, got " * string(split(feat.data, "|")[1]) * "."))
    strs = string(split(feat.data, "|")[2])
    if occursin("PRON", strs)
        suff_feat = split(strs, ":")

        if length(suff_feat[2]) === 2
            if string(suff_feat[2][1]) === "P"
                person = Symbol(suff_feat[2][1] * "L")
            else
                person = Symbol(suff_feat[2][1])
            end
            number = Symbol(suff_feat[2][2])
            other_feats = [PARTOFSPEECH[person], PARTOFSPEECH[number]]
        else
            if string(suff_feat[2][1]) === "P"
                person = Symbol(suff_feat[2][1] * "L")
            else
                person = Symbol(suff_feat[2][1])
            end
            gender = Symbol(suff_feat[2][2])
            number = Symbol(suff_feat[2][3])
            other_feats = [PARTOFSPEECH[person], PARTOFSPEECH[gender], PARTOFSPEECH[number]]
        end
        return Suffix(Symbol(suff_feat[1]), PARTOFSPEECH[Symbol(suff_feat[1])], other_feats)
    else
        return Suffix(Symbol(strs), PARTOFSPEECH[Symbol(strs)], AbstractFeature[])
    end
end

"""
    parse(::Type{Features}, f::AbstractString)

Extract the features of a morphological `Feature` object.

# Examples
```julia-repl
julia> data = QuranData()
julia> crps, tnzl = load(data)
julia> crpsdata = table(crps)
julia> tnzldata = table(tnzl)
julia> parse(Features, select(crpsdata.data, :features)[53])
Stem(:NEG, NEG, AbstractFeature[Lemma("laA"), Special("<in~")])
```
"""
function parse(::Type{Features}, f::AbstractString)
    try
        return prefix(Features(f))
    catch
        try
            return stem(Features(f))
        catch
            return suffix(Features(f))
        end
    end
end

"""
    isfeature(feat::Features, pos::Type{<:AbstractFeature})

Check if the morphological `Feature` object is a type of `pos`.

# Examples
```julia-repl
julia> data = QuranData()
julia> crps, tnzl = load(data)
julia> crpsdata = table(crps)
julia> tnzldata = table(tnzl)
julia> isfeature(parse(Features, select(crpsdata[1].data, :features)[2]), Stem)
true
```
"""
function isfeature(feat::AbstractFeature, pos::Type{<:AbstractFeature})
    try
        feats = vcat(feat, feat.pos, feat.feats)
        out = findfirst(x -> x isa pos, feats)
        return out isa Nothing ? false : true
    catch
        feats = vcat(feat, feat.pos)
        out = findfirst(x -> x isa pos, feats)
        return out isa Nothing ? false : true
    end
end

"""
    root(feat::AbstractFeature)

Extract the root of the feature.

# Examples
```julia-repl
julia> data = QuranData()
julia> crps, tnzl = load(data)
julia> crpsdata = table(crps)
julia> tnzldata = table(tnzl)
julia> root(parse(Features, select(crpsdata[112].data, :features)[1]))
"qwl"
```
"""
function root(feat::AbstractFeature)
    try
        idx = findfirst(x -> x isa Root, feat.feats)
        return idx isa Nothing ? missing : feat.feats[idx].data
    catch
        return missing
    end
end

function parse(::Type{Lemma}, l::AbstractString)
    idx = findfirst(r"\d", l)
    if idx isa Nothing
        return Lemma(l)
    else
        out = AbstractFeature[]
        push!(out, Lemma(l[1:idx.start-1]))
        push!(out, PARTOFSPEECH[Symbol(l[idx])])
        return out
    end
end

"""
    lemma(feat::AbstractFeature)

Extract the lemma of the feature.

# Examples
```julia-repl
julia> data = QuranData()
julia> crps, tnzl = load(data)
julia> crpsdata = table(crps)
julia> tnzldata = table(tnzl)
julia> lemma(parse(Features, select(crpsdata[112].data, :features)[1]))
"qaAla"
```
"""
function lemma(feat::AbstractFeature)
    try
        idx = findfirst(x -> x isa Lemma, feat.feats)
        if idx isa Nothing
            return missing
        else
            lem = parse(Lemma, feat.feats[idx].data)
            if lem isa Lemma
                lem = lem.data
            else
                lem[1].data
            end
        end
    catch
        return missing
    end
end

"""
    special(feat::AbstractFeature)

Extract the special feature of the token.

# Examples
```julia-repl
julia> data = QuranData()
julia> crps, tnzl = load(data)
julia> crpsdata = table(crps)
julia> tnzldata = table(tnzl)
julia> special(parse(Features, select(crpsdata.data, :features)[53]))
"<in~"
```
"""
function special(feat::AbstractFeature)
    try
        idx = findfirst(x -> x isa Special, feat.feats)
        return idx isa Nothing ? missing : feat.feats[idx].data
    catch
        return missing
    end
end

ext: jl


text:
include("IOM/src/share/constants.jl")
include("IOM/src/share/PolelikeCoordinate.jl")
include("IOM/src/share/BasicMatrixOperators.jl")
include("IOM/src/share/AdvancedMatrixOperators.jl")

using .PolelikeCoordinate
using NCDatasets
using ArgParse, JSON

output_file = "output.nc"

# Load domain and construct operators
domain_file = "CESM_domains/domain.ocn.gx1v6.090206.nc"
gf = PolelikeCoordinate.CurvilinearSphericalGridFile(
    domain_file;
    R  = Re,
    Ω  = Ω,
)

mask_sT = reshape(gf.mask, 1, size(gf.mask)...)
gd_slab = PolelikeCoordinate.genGrid(gf, [0, -50.0])

println("Constructing operators...")
@time amo_slab = AdvancedMatrixOperators(;
    gd     = gd_slab,
    mask_T     = mask_sT,
    deepmask_T = mask_sT,
)

u0 = 1.0
v0 = 1.0
a = 2.0
b = 3.0

ϕ_sT = gd_slab.ϕ_T
λ_sT = gd_slab.λ_T
u_east_sT  = u0 * cos.(ϕ_sT) .* sin.(a * λ_sT)
v_north_sT = v0 * cos.(ϕ_sT).^2 .* cos.(b * λ_sT)

ζ_ana_sT = 1 / gd_slab.R * ( 2 * u0 * sin.(ϕ_sT) .* sin.(a * λ_sT) - v0 * b * cos.(ϕ_sT) .* sin.(b * λ_sT) )

u_sT, v_sT = PolelikeCoordinate.project(
    gd_slab,
    u_east_sT,
    v_north_sT,
    direction=:Forward,
    grid=:T,
)

ζ_num_sT = reshape(
    amo_slab.T_CURLx_T * view(v_sT, :) + amo_slab.T_CURLy_T * view(u_sT, :),
    gd_slab.Nx, gd_slab.Ny,
)

ζ_ana_sT = reshape(
    amo_slab.T_bordermask_T * view(ζ_ana_sT, :),
    gd_slab.Nx, gd_slab.Ny,
)

ζ_dif_sT = ζ_num_sT - ζ_ana_sT
ζ_dif_sT[gf.mask .== 0.0] .= NaN

Dataset(output_file, "c") do ds

    defDim(ds, "Nx", gd_slab.Nx)
    defDim(ds, "Ny", gd_slab.Ny)
    for (varname, vardata, vardim, attrib) in [
        ("zeta_num", ζ_num_sT, ("Nx", "Ny",), Dict()),
        ("zeta_ana", ζ_ana_sT, ("Nx", "Ny",), Dict()),
        ("zeta_dif", ζ_dif_sT, ("Nx", "Ny",), Dict()),

    ]

        var = defVar(ds, varname, Float64, vardim)
        var.attrib["_FillValue"] = 1e20

        for (k, v) in attrib
            var.attrib[k] = v
        end

        rng = []
        for i in 1:length(vardim)-1
            push!(rng, Colon())
        end
        push!(rng, 1:size(vardata)[end])
        var[rng...] = vardata

    end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
# function add_ref_dcgrid!(ref::Dict{Symbol,<:Any}, data::Dict{String,<:Any})
#     for (n, nw_ref) in ref[:nw]
#         if haskey(nw_ref, :convdc)
#             #Filter converters & DC branches with status 0 as well as wrong bus numbers
#             nw_ref[:convdc] = Dict([x for x in nw_ref[:convdc] if (x.second["status"] == 1 && x.second["busdc_i"] in keys(nw_ref[:busdc]) && x.second["busac_i"] in keys(nw_ref[:bus]))])
#             nw_ref[:branchdc] = Dict([x for x in nw_ref[:branchdc] if (x.second["status"] == 1 && x.second["fbusdc"] in keys(nw_ref[:busdc]) && x.second["tbusdc"] in keys(nw_ref[:busdc]))])
#
#             # DC grid arcs for DC grid branches
#             nw_ref[:arcs_dcgrid_from] = [(i,branch["fbusdc"],branch["tbusdc"]) for (i,branch) in nw_ref[:branchdc]]
#             nw_ref[:arcs_dcgrid_to]   = [(i,branch["tbusdc"],branch["fbusdc"]) for (i,branch) in nw_ref[:branchdc]]
#             nw_ref[:arcs_dcgrid] = [nw_ref[:arcs_dcgrid_from]; nw_ref[:arcs_dcgrid_to]]
#             nw_ref[:arcs_conv_acdc] = [(i,conv["busac_i"],conv["busdc_i"]) for (i,conv) in nw_ref[:convdc]]
#             #bus arcs of the DC grid
#             bus_arcs_dcgrid = Dict([(bus["busdc_i"], []) for (i,bus) in nw_ref[:busdc]])
#             for (l,i,j) in nw_ref[:arcs_dcgrid]
#                 push!(bus_arcs_dcgrid[i], (l,i,j))
#             end
#             nw_ref[:bus_arcs_dcgrid] = bus_arcs_dcgrid
#
#                 # bus_convs for AC side power injection of DC converters
#             bus_convs_ac = Dict([(i, []) for (i,bus) in nw_ref[:bus]])
#             for (i,conv) in nw_ref[:convdc]
#                 push!(bus_convs_ac[conv["busac_i"]], i)
#             end
#             nw_ref[:bus_convs_ac] = bus_convs_ac
#
#             # bus_convs for AC side power injection of DC converters
#             bus_convs_dc = Dict([(bus["busdc_i"], []) for (i,bus) in nw_ref[:busdc]])
#             for (i,conv) in nw_ref[:convdc]
#                 push!(bus_convs_dc[conv["busdc_i"]], i)
#             end
#             nw_ref[:bus_convs_dc] = bus_convs_dc
#             # Add DC reference buses
#             ref_buses_dc = Dict{String, Any}()
#             for (k,v) in nw_ref[:convdc]
#                 if v["type_dc"] == 2
#                     ref_buses_dc["$k"] = v
#                 end
#             end
#
#             if length(ref_buses_dc) == 0
#                 for (k,v) in nw_ref[:convdc]
#                     if v["type_ac"] == 2
#                         ref_buses_dc["$k"] = v
#                     end
#                 end
#                 Memento.warn(_PM._LOGGER, "no reference DC bus found, setting reference bus based on AC bus type")
#             end
#
#             for (k,conv) in nw_ref[:convdc]
#                 conv_id = conv["index"]
#                 if conv["type_ac"] == 2 && conv["type_dc"] == 1
#                     Memento.warn(_PM._LOGGER, "For converter $conv_id is chosen P is fixed on AC and DC side. This can lead to infeasibility in the PF problem.")
#                 elseif conv["type_ac"] == 1 && conv["type_dc"] == 1
#                     Memento.warn(_PM._LOGGER, "For converter $conv_id is chosen P is fixed on AC and DC side. This can lead to infeasibility in the PF problem.")
#                 end
#                 convbus_ac = conv["busac_i"]
#                 if conv["Vmmax"] < nw_ref[:bus][convbus_ac]["vmin"]
#                     Memento.warn(_PM._LOGGER, "The maximum AC side voltage of converter $conv_id is smaller than the minimum AC bus voltage")
#                 end
#                 if conv["Vmmin"] > nw_ref[:bus][convbus_ac]["vmax"]
#                     Memento.warn(_PM._LOGGER, "The miximum AC side voltage of converter $conv_id is larger than the maximum AC bus voltage")
#                 end
#             end
#
#             if length(ref_buses_dc) > 1
#                 ref_buses_warn = ""
#                 for (rb) in keys(ref_buses_dc)
#                     ref_buses_warn = ref_buses_warn*rb*", "
#                 end
#                 Memento.warn(_PM._LOGGER, "multiple reference buses found, i.e. "*ref_buses_warn*"this can cause infeasibility if they are in the same connected component")
#             end
#
#
#             # if haskey(pm.setting, "find_all_ac_grids") && pm.tetting["find_all_ac_grids"] == true
#             #     ACgrids = find_all_ac_grids(nw_ref[:branch], nw_ref[:bus])
#             #     for (i, grid) in ACgrids
#             #         a = 0
#             #         for (j, bus) in nw_ref[:ref_buses]
#             #             if (bus["bus_i"] in grid["Buses"])
#             #                 a = 1
#             #             end
#             #         end
#             #         if a == 0
#             #             Memento.warn(_PM._LOGGER, "Grid $i does not have any voltage reference bus, this might cause infeasibility")
#             #         end
#             #     end
#             # end
#             nw_ref[:ref_buses_dc] = ref_buses_dc
#             nw_ref[:buspairsdc] = buspair_parameters_dc(nw_ref[:arcs_dcgrid_from], nw_ref[:branchdc], nw_ref[:busdc])
#         else
#             nw_ref[:convdc] = Dict{String, Any}()
#             nw_ref[:busdc] = Dict{String, Any}()
#             nw_ref[:branchdc] = Dict{String, Any}()
#             # DC grid arcs for DC grid branches
#             nw_ref[:arcs_dcgrid] = Dict{String, Any}()
#             nw_ref[:arcs_dcgrid_from] = Dict{String, Any}()
#             nw_ref[:arcs_dcgrid_to] = Dict{String, Any}()
#             nw_ref[:arcs_conv_acdc] = Dict{String, Any}()
#             nw_ref[:bus_arcs_dcgrid] = Dict{String, Any}()
#             nw_ref[:arcs_reserves_syn] = Dict{String, Any}()
#             bus_convs_ac = Dict([(i, []) for (i,bus) in nw_ref[:bus]])
#             for (i,conv) in nw_ref[:convdc]
#                 push!(bus_convs_ac[conv["busac_i"]], i)
#             end
#             nw_ref[:bus_convs_ac] = bus_convs_ac
#             nw_ref[:bus_convs_dc] = Dict{String, Any}()
#             nw_ref[:ref_buses_dc] = Dict{String, Any}()
#             nw_ref[:buspairsdc] = Dict{String, Any}()
#         end
#     end
# end
function add_ref_dcgrid!(ref::Dict{Symbol,<:Any}, data::Dict{String,<:Any})
    for (n, nw_ref) in ref[:nw]
        if haskey(nw_ref, :convdc)
            #Filter converters & DC branches with status 0 as well as wrong bus numbers
            nw_ref[:convdc] = Dict([x for x in nw_ref[:convdc] if (x.second["status"] == 1 && x.second["busdc_i"] in keys(nw_ref[:busdc]) && x.second["busac_i"] in keys(nw_ref[:bus]))])
            nw_ref[:branchdc] = Dict([x for x in nw_ref[:branchdc] if (x.second["status"] == 1 && x.second["fbusdc"] in keys(nw_ref[:busdc]) && x.second["tbusdc"] in keys(nw_ref[:busdc]))])
            # DC grid arcs for DC grid branches
            nw_ref[:arcs_dcgrid_from] = [(i,branch["fbusdc"],branch["tbusdc"]) for (i,branch) in nw_ref[:branchdc]]
            nw_ref[:arcs_dcgrid_to]   = [(i,branch["tbusdc"],branch["fbusdc"]) for (i,branch) in nw_ref[:branchdc]]
            nw_ref[:arcs_dcgrid] = [nw_ref[:arcs_dcgrid_from]; nw_ref[:arcs_dcgrid_to]]
            nw_ref[:arcs_conv_acdc] = [(i,conv["busac_i"],conv["busdc_i"]) for (i,conv) in nw_ref[:convdc]]
            #bus arcs of the DC grid
            bus_arcs_dcgrid = Dict([(bus["busdc_i"], []) for (i,bus) in nw_ref[:busdc]])
            for (l,i,j) in nw_ref[:arcs_dcgrid]
                push!(bus_arcs_dcgrid[i], (l,i,j))
            end

            nw_ref[:bus_arcs_dcgrid] = bus_arcs_dcgrid
        # for OPF OPF OPF OPF##########

        bus_convs_dirct =  Dict([(reserve["syncarea"], []) for (i,reserve) in nw_ref[:reserves]])
        #to determine which converters are connected to particular dc line contingency
        for (b,br) in nw_ref[:branchdc]
            if br["rateA"] == 0
                for (c, cv) in nw_ref[:convdc]
                    if cv["busdc_i"] == br["fbusdc"] ||  cv["busdc_i"] == br["tbusdc"]
                        push!(bus_convs_dirct[cv["syncarea"]], c)
                    end
                end
            end
        end

        nw_ref[:bus_convs_dirct] =   bus_convs_dirct
        display("bus_convs_dirct:$n")
        display(bus_convs_dirct)

            # display("nw)ref")
            # display(nw_ref)
            if haskey(nw_ref, :reserves)
            # display("Multinetwork:reserves")

            nw_ref[:arcs_sync_conv] = [(i,conv["busac_i"],conv["syncarea"]) for (i,conv) in nw_ref[:convdc]]
            # display("arcs_sync_conv")
            # display(nw_ref[:load])
            nw_ref[:reserves] = Dict([x for x in nw_ref[:reserves]])

            bus_arcs_load = Dict([(reserve["syncarea"], []) for (i,reserve) in nw_ref[:reserves]])

            for (la,load) in nw_ref[:load]
                    ac_bus = load["load_bus"]
                    # display("la:$la")
                    # display("ac_bus:$ac_bus")
                    for (c,conv) in nw_ref[:convdc]
                        # display("c:$c")
                        if conv["busac_i"] == ac_bus
                            # display("busac_i")
                            # display(conv["busac_i"])
                            # display("syncarea")
                            # display(conv["syncarea"])
                            # display("la")
                            # display(la)
                        push!(bus_arcs_load[conv["syncarea"]], la)
                        end
                    end
                end

            nw_ref[:bus_arcs_load] = bus_arcs_load
            # display("bus_arcs_load")
            #
            # display(bus_arcs_load)

            bus_arcs_conv = Dict([(reserve["syncarea"], []) for (i,reserve) in nw_ref[:reserves]])
            for (l,i,j) in nw_ref[:arcs_sync_conv]
                push!(bus_arcs_conv[j], l)
            end
            nw_ref[:bus_arcs_conv] = bus_arcs_conv
            # display("bus_arcs_conv")
            # display(bus_arcs_conv)

            bus_arcs_sync = Dict([(reserve["syncarea"], []) for (i,reserve) in nw_ref[:reserves]])
            for (l,i,j) in nw_ref[:arcs_sync_conv]
                push!(bus_arcs_sync[j], i)
            end
            nw_ref[:bus_arcs_sync] = bus_arcs_sync
            # display(bus_arcs_sync)

            arcs_reserves_syn = Dict([(i, []) for (i,reserve) in nw_ref[:reserves]])
            for (i,reserve) in nw_ref[:reserves]
                push!(arcs_reserves_syn[reserve["syncarea"]], i)
            end
            nw_ref[:arcs_reserves_syn] = arcs_reserves_syn
            # display("arcs_reserves_syn")
            # display(arcs_reserves_syn)

            bus_syc_arc = Dict([(i, []) for (i,bus) in nw_ref[:bus]])
            for (i,conv) in nw_ref[:convdc]
                 # display(conv)
                 push!(bus_syc_arc[conv["busac_i"]],conv["syncarea"])
            end
            # display(bus_syc_arc)
        end
    # end of  of OPF OPF OPF OPF

             # bus_convs for AC side power injection of DC converters
            bus_convs_ac = Dict([(i, []) for (i,bus) in nw_ref[:bus]])
            for (i,conv) in nw_ref[:convdc]
                push!(bus_convs_ac[conv["busac_i"]], i)
            end
            nw_ref[:bus_convs_ac] = bus_convs_ac

            # bus_convs for AC side power injection of DC converters
            bus_convs_dc = Dict([(bus["busdc_i"], []) for (i,bus) in nw_ref[:busdc]])
            for (i,conv) in nw_ref[:convdc]
                push!(bus_convs_dc[conv["busdc_i"]], i)
            end
            nw_ref[:bus_convs_dc] = bus_convs_dc
            # Add DC reference buses
            ref_buses_dc = Dict{String, Any}()
            for (k,v) in nw_ref[:convdc]
                if v["type_dc"] == 2
                    ref_buses_dc["$k"] = v
                end
            end

            if length(ref_buses_dc) == 0
                for (k,v) in nw_ref[:convdc]
                    if v["type_ac"] == 2
                        ref_buses_dc["$k"] = v
                    end
                end
                Memento.warn(_PM._LOGGER, "no reference DC bus found, setting reference bus based on AC bus type")
            end

            for (k,conv) in nw_ref[:convdc]
                conv_id = conv["index"]
                if conv["type_ac"] == 2 && conv["type_dc"] == 1
                    Memento.warn(_PM._LOGGER, "For converter $conv_id is chosen P is fixed on AC and DC side. This can lead to infeasibility in the PF problem.")
                elseif conv["type_ac"] == 1 && conv["type_dc"] == 1
                    Memento.warn(_PM._LOGGER, "For converter $conv_id is chosen P is fixed on AC and DC side. This can lead to infeasibility in the PF problem.")
                end
                convbus_ac = conv["busac_i"]
                if conv["Vmmax"] < nw_ref[:bus][convbus_ac]["vmin"]
                    Memento.warn(_PM._LOGGER, "The maximum AC side voltage of converter $conv_id is smaller than the minimum AC bus voltage")
                end
                if conv["Vmmin"] > nw_ref[:bus][convbus_ac]["vmax"]
                    Memento.warn(_PM._LOGGER, "The miximum AC side voltage of converter $conv_id is larger than the maximum AC bus voltage")
                end
            end

            if length(ref_buses_dc) > 1
                ref_buses_warn = ""
                for (rb) in keys(ref_buses_dc)
                    ref_buses_warn = ref_buses_warn*rb*", "
                end
                Memento.warn(_PM._LOGGER, "multiple reference buses found, i.e. "*ref_buses_warn*"this can cause infeasibility if they are in the same connected component")
            end

            # if haskey(pm.setting, "find_all_ac_grids") && pm.tetting["find_all_ac_grids"] == true
            #     ACgrids = find_all_ac_grids(nw_ref[:branch], nw_ref[:bus])
            #     for (i, grid) in ACgrids
            #         a = 0
            #         for (j, bus) in nw_ref[:ref_buses]
            #             if (bus["bus_i"] in grid["Buses"])
            #                 a = 1
            #             end
            #         end
            #         if a == 0
            #             Memento.warn(_PM._LOGGER, "Grid $i does not have any voltage reference bus, this might cause infeasibility")
            #         end
            #     end
            # end
            nw_ref[:ref_buses_dc] = ref_buses_dc
            nw_ref[:buspairsdc] = buspair_parameters_dc(nw_ref[:arcs_dcgrid_from], nw_ref[:branchdc], nw_ref[:busdc])
        else
            nw_ref[:convdc] = Dict{String, Any}()
            nw_ref[:busdc] = Dict{String, Any}()
            nw_ref[:branchdc] = Dict{String, Any}()
            nw_ref[:arcs_reserves_syn] = Dict{String, Any}()
            # DC grid arcs for DC grid branches
            nw_ref[:arcs_dcgrid] = Dict{String, Any}()
            nw_ref[:arcs_dcgrid_from] = Dict{String, Any}()
            nw_ref[:arcs_dcgrid_to] = Dict{String, Any}()
            nw_ref[:arcs_conv_acdc] = Dict{String, Any}()
            nw_ref[:bus_arcs_dcgrid] = Dict{String, Any}()
            bus_convs_ac = Dict([(i, []) for (i,bus) in nw_ref[:bus]])
            for (i,conv) in nw_ref[:convdc]
                push!(bus_convs_ac[conv["busac_i"]], i)
            end
            nw_ref[:bus_convs_ac] = bus_convs_ac
            nw_ref[:bus_convs_dc] = Dict{String, Any}()
            nw_ref[:ref_buses_dc] = Dict{String, Any}()
            nw_ref[:buspairsdc] = Dict{String, Any}()
        end
    end
end

"compute bus pair level structures"
function buspair_parameters_dc(arcs_dcgrid_from, branches, buses)
    buspair_indexes = collect(Set([(i,j) for (l,i,j) in arcs_dcgrid_from]))

    bp_branch = Dict([(bp, Inf) for bp in buspair_indexes])

    for (l,branch) in branches
        i = branch["fbusdc"]
        j = branch["tbusdc"]

        bp_branch[(i,j)] = min(bp_branch[(i,j)], l)
    end

    buspairs = Dict([((i,j), Dict(
    "branch"=>bp_branch[(i,j)],
    "vm_fr_min"=>buses[i]["Vdcmin"],
    "vm_fr_max"=>buses[i]["Vdcmax"],
    "vm_to_min"=>buses[j]["Vdcmin"],
    "vm_to_max"=>buses[j]["Vdcmax"]
    )) for (i,j) in buspair_indexes])

    return buspairs
end

function find_all_ac_grids(branches_ac, buses_ac)
    ACgrids = Dict{String, Any}()

    if isempty(branches_ac)
        for (i, bus) in buses_ac
            ACgrids["$i"] = Dict{String, Any}()
            ACgrids["$i"]["Buses"] = bus["index"]
        end
    else
        ACgrids["1"] = Dict{String, Any}()
        ACgrids["1"]["Buses"] = [branches_ac[1]["f_bus"] branches_ac[1]["t_bus"]]
        closed_buses = [branches_ac[1]["f_bus"] branches_ac[1]["t_bus"]]
        closed_branches = [1]
        connections = []
        buses = []
        for (i, bus) in buses_ac
            if VERSION < v"0.7.0-"
                buses = cat(1,buses,bus["index"])
            else
                buses = cat(buses, bus["index"], dims = 1)
            end
        end
        grid_id = 1
        iter_id = 1
        branch_iter = 1
        while length(closed_buses) != length(buses) && iter_id < 10
            while branch_iter <= length(branches_ac)
                for (i, branch) in branches_ac
                    for (index, grid) in ACgrids
                        if (branch["t_bus"] in grid["Buses"]) && (branch["f_bus"] in grid["Buses"])
                            if !(branch["index"] in closed_branches)
                                closed_branches = [closed_branches branch["index"]]
                            end
                        elseif (branch["f_bus"] in grid["Buses"])
                            if !(branch["t_bus"] in grid["Buses"])
                                ACgrids["$index"]["Buses"] = [grid["Buses"] branch["t_bus"]]
                                closed_buses = [closed_buses branch["t_bus"]]
                                closed_branches = [closed_branches branch["index"]]
                            end
                        elseif (branch["t_bus"] in grid["Buses"])
                            if !(branch["f_bus"] in grid["Buses"])
                                ACgrids["$index"]["Buses"] = [grid["Buses"] branch["f_bus"]]
                                closed_buses = [closed_buses branch["f_bus"]]
                                closed_branches = [closed_branches branch["index"]]
                            end
                        end
                    end
                end
                branch_iter = branch_iter + 1
            end
            if length(closed_branches) < length(branches_ac)
                grid_id = grid_id + 1
                branch_iter = 1
                ACgrids["$grid_id"] = Dict{String, Any}()
                for (i, branch) in branches_ac
                    if !(branch["index"] in closed_branches) && isempty(ACgrids["$grid_id"])
                        ACgrids["$grid_id"]["Buses"] = [branch["f_bus"] branch["t_bus"]]
                        closed_branches = [closed_branches branch["index"]]
                    end
                end
            end
            iter_id = iter_id + 1 # to avoid infinite loop -> if not all subgrids detected
        end
    end
    return ACgrids
end

function add_candidate_dcgrid!(ref::Dict{Symbol,<:Any}, data::Dict{String,<:Any})
    for (n, nw_ref) in ref[:nw]
        if haskey(nw_ref, :convdc_ne)
            nw_ref[:arcs_dcgrid_from_ne] = [(i,branch["fbusdc"],branch["tbusdc"]) for (i,branch) in nw_ref[:branchdc_ne]]
            nw_ref[:arcs_dcgrid_to_ne]   = [(i,branch["tbusdc"],branch["fbusdc"]) for (i,branch) in nw_ref[:branchdc_ne]]
            nw_ref[:arcs_dcgrid_ne] = [nw_ref[:arcs_dcgrid_from_ne]; nw_ref[:arcs_dcgrid_to_ne]]
            nw_ref[:arcs_conv_acdc_ne] = [(i,conv["busac_i"],conv["busdc_i"]) for (i,conv) in nw_ref[:convdc_ne]]
            nw_ref[:arcs_conv_acdc_acbus_ne] = [(i,conv["busac_i"]) for (i,conv) in nw_ref[:convdc_ne]]
            nw_ref[:arcs_incident_ne] = [(i, incident["convdc"]) for (i,incident) in nw_ref[:incident]]

            #bus arcs of the DC grid
            bus_arcs_dcgrid_ne = Dict([(bus["busdc_i"], []) for (i,bus) in nw_ref[:busdc_ne]])
            for (l,i,j) in nw_ref[:arcs_dcgrid_ne]
                if haskey(bus_arcs_dcgrid_ne, i)
                    push!(bus_arcs_dcgrid_ne[i], (l,i,j))
                elseif  haskey(nw_ref[:bus_arcs_dcgrid], i)
                    bus_arcs_dcgrid_ne[i] = []
                    push!(bus_arcs_dcgrid_ne[i], (l,i,j))
                else
                    bus_arcs_dcgrid_ne[i] = []
                end
            end
            nw_ref[:bus_arcs_dcgrid_ne] = bus_arcs_dcgrid_ne
            # bus_convs for AC side power injection of DC converters
            bus_convs_ac = Dict([(i, []) for (i,bus) in nw_ref[:bus]])
            for (i,conv) in nw_ref[:convdc_ne]
                push!(bus_convs_ac[conv["busac_i"]], i)
            end
            nw_ref[:bus_convs_ac_ne] = bus_convs_ac

            # add new converters to existting DC buses
            bus_convs_dc_ne = Dict([(bus["busdc_i"], []) for (i,bus) in nw_ref[:busdc]])
            for (i,conv) in nw_ref[:convdc_ne]
                if haskey(bus_convs_dc_ne, conv["busdc_i"])
                    push!(bus_convs_dc_ne, i)
                end
            end
            nw_ref[:bus_convs_dc_ne] = bus_convs_dc_ne
            # Bus converters for candidate DC buses
            bus_ne_convs_dc_ne = Dict([(bus["busdc_i"], []) for (i,bus) in nw_ref[:busdc_ne]])
            for (i,conv) in nw_ref[:convdc_ne]
                if haskey(bus_ne_convs_dc_ne, conv["busdc_i"])
                    push!(bus_ne_convs_dc_ne[conv["busdc_i"]], i)
                end
            end
            nw_ref[:bus_ne_convs_dc_ne] = bus_ne_convs_dc_ne
            nw_ref[:ref_buses_dc_ne] = Dict{String, Any}()
            nw_ref[:buspairsdc_ne] = buspair_parameters_dc_ne(nw_ref[:arcs_dcgrid_from_ne], nw_ref[:branchdc_ne], nw_ref[:busdc_ne], nw_ref[:busdc])
        else
            nw_ref[:convdc_ne] = Dict{String, Any}()
            nw_ref[:busdc_ne] = Dict{String, Any}()
            nw_ref[:branchdc_ne] = Dict{String, Any}()
            # DC grid arcs for DC grid branches
            nw_ref[:arcs_dcgrid_ne] = Dict{String, Any}()
            nw_ref[:arcs_conv_acdc_ne] = Dict{String, Any}()
            nw_ref[:bus_arcs_dcgrid_ne] = Dict{String, Any}()
            bus_convs_ac = Dict([(i, []) for (i,bus) in nw_ref[:bus]])
            for (i,conv) in nw_ref[:convdc_ne]
                push!(bus_convs_ac[conv["busac_i"]], i)
            end
            nw_ref[:bus_convs_ac_ne] = bus_convs_ac
            nw_ref[:bus_convs_dc_ne] = Dict{String, Any}()
            nw_ref[:ref_buses_dc_ne] = Dict{String, Any}()
            nw_ref[:buspairsdc_ne] = Dict{String, Any}()
        end
    end
end

function add_ref_frequency!(ref::Dict{Symbol,<:Any}, data::Dict{String,<:Any})
    for (n, nw_ref) in ref[:nw]
        nw_ref[:reserves] = Dict([x for x in nw_ref[:reserves]])

        bus_convs_dirct =  Dict([(reserve["syncarea"], []) for (i,reserve) in nw_ref[:reserves]])
            #to determine which converters are connected to particular dc line contingency - FS strategies
            for (b,br) in nw_ref[:branchdc]
                if br["rateA"] == 0
                    for (c, cv) in nw_ref[:convdc]
                        if cv["busdc_i"] == br["fbusdc"] ||  cv["busdc_i"] == br["tbusdc"]
                            push!(bus_convs_dirct[cv["syncarea"]], c)
                        end
                    end
                end
            end
         nw_ref[:bus_convs_dirct] =   bus_convs_dirct

         if isempty(bus_convs_dirct) == true
             Memento.error(_LOGGER, "no contigency defined for frequency constraint. set rateA to zero for the line that has contigency")
         end

        nw_ref[:arcs_sync_conv] = [(i,conv["busac_i"],conv["syncarea"]) for (i,conv) in nw_ref[:convdc]]

        bus_arcs_load = Dict([(reserve["syncarea"], []) for (i,reserve) in nw_ref[:reserves]])
         for (la,load) in nw_ref[:load]
            ac_bus = load["load_bus"]
              for (c,conv) in nw_ref[:convdc]
                 if conv["busac_i"] == ac_bus
                 push!(bus_arcs_load[conv["syncarea"]], la)
                 end
             end
         end
         nw_ref[:bus_arcs_load] = bus_arcs_load

         bus_arcs_conv = Dict([(reserve["syncarea"], []) for (i,reserve) in nw_ref[:reserves]])
         for (l,i,j) in nw_ref[:arcs_sync_conv]
             push!(bus_arcs_conv[j], l)
         end
         nw_ref[:bus_arcs_conv] = bus_arcs_conv

         #which reserves belongs to which synchronous area
         arcs_reserves_syn = Dict([(i, []) for (i,reserve) in nw_ref[:reserves]])
         for (i,reserve) in nw_ref[:reserves]
           push!(arcs_reserves_syn[reserve["syncarea"]], i)
         end
         nw_ref[:arcs_reserves_syn] = arcs_reserves_syn
     end
end

function add_ref_frequency_candidates!(ref::Dict{Symbol,<:Any}, data::Dict{String,<:Any})
    for (n, nw_ref) in ref[:nw]
    bus_convs_dirct_ne =  Dict([(reserve["syncarea"], []) for (i,reserve) in nw_ref[:reserves]])
    #to determine which converters are connected to particular dc line contingency
    for (b,br) in nw_ref[:branchdc_ne]
        if br["rateA"] == 0
            for (c, cv) in nw_ref[:convdc_ne]
                if cv["busdc_i"] == br["fbusdc"] ||  cv["busdc_i"] == br["tbusdc"]
                    push!(bus_convs_dirct_ne[cv["syncarea"]], c)
                end
            end
        end
    end
        nw_ref[:bus_convs_dirct_ne] =   bus_convs_dirct_ne

                nw_ref[:arcs_sync_conv_ne] = [(i,conv["busac_i"],conv["syncarea"]) for (i,conv) in nw_ref[:convdc_ne]]
                nw_ref[:reserves] = Dict([x for x in nw_ref[:reserves]])
                bus_arcs_load_ne = Dict([(reserve["syncarea"], []) for (i,reserve) in nw_ref[:reserves]])
                for (la,load) in nw_ref[:load]
                        ac_bus = load["load_bus"]
                        for (c,conv) in nw_ref[:convdc_ne]
                            if conv["busac_i"] == ac_bus
                             push!(bus_arcs_load_ne[conv["syncarea"]], la)
                             end
                        end
                end
                nw_ref[:bus_arcs_load_ne] = bus_arcs_load_ne

                #old
                bus_arcs_conv_ne = Dict([(reserve["syncarea"], []) for (i,reserve) in nw_ref[:reserves]])
                for (l,i,j) in nw_ref[:arcs_sync_conv_ne]
                    push!(bus_arcs_conv_ne[j], l)
                end
                nw_ref[:bus_arcs_conv_ne] = bus_arcs_conv_ne

                #new
                bus_arcs_incident_ne = Dict([(i, []) for (i,branch) in nw_ref[:branchdc_ne]])
                for (i,incident) in nw_ref[:incident]
                    for b = 1: length(incident)-3
                        if incident["br$(b)"] == 1
                        push!(bus_arcs_incident_ne[b], incident["convdc"])
                        end
                    end
                end
                nw_ref[:bus_arcs_incident_ne] = bus_arcs_incident_ne

                bus_arcs_sync_ne = Dict([(reserve["syncarea"], []) for (i,reserve) in nw_ref[:reserves]])
                for (l,i,j) in nw_ref[:arcs_sync_conv_ne]
                    push!(bus_arcs_sync_ne[j], i)
                end
                nw_ref[:bus_arcs_sync_ne] = bus_arcs_sync_ne

                arcs_reserves_syn_ne = Dict([(i, []) for (i,reserve) in nw_ref[:reserves]])
                for (i,reserve) in nw_ref[:reserves]
                    push!(arcs_reserves_syn_ne[reserve["syncarea"]], i)
                end
                nw_ref[:arcs_reserves_syn_ne] = arcs_reserves_syn_ne

                bus_syc_arc_ne = Dict([(i, []) for (i,bus) in nw_ref[:bus]])
                for (i,conv) in nw_ref[:convdc_ne]
                     push!(bus_syc_arc_ne[conv["busac_i"]],conv["syncarea"])
                end
            end
end

"compute bus pair level structures"
function buspair_parameters_dc_ne(arcs_dcgrid_from, branches, busesdc_ne, busesdc)
    buspair_indexes = collect(Set([(i,j) for (l,i,j) in arcs_dcgrid_from]))

    bp_branch = Dict([(bp, Inf) for bp in buspair_indexes])
    for (l,branch) in branches
        i = branch["fbusdc"]
        j = branch["tbusdc"]
        bp_branch[(i,j)] = min(bp_branch[(i,j)], l)
    end
    buspairs = Dict([((i,j), Dict{String,Any}(
    "branch" => nothing,
    "vm_fr_min" => nothing,
    "vm_fr_max" => nothing,
    "vm_to_min" => nothing,
    "vm_to_max" => nothing
    )) for (i,j) in buspair_indexes])

    for (i,j) in buspair_indexes
        buspairs[(i,j)]["branch"] = bp_branch[(i,j)]
        for (k, bus) in busesdc_ne
            if i == bus["busdc_i"]
                buspairs[(i,j)]["vm_fr_min"] = bus["Vdcmin"]
                buspairs[(i,j)]["vm_fr_max"] = bus["Vdcmax"]
            end
            if j == bus["busdc_i"]
                buspairs[(i,j)]["vm_to_min"] = bus["Vdcmin"]
                buspairs[(i,j)]["vm_to_max"] = bus["Vdcmax"]
            end
        end
        for (k, bus) in busesdc
            if i == bus["busdc_i"]
                buspairs[(i,j)]["vm_fr_min"] = bus["Vdcmin"]
                buspairs[(i,j)]["vm_fr_max"] = bus["Vdcmax"]
            end
            if j == bus["busdc_i"]
                buspairs[(i,j)]["vm_to_min"] = bus["Vdcmin"]
                buspairs[(i,j)]["vm_to_max"] = bus["Vdcmax"]
            end
        end
    end

    return buspairs
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

@testset "graph.connectivity" begin
    emp = plaingraph(5, Tuple{Int,Int}[])
    @test length(connected_components(emp)) == 5
    @test isempty(bridges(emp))
    @test isempty(cutvertices(emp))
    @test length(biconnected_components(emp)) == 0 # no edges
    @test length(two_edge_connected(emp)) == 5

    pg = pathgraph(5)
    @test length(bridges(pg)) == 4
    @test length(connected_components(pg)) == 1
    @test length(biconnected_components(pg)) == 4
    @test length(two_edge_connected(pg)) == 5

    disconn = plaingraph(5, [(1, 2), (2, 3), (4, 5)])
    @test length(connected_components(disconn)) == 2

    cyc = cyclegraph(5)
    @test isempty(bridges(cyc))
    @test length(biconnected_components(cyc)) == 1
    @test length(two_edge_connected(cyc)) == 1

    fused = plaingraph(8, [
        (1, 2), (2, 3), (1, 3), (3, 4), (4, 5),
        (5, 6), (4, 6), (5, 7), (7, 8), (8, 6)
    ])
    bf = bridges(fused)
    @test length(bf) == 1
    @test bf[1] == 4
    @test length(cutvertices(fused)) == 2
    @test length(biconnected_components(fused)) == 3
    @test length(two_edge_connected(fused)) == 2

    branched = plaingraph(9, [
        (1, 2), (2, 3), (3, 4), (4, 5), (4, 6),
        (3, 7), (7, 8), (8, 9), (9, 1)
    ])
    @test length(bridges(branched)) == 3
    @test length(cutvertices(branched)) == 2

    spiro = plaingraph(10, [
        (1, 2), (2, 3), (3, 4), (3, 5), (5, 6),
        (6, 7), (7, 5), (5, 8), (8, 9), (9, 5), (5, 10)
    ])
    @test length(bridges(spiro)) == 5
    @test length(cutvertices(spiro)) == 3
    @test length(biconnected_components(spiro)) == 7
    @test length(two_edge_connected(spiro)) == 6
end # graph.connectivity


text:
using Base.Test

include("dpmm.jl")
include("smc_samplers.jl")

@testset "CRP tests" begin

    @testset "assignment book-keeping" begin
        crp = CRPState()

        incorporate!(crp, 1, next_new_cluster(crp))
        c1 = assignment(crp, 1)
        @test c1 == 1

        incorporate!(crp, 2, next_new_cluster(crp))
        c2 = assignment(crp, 2)
        @test c2 == 2

        incorporate!(crp, 3, c2)
        @test counts(crp, c1) == 1
        @test counts(crp, c2) == 2
        @test next_new_cluster(crp) == c2 + 1
        @test crp.next_cluster == c2 + 2

        unincorporate!(crp, 1)
        @test !has_cluster(crp, c1)
        @test counts(crp, c2) == 2
        @test next_new_cluster(crp) == c1

        unincorporate!(crp, 2)
        @test counts(crp, c2) == 1

        unincorporate!(crp, 3)
        @test !has_cluster(crp, c2)
        @test next_new_cluster(crp) == c2

        # this is never decremented
        @test crp.next_cluster == c2 + 2
    end

    @testset "log_probability" begin

        alpha = 0.25
        crp = CRPState()
        @test log_probability(crp, alpha) == 0.0

        # [1]
        incorporate!(crp, 1, next_new_cluster(crp))
        @test log_probability(crp, alpha) == 0.0

        # join the same table as 1: [1,2]
        incorporate!(crp, 2, assignment(crp, 1))
        actual = log_probability(crp, alpha)
        expected = log(1. / (1. + alpha))
        @test isapprox(actual, expected)

        # join a new table: [1] [2]
        unincorporate!(crp, 2)
        incorporate!(crp, 2, next_new_cluster(crp))
        actual = log_probability(crp, alpha)
        expected = log(alpha / (1. + alpha))
        @test isapprox(actual, expected)

        # join same table as 1: [1, 3] [2]
        incorporate!(crp, 3, assignment(crp, 1))
        actual = log_probability(crp, alpha)
        expected = log((alpha / (1.+alpha)) *
                       (1./(2.+alpha)))
        @test isapprox(actual, expected)

        # join same table as 1 [1, 3, 4] [2]
        incorporate!(crp, 4, assignment(crp, 1))
        actual = log_probability(crp, alpha)
        expected = log((alpha / (1.+alpha)) *
                       (1./(2.+alpha)) *
                       (2./(3.+alpha)))
        @test isapprox(actual, expected)
    end

end

@testset "NIGN" begin

    # arbitrary values
    m = 1.23
    r = 2.41
    nu = 5.44
    s = 4.2
    alpha = nu/2.
    beta = s/2.
    params = NIGNParams(m, r, nu, s)

    nign = NIGN()
    actual = log_probability_density(nign, params)
    @test abs(actual) < 1e-15

    # https://www.cs.ubc.ca/~murphyk/Papers/bayesGauss.pdf
    x1 = -4.1
    x2 = 2.22
    incorporate!(nign, x1)
    incorporate!(nign, x2)
    (mu_actual, r_actual, s_actual, nu_actual) = posterior_params(nign, params)
    # m, r, s, nu
    avg = mean([x1, x2])
    sum_sq = sum(([x1, x2] - avg) .* ([x1,x2] - avg))
    n = 2
    mu_expect = (r * m + n * avg) / (r + n)
    r_expect = r + n
    alpha_expect = alpha + n/2.
    beta_expect = beta + 0.5 * sum_sq  + r*n*((avg-m)*(avg-m))/(2*(r+n))
    nu_expect = alpha_expect*2.
    s_expect = beta_expect*2.
    @test isapprox(mu_actual, mu_expect)
    @test isapprox(r_actual, r_expect)
    @test isapprox(s_actual, s_expect)
    @test isapprox(nu_actual, nu_expect)

    # should go back to empty again
    unincorporate!(nign, x1)
    unincorporate!(nign, x2)
    actual = log_probability_density(nign, params)
    @test abs(actual) < 1e-15

end

@testset "DPMM SMC" begin

    @testset "extreme easy two-cluster dataset" begin
        srand(1)

        # generate dataset
        num_data = 50
        data = Array{Float64,1}(num_data)
        cluster_1_members = []
        cluster_2_members = []
        for i=1:num_data
            if rand() < 0.5
                mu = -1.0
                push!(cluster_1_members, i)
            else
                mu = 1.0
                push!(cluster_2_members, i)
            end
            std = 0.01
            data[i] = rand(Normal(mu, std))
        end

        # set custom values of hypers to match extreme dataset
        m_prior = Normal(0, 1) # mean cluster center
        r_prior = Gamma(1, 1) # r * data-precision is the precision of cluster centers
        nu_prior = Gamma(1, 1)# expected data-precision of cluster is nu/s
        s_prior = Gamma(1, 1)
        nign_params_prior = NIGNParamsPrior(m_prior, r_prior, nu_prior, s_prior)
        alpha_prior = Gamma(1, 1)
        hypers = Hyperparameters(nign_params_prior, alpha_prior)

        gibbs_sweeps_per_datum = 1
        num_particles = 1
        scheme = make_optimal_proposal_smc_scheme(data, hypers,
                               gibbs_sweeps_per_datum, num_particles)
        state, lml_estimate = smc(scheme)

        # check that there are two clusters, and that clustering
        # is correct
        @test length(state.components) == 2
        cluster_1 = assignment(state.crp, cluster_1_members[1])
        cluster_2 = assignment(state.crp, cluster_2_members[2])
        @test cluster_1 != cluster_2
        for i in cluster_1_members[2:end]
            @test assignment(state.crp, i) == cluster_1
        end
        for i in cluster_2_members[2:end]
            @test assignment(state.crp, i) == cluster_2
        end

    end

    @testset "prior proposal weight" begin
        alpha = 0.01
        m = 0.9
        r = 1.1
        nu = 1.2
        s = 1.3
        nign_params = NIGNParams(m, r, nu, s)
        alpha_prior = Gamma(1, 1)
        nign_params_prior = NIGNParamsPrior(Normal(0, 1), Gamma(1, 1), Gamma(1, 1), Gamma(1, 1))
        hypers = Hyperparameters(nign_params_prior, alpha_prior)
	    state = DPMMState(alpha, CRPState(), OrderedDict{Int,NIGN}(),
                      nign_params, hypers)

        srand(1)

        x1 = 4.2
        ll1 = add_datum_prior!(state, 1, x1)
        ll1_rev = remove_datum_prior!(state, 1, x1)
        @test isapprox(ll1, ll1_rev)

        x2 = -3.1
        ll1_new = add_datum_prior!(state, 1, x1)
        @test isapprox(ll1, ll1_new)
        ll2 = add_datum_prior!(state, 2, x2)
        # NOTE: we seeded so this should remain true, if it doesn't change the
        # seed or make alpha even smaller
        @test assignment(state.crp, 1) == assignment(state.crp, 2)
        component = NIGN()
        incorporate!(component, x1)
        incorporate!(component, x2)
        total_logp = log_probability_density(component, nign_params)
        @test isapprox(total_logp, ll1_new + ll2)

        ll1_rev = remove_datum_prior!(state, 1, x1)
        ll2_rev = remove_datum_prior!(state, 2, x2)
        @test isapprox(total_logp, ll1_rev + ll2_rev)
    end

    @testset "optimal proposal weight" begin
        alpha = 2.45
        m = 0.9
        r = 1.1
        nu = 1.2
        s = 1.3
        nign_params = NIGNParams(m, r, nu, s)
        alpha_prior = Gamma(1, 1)
        nign_params_prior = NIGNParamsPrior(Normal(0, 1), Gamma(1, 1), Gamma(1, 1), Gamma(1, 1))
        hypers = Hyperparameters(nign_params_prior, alpha_prior)
	    state = DPMMState(alpha, CRPState(), OrderedDict{Int,NIGN}(),
                      nign_params, hypers)

        srand(1)

        x1 = 4.2
        ll1 = add_datum_conditional!(state, 1, x1)
        ll1_rev = remove_datum_conditional!(state, 1, x1)
        @test isapprox(ll1, ll1_rev)
        component = NIGN()
        incorporate!(component, x1)
        @test isapprox(ll1, log_probability_density(component, nign_params))

        ll1_new = add_datum_conditional!(state, 1, x1)
        x2 = 1.5
        @test isapprox(ll1, ll1_new)
        ll2 = add_datum_conditional!(state, 2, x2)
        ll2_rev = remove_datum_conditional!(state, 2, x2)
        @test isapprox(ll2, ll2_rev)

        # explicitly sum over the two possible assignments
        # p(x_2 = alone | x_1) p(y_2 | x_1, x_2 = alone, y_1)
        component_alone = NIGN()
        incorporate!(component_alone, x2)
        ll2_alone = log_probability_density(component_alone, nign_params)
        ll2_alone += log(alpha / (1. + alpha))

        # p(x_2 = together | x_1) p(y_2 | x_1, x_2 = together, y_1)
        component_together = NIGN()
        incorporate!(component_together, x1)
        component_old_logpdf = log_probability_density(component_together, nign_params)
        incorporate!(component_together, x2)
        ll2_together = log_probability_density(component_together, nign_params) - component_old_logpdf
        ll2_together += log(1. / (1. + alpha))

        # p(y_2 | x_1, y_1) = sum of two terms above
        ll2_expected = logsumexp([ll2_alone, ll2_together])
        @test isapprox(ll2, ll2_expected)
    end

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

let nextidx = 0
    global nextproc
    function nextproc()
        p = -1
        if p == -1
            p = workers()[(nextidx % nworkers()) + 1]
            nextidx += 1
        end
        p
    end
end

spawnat(p, thunk) = sync_add(remotecall(thunk, p))

spawn_somewhere(thunk) = spawnat(nextproc(),thunk)

macro spawn(expr)
    thunk = esc(:(()->($expr)))
    :(spawn_somewhere($thunk))
end

macro spawnat(p, expr)
    thunk = esc(:(()->($expr)))
    :(spawnat($(esc(p)), $thunk))
end

"""
    @fetch

Equivalent to `fetch(@spawn expr)`.
See [`fetch`](@ref) and [`@spawn`](@ref).
"""
macro fetch(expr)
    thunk = esc(:(()->($expr)))
    :(remotecall_fetch($thunk, nextproc()))
end

"""
    @fetchfrom

Equivalent to `fetch(@spawnat p expr)`.
See [`fetch`](@ref) and [`@spawnat`](@ref).
"""
macro fetchfrom(p, expr)
    thunk = esc(:(()->($expr)))
    :(remotecall_fetch($thunk, $(esc(p))))
end

# extract a list of modules to import from an expression
extract_imports!(imports, x) = imports
function extract_imports!(imports, ex::Expr)
    if Meta.isexpr(ex, (:import, :using))
        return push!(imports, ex.args[1])
    elseif Meta.isexpr(ex, :let)
        return extract_imports!(imports, ex.args[1])
    elseif Meta.isexpr(ex, (:toplevel, :block))
        for i in eachindex(ex.args)
            extract_imports!(imports, ex.args[i])
        end
    end
    return imports
end
extract_imports(x) = extract_imports!(Symbol[], x)

"""
    @everywhere expr

Execute an expression under `Main` everywhere. Equivalent to calling
`eval(Main, expr)` on all processes. Errors on any of the processes are collected into a
`CompositeException` and thrown. For example :

    @everywhere bar=1

will define `Main.bar` on all processes.

Unlike [`@spawn`](@ref) and [`@spawnat`](@ref),
`@everywhere` does not capture any local variables. Prefixing
`@everywhere` with [`@eval`](@ref) allows us to broadcast
local variables using interpolation :

    foo = 1
    @eval @everywhere bar=\$foo

The expression is evaluated under `Main` irrespective of where `@everywhere` is called from.
For example :

    module FooBar
        foo() = @everywhere bar()=myid()
    end
    FooBar.foo()

will result in `Main.bar` being defined on all processes and not `FooBar.bar`.
"""
macro everywhere(ex)
    imps = [Expr(:import, m) for m in extract_imports(ex)]
    quote
        $(isempty(imps) ? nothing : Expr(:toplevel, imps...))
        sync_begin()
        for pid in workers()
            async_run_thunk(()->remotecall_fetch(eval_ew_expr, pid, $(Expr(:quote,ex))))
            yield() # ensure that the remotecall_fetch has been started
        end

        # execute locally last as we do not want local execution to block serialization
        # of the request to remote nodes.
        if nprocs() > 1
            async_run_thunk(()->eval_ew_expr($(Expr(:quote,ex))))
        end

        sync_end()
    end
end

eval_ew_expr(ex) = (eval(Main, ex); nothing)

# Statically split range [1,N] into equal sized chunks for np processors
function splitrange(N::Int, np::Int)
    each = div(N,np)
    extras = rem(N,np)
    nchunks = each > 0 ? np : extras
    chunks = Vector{UnitRange{Int}}(nchunks)
    lo = 1
    for i in 1:nchunks
        hi = lo + each - 1
        if extras > 0
            hi += 1
            extras -= 1
        end
        chunks[i] = lo:hi
        lo = hi+1
    end
    return chunks
end

function preduce(reducer, f, R)
    N = length(R)
    chunks = splitrange(N, nworkers())
    all_w = workers()[1:length(chunks)]

    w_exec = Task[]
    for (idx,pid) in enumerate(all_w)
        t = Task(()->remotecall_fetch(f, pid, reducer, R, first(chunks[idx]), last(chunks[idx])))
        schedule(t)
        push!(w_exec, t)
    end
    reduce(reducer, [wait(t) for t in w_exec])
end

function pfor(f, R)
    [@spawn f(R, first(c), last(c)) for c in splitrange(length(R), nworkers())]
end

function make_preduce_body(var, body)
    quote
        function (reducer, R, lo::Int, hi::Int)
            $(esc(var)) = R[lo]
            ac = $(esc(body))
            if lo != hi
                for $(esc(var)) in R[(lo+1):hi]
                    ac = reducer(ac, $(esc(body)))
                end
            end
            ac
        end
    end
end

function make_pfor_body(var, body)
    quote
        function (R, lo::Int, hi::Int)
            for $(esc(var)) in R[lo:hi]
                $(esc(body))
            end
        end
    end
end

"""
    @parallel

A parallel for loop of the form :

    @parallel [reducer] for var = range
        body
    end

The specified range is partitioned and locally executed across all workers. In case an
optional reducer function is specified, `@parallel` performs local reductions on each worker
with a final reduction on the calling process.

Note that without a reducer function, `@parallel` executes asynchronously, i.e. it spawns
independent tasks on all available workers and returns immediately without waiting for
completion. To wait for completion, prefix the call with [`@sync`](@ref), like :

    @sync @parallel for var = range
        body
    end
"""
macro parallel(args...)
    na = length(args)
    if na==1
        loop = args[1]
    elseif na==2
        reducer = args[1]
        loop = args[2]
    else
        throw(ArgumentError("wrong number of arguments to @parallel"))
    end
    if !isa(loop,Expr) || loop.head !== :for
        error("malformed @parallel loop")
    end
    var = loop.args[1].args[1]
    r = loop.args[1].args[2]
    body = loop.args[2]
    if na==1
        thecall = :(pfor($(make_pfor_body(var, body)), $(esc(r))))
    else
        thecall = :(preduce($(esc(reducer)), $(make_preduce_body(var, body)), $(esc(r))))
    end
    thecall
end

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
function Preprocess(VideoString::String)
    AV = VideoIO.open(VideoString)
    videoObj = VideoIO.openvideo(AV)
    raw = VideoIO.read(videoObj)
    mytype = typeof(Gray.(raw))
    img = convert(Array{Float64},Gray.(raw)).*255

    ImgSet = [img]
    counter = 0
    while !eof(videoObj)
        read!(videoObj, raw)
        img = convert(Array{Float64},Gray.(raw)).*255
        push!(ImgSet,img)
        counter+=1
    end
    close(videoObj)

    return ImgSet
end

export Preprocess


text:
using ClusterManagers, DSGE, ModelConstructors, FileIO

fn = dirname(@__FILE__)

# This script auto-generates a packet of
# plots and figures to analyze
# estimation, forecast, and impulse response results.
# This script also serves as an example of how we recommend
# structuring "master" files for running forecasts, i.e.
# this file has a recommended workflow for the
# entirety of a forecast.

# A note also on the options. Generally, the only memory-intensive
# operation is the computation of the shock decompositions,
# which involves handling many large matrices. If memory
# needs to be pre-allocated, then 1GB is often sufficient
# if you do not run the shock decompositions portion of the code.

# What do you want to do?
run_full_forecast = true  # Run full distribution forecast
do_histforecast   = true  # Compute history and forecast of observables and pseudo-observables
do_shockdecs      = true  # Compute shock decompositions
do_irfs           = true  # Compute impulse response functions
make_packet       = true  # Create packet w/info on estimation, forecasts, and IRFs
make_plots        = true  # Make plots
add_workers       = false # Run in parallel
n_workers         = 10

# Initialize model objects and desired settings
m = Model1002("ss10")
m <= Setting(:sampling_method, :MH)
usual_settings!(m, "181115"; cdvt = "181115", fcast_date = quartertodate("2018-Q4"))
m <= Setting(:use_population_forecast, false) # We do not make this data available
m <= Setting(:saveroot, "$(fn)/../../save/") # set save root
m <= Setting(:dataroot, "$(fn)/../../save/input_data/") # set data root
m <= Setting(:date_forecast_end, quartertodate("2030-Q1")) # When to stop forecast
forecast_string = "" # to identify different forecasts from each other
m <= Setting(:forecast_block_size, 40) # this depends on the number of samples from estimations
                                      # and thinning during forecast (determined by
                                      # the setting :forecast_jstep)
if get_setting(m, :sampling_method) == :SMC
    # If SMC is the estimation method, then we need t adjust some default settings
    m <= Setting(:forecast_jstep, 1)
end

# Load data to create data set. There will be empty columns
# since not all the data used is available. This requires
# setting the check_empty_columns flag to false.
df = load_data(m; check_empty_columns = false)

# Use overrides to give the correct file path for saved estimation output
overrides = forecast_input_file_overrides(m)
overrides[:full] = "$(fn)/../../test/reference/mhsave_vint=181115.h5"

# Full-distribution forecast
if run_full_forecast

    output_vars = Vector{Symbol}(undef,0)
    if do_histforecast
        # Write data to create historical and forecast output
        output_vars = vcat(output_vars, [:histpseudo, :histobs, :histstdshocks,
                                         :hist4qpseudo, :hist4qobs, :histutpseudo,
                                         :forecastpseudo, :forecastobs, :forecastutpseudo,
                                         :forecast4qpseudo, :forecast4qobs, :forecaststdshocks])
    end

    if do_shockdecs
        # Shock decompositions of forecasts
        output_vars = vcat(output_vars, [:dettrendobs, :dettrendpseudo, :trendobs,
                                         :trendpseudo, :shockdecpseudo, :shockdecobs])
    end

    if do_irfs
        # Impulse response to all shocks, including for endogenous states
        output_vars = vcat(output_vars, [:irfstates, :irfobs, :irfpseudo])
    end

    if add_workers
        my_procs = addprocs(n_workers)
        @everywhere using DSGE, OrderedCollections
    end

    usual_forecast(m, :full, :none, output_vars,
                   est_override = overrides[:full],
                   forecast_string = forecast_string,
                   density_bands = [.5, .6, .68, .7, .8, .9],
                   check_empty_columns = false)
end

# Make plots and packet
output_vars = [:forecastobs, :forecastpseudo]
if do_shockdecs
    output_vars = vcat(output_vars, [:shockdecobs, :shockdecpseudo])
end

if do_irfs
    # we do not plot irfs for endogenous states b/c too many!
    output_vars = vcat(output_vars, [:irfobs, :irfpseudo])
    sections = [:estimation, :forecast, :irf]
else
    sections = [:estimation, :forecast]
end

if make_plots
    plot_standard_packet(m, :full, :none, output_vars,
                         forecast_string = forecast_string,
                         sections = sections)
end

if make_packet
    # Choose forecast centered or the standard packet.
    # Note that these functions write to the same tex file,
    # so you can use only one for a given data vintage!
    # write_forecast_centric_packet(m, :full, :none, output_vars,
    #                               sections = sections, forecast_string = forecast_string)
    write_standard_packet(m, :full, :none, output_vars,
                          sections = sections, forecast_string = forecast_string)
    moment_tables(m)
end

if add_workers
    rmprocs(my_procs)
end

nothing


text:
import Xfoil
import LiftProps
using Plots
pyplot()

airfoil_file = "naca2412.dat"

open(airfoil_file,"r") do f
  global x = Float64[]
  global y = Float64[]
  for line in eachline(f)
    x = append!(x,Meta.parse(split(chomp(line))[1]))
    y = append!(y,Meta.parse(split(chomp(line))[2]))
  end
end

aoa = range(-15,stop=20,length=45)

cl,cd,cdp,cm,converged = Xfoil.xfoilsweep(x,y,aoa,400000.0,printdata=true,clminstop=true,clmaxstop=true)

aoa = aoa[findall(converged)]
cl = cl[findall(converged)]
cd = cd[findall(converged)]
cdp = cdp[findall(converged)]
cm = cm[findall(converged)]

liftslope,zeroliftangle,aoafit,clfit = LiftProps.fitliftslope(aoa,cl)
_,clmax1 = LiftProps.findclmax(aoa,cl)
_,clmax2 = LiftProps.findclmax(aoa,cl,2)
_,clmax3 = LiftProps.findclmax(aoa,cl,1.0)
_,clmaxlinear1 = LiftProps.findclmaxlinear(aoa,cl,liftslope,zeroliftangle,interpolate=false)
_,clmaxlinear2 = LiftProps.findclmaxlinear(aoa,cl,liftslope,zeroliftangle,interpolate=true)
_,clmin1 = LiftProps.findclmin(aoa,cl)
_,clmin2 = LiftProps.findclmin(aoa,cl,2)
_,clmin3 = LiftProps.findclmin(aoa,cl,1.0)

# LIFT PLOT
plot(overwrite_figure = false)
myfont = font(12)
scatter!(aoa*180/pi,cl,label = "",markersize=5)
plot!(ylims = [min(minimum(cl),-0.1), maximum(cl)*1.1])
plot!(xlabel = "\$\\alpha\$ (degrees)",ylabel = "\$c_{l}\$")
plot!(legendfont = myfont,guidefont = myfont,tickfont = myfont)
# Plot XFOIL Data
scatter!(aoa*180/pi,cl,label = "XFOIL Data",markersize=5,markercolor=:red)
# Plot Linear Lift Slope Line
indclmin = argmin(abs.(cl.-minimum(cl)))
liftslopelinex = range(min(zeroliftangle-2.5*pi/180, aoa[indclmin]),stop=maximum(aoa),length=1000)
liftslopeliney = liftslope*(liftslopelinex.-zeroliftangle)
minAoAidx = argmin(abs.(minimum(cl) .- liftslopeliney))
maxAoAidx = argmin(abs.(maximum(cl) .- liftslopeliney))
minAoA = liftslopelinex[minAoAidx]
maxAoA = liftslopelinex[maxAoAidx]
liftslopelinex = range(min(zeroliftangle.-2.5*pi/180,minAoA), stop=maxAoA, length=1000)
liftslopeliney = liftslope*(liftslopelinex.-zeroliftangle)
plot!(liftslopelinex*180/pi,liftslopeliney, linewidth=2,
    label="Linear Lift Slope Line", linecolor=:blue)
# Plot Linear Lift Slope Line Data
scatter!(aoafit*180/pi,clfit,label = "Least Squares Data",
    markersize = 8,markercolor = :red,markershape = :x)
# Plot Zero Lift Angle
plot!([zeroliftangle*180/pi],seriestype = [:vline],label = "\$\\alpha_{0}\$",
    linestyle = :dashdot,linewidth = 1.5,linecolor = :brown)
# Plot clmax
plot!([clmax1],seriestype = [:hline],label = "\$c_{l,max,1}\$",
    linestyle = :dash,linewidth = 1,linecolor = :green)
plot!([clmax2],seriestype = [:hline],label = "\$c_{l,max,2}\$",
    linestyle = :dashdot,linewidth = 2,linecolor = :green)
plot!([clmax3],seriestype = [:hline],label = "\$c_{l,max,3}\$",
    linestyle = :dot,linewidth = 3,linecolor = :green)
plot!(show = true)
# Plot clmaxlinear
plot!([clmaxlinear1],seriestype = [:hline],label = "\$c_{l,max,linear,1}\$",
    linestyle = :dash,linewidth = 1,linecolor = :orange)
plot!([clmaxlinear2],seriestype = [:hline],label = "\$c_{l,max,linear,2}\$",
    linestyle = :dashdot,linewidth = 2,linecolor = :orange)
# Plot clmin
plot!([clmin1],seriestype = [:hline],label = "\$c_{l,min,1}\$",
    linestyle = :dash,linewidth = 1,linecolor = :purple)
plot!([clmin2],seriestype = [:hline],label = "\$c_{l,min,2}\$",
    linestyle = :dashdot,linewidth = 2,linecolor = :purple)
plot!([clmin3],seriestype = [:hline],label = "\$c_{l,min,3}\$",
    linestyle = :dot,linewidth = 3,linecolor = :purple)
plot!(show = true)


text:
#############################################################################
#############################################################################
#
# This file defines the polynomial type with several operations
#
#############################################################################
#############################################################################

####################################
# Polynomial type and construction #
####################################

"""
A Polynomial type - designed to be for polynomials with integer coefficients.
"""
struct Polynomial
    terms::Array{Any,1}
        #The terms in the stack need to satisfy:
            # Will never have terms with 0 coefficient
            # Will never have two terms with same coefficient
        #An empty terms stack means that the polynomial is zero
    Polynomial() = new(Array{Any,1}[])
    #Inner constructor
    Polynomial(h::Array{Any,1}) = new(h)
end

struct PolynomialModP
    terms::Polynomial
    prime::Int

    #Constructor for the case where only a prime is given, initialize an empty polynomial and enforce mod(prime) conditions
    PolynomialModP(prime::Int) = new(Polynomial(), prime)
    #Inner constructor
    PolynomialModP(poly::Polynomial, prime::Int) = new(poly, prime)
end

"""
Construct a polynomial with a single term.
"""
function Polynomial(t::Term)
    terms = Array{Any,1}([])
    t.coeff != 0 && push!(terms, t)
    return Polynomial(terms)
end

"""
Construct a polynomial with a vector of terms.
"""
function Polynomial(tv::Vector{Term})
    terms = Array{Any,1}([])
    for t in tv
        t.coeff != 0 && push!(terms,t)
    end
    return Polynomial(terms)
end

"""
Construct a polynomial of the form x^p-x.
"""
cyclotonic_polynomial(p::Int) = Polynomial([Term(1,p), Term(-1,0)])

"""
Construct a polynomial of the form x-n.
"""
linear_monic_polynomial(n::Int) = Polynomial([Term(1,1), Term(-n,0)])

"""
Construct a polynomial of the form x.
"""
x_poly() = Polynomial(Term(1,1))

"""
Creates the zero polynomial.
"""
zero(::Type{Polynomial})::Polynomial = Polynomial()

"""
Creates the unit polynomial.
"""
one(::Type{Polynomial})::Polynomial = Polynomial(one(Term))
one(p::Polynomial) = one(typeof(p))

"""
Generates a random polynomial.
"""
function rand(::Type{Polynomial} ;
                degree::Int = -1,
                terms::Int = -1,
                max_coeff::Int = 100,
                mean_degree::Float64 = 5.0,
                prob_term::Float64  = 0.7,
                monic = false,
                condition = (p)->true)

    while true
        _degree = degree == -1 ? rand(Poisson(mean_degree)) : degree
        _terms = terms == -1 ? rand(Binomial(_degree,prob_term)) : terms
        degrees = vcat(sort(sample(0:_degree-1,_terms,replace = false)),_degree)
        coeffs = rand(1:max_coeff,_terms+1)
        monic && (coeffs[end] = 1)
        p = Polynomial( [Term(coeffs[i],degrees[i]) for i in 1:length(degrees)] )
        condition(p) && return p
    end
end

###########
# Display #
###########

"""
Show a polynomial.
"""
function show(io::IO, p::Polynomial)
    p = deepcopy(p)
    if iszero(p)
        print(io,"0")
    else
        n = length(p.terms)
        for (i,t) in enumerate(p.terms)
            print(io, t, i != n ? " + " : "")
        end
    end
end

"""
Show a polynomial Modulo P.
"""
function show(io::IO, p::PolynomialModP)
    input_poly = deepcopy(p.terms)
    if iszero(input_poly)
        print(io."0")
    else
        n = length(input_poly.terms)
        for (i,t) in enumerate(input_poly.terms)
            print(io, t, i != n ? " + " : "")
        end
        print(" mod ", p.prime)
    end
end

##############################################
# Iteration over the terms of the polynomial #
##############################################

"""
Allows to do iteration over the terms of the polynomial. The iteration is in an arbitrary order.
"""
iterate(p::Polynomial, state=1) = iterate(p.terms, state)

##############################
# Queries about a polynomial #
##############################

"""
The number of terms of the polynomial.
"""
length(p::Polynomial) = length(p.terms)

"""
The leading term of the polynomial.
"""
#Assuming every array of polynomials given to the struct is in ascending order of the term degrees.
leading(p::Polynomial)::Term = isempty(p.terms) ? zero(Term) : last(p.terms)

"""
The trailing term of the polynomial.
"""
#Assuming every array of polynomials given to the struct is in ascending order of the term degrees.
trailing(p::Polynomial)::Term = isempty(p.terms) ? zero(Term) : first(p.terms)

"""
Returns the coefficients of the polynomial.
"""
coeffs(p::Polynomial)::Vector{Int} = [t.coeff for t in p]

"""
The degree of the polynomial.
"""
degree(p::Polynomial)::Int = leading(p).degree

"""
The content of the polynomial is the GCD of its coefficients.
"""
content(p::Polynomial)::Int = euclid_alg(coeffs(p))

"""
Evaluate the polynomial at a point `x`.
"""
evaluate(f::Polynomial, x::T) where T <: Number = sum(evaluate(t,x) for t in f)

################################
# Pushing and popping of terms #
################################

"""
Push a new term into the polynomial.
"""
#Note that ideally this would throw and error if pushing another term of degree that is already in the polynomial
function push!(p::Polynomial, t::Term)
    iszero(t) && return #don't push a zero
    push!(p.terms,t)
end

"""
Pop the trailing term out of the polynomial.
"""
pop!(p::Polynomial)::Term = pop!(p.terms)

"""
Pop the leading term out of the polynomial.
"""
popfirst!(p::Polynomial)::Term = popfirst!(p.terms)

"""
Check if the polynomial is zero.
"""
iszero(p::Polynomial)::Bool = isempty(p.terms)

#################################################################
# Transformation of the polynomial to create another polynomial #
#################################################################

"""
The negative of a polynomial.
"""
-(p::Polynomial) = Polynomial(map((pt)->-pt, p.terms))

"""
Create a new polynomial which is the derivative of the polynomial.
"""
function derivative(p::Polynomial)::Polynomial
    der_p = Polynomial()
    for term in p
        der_term = derivative(term)
        !iszero(der_term) && push!(der_p,der_term)
    end
    return der_p
end

"""
The prim part (multiply a polynomial by the inverse of its content).
"""
prim_part(p::Polynomial) = p ÷ content(p)

"""
A square free polynomial.
"""
square_free(p::Polynomial, prime::Int)::Polynomial = (PolynomialModP(p, prime) ÷ PolynomialModP(gcd(PolynomialModP(p, prime),PolynomialModP(derivative(p), prime)), prime))

#################################
# Queries about two polynomials #
#################################

"""
Check if two polynomials are the same
"""
==(p1::Polynomial, p2::Polynomial)::Bool = p1.terms == p2.terms

"""
Check if a polynomial is equal to 0.
"""
#Note that in principle there is a problem here. E.g The polynomial 3 will return true to equalling the integer 2.
==(p::Polynomial, n::T) where T <: Real = iszero(p) == iszero(n)

##################################################################
# Operations with two objects where at least one is a polynomial #
##################################################################

"""
Subtraction of two polynomials.
"""
-(p1::Polynomial, p2::Polynomial)::Polynomial = p1 + (-p2)

"""
Multiplication of polynomial and term.
"""
*(t::Term,p1::Polynomial)::Polynomial = iszero(t) ? Polynomial() : Polynomial(map((pt)->t*pt, p1.terms))
*(p1::Polynomial, t::Term)::Polynomial = t*p1

"""
Multiplication of polynomial and an integer.
"""
*(n::Int,p::Polynomial)::Polynomial = p*Term(n,0)
*(p::Polynomial,n::Int)::Polynomial = n*p

"""
Integer division of a polynomial by an integer.

Warning this may not make sense if n does not divide all the coefficients of p.
"""
÷(p::Polynomial,n::Int) = (prime)->Polynomial(map((pt)->((pt ÷ n)(prime)), p.terms))

"""
Integer division of a polynomial mod(P) by an integer.

Warning this may not make sense if n does not divide all the coefficients of p.
"""
÷(p::PolynomialModP,n::Int) = PolynomialModP(Polynomial(map((pt)->((pt ÷ n)(p.prime)), p.terms)), p.prime)

"""
Take the smod of a polynomial with an integer.
"""
function mod(f::Polynomial, p::Int)::Polynomial
    p_out = Polynomial()
    for t in f
        push!(p_out, mod(t, p)) #if coeff reduced to zero, push! will handle it
    end
    return p_out
end

"""
Power of a polynomial mod prime.
"""
# function pow_mod(p::Polynomial, n::Int, prime::Int)
#     n < 0 && error("No negative power")
#     out = one(p)
#     for _ in 1:n
#         out *= p
#         out = mod(out, prime)
#     end
#     return out
# end

"""
Power of a polynomial mod prime.
"""
function ^(p::PolynomialModP, n::Int)
    n < 0 && error("No negative power")
    out = one(p.terms)
    for _ in 1:n
        out *= p.terms
        out = mod(out, p.prime)
    end
    return out
end


text:
using Test
using YaoCompiler
using YaoLocations
using YaoArrayRegister
using YaoCompiler.Intrinsics

@device function test_intrinsic(theta, phi)
    # syntax sugar
    1 => X
    @gate 2 => Z
    @ctrl 1 4 => Rx(theta)
    @ctrl 2 4 => Ry(phi)
    return
end

buf = IOBuffer()
show(buf, MIME"text/plain"(), test_intrinsic)
@test occursin("test_intrinsic (generic routine with 1 methods)", String(take!(buf)))

show(buf, MIME"text/plain"(), X)
@test occursin("X (intrinsic operation)", String(take!(buf)))

show(buf, MIME"text/plain"(), Rx(2.0))
@test occursin("Rx(2.0) (intrinsic operation)", String(take!(buf)))

show(buf, MIME"text/plain"(), typeof(test_intrinsic))
@test occursin("typeof(test_intrinsic)", String(take!(buf)))

show(buf, AnyReg())
@test occursin("AnyReg()", String(take!(buf)))


text:
# Model
include_model("cybertruck")

# Horizon
T = 26

# Time step
tf = 1.0
h = tf / (T-1)

_uu = Inf * ones(model.m)
_uu[model.idx_u] = [Inf; 1.0]
_ul = zeros(model.m)
_ul[model.idx_u] .= [0.0; -1.0]

ul, uu = control_bounds(model, T, _ul, _uu)

# Initial and final states
q1 = [0.0, 1.0, 0.0, -1.0 * pi / 2.0]
qT = [3.0, 0.0, 0.0, pi / 2.0]

x1 = [q1; q1]
xT = [qT; qT]

xl, xu = state_bounds(model, T, x1 = x1)

# Objective
Qq = Diagonal(ones(model.nq))
Q = cat(0.5 * Qq, 0.5 * Qq, dims = (1, 2))
QT = cat(10000.0 * Qq, 10000.0 * Diagonal(ones(model.nq)), dims = (1, 2))
R = Diagonal([1.0e-3 * ones(model.nu)..., zeros(model.m - model.nu)...])

obj_tracking = quadratic_tracking_objective(
    [t < T ? Q : QT for t = 1:T],
    [R for t = 1:T-1],
    [xT for t = 1:T],
    [zeros(model.m) for t = 1:T])

obj_penalty = PenaltyObjective(1.0e5, model.m)
obj = MultiObjective([obj_tracking, obj_penalty])

# Constraints
include_constraints(["obstacles", "control_complementarity", "contact"])
p_car1 = [3.0, 0.65]
p_car2 = [3.0, -0.65]

function obstacles!(c, x)
    c[1] = circle_obs(x[1], x[2], p_car1[1], p_car1[2], 0.5)
    c[2] = circle_obs(x[1], x[2], p_car2[1], p_car2[2], 0.5)
    nothing
end

n_obs_stage = 2
n_obs_con = n_obs_stage * T
con_obstacles = ObstacleConstraints(n_obs_con, (1:n_obs_con), n_obs_stage)

n_cc_stage = 2
n_cc_con = n_cc_stage * (T - 1)
con_ctrl_comp = ControlComplementarity(n_cc_con, (1:n_cc_con), n_cc_stage)

con_contact = contact_constraints(model, T)

con = multiple_constraints([con_contact, con_ctrl_comp, con_obstacles])

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# Trajectory initialization
q_ref = linear_interpolation(q1, qT, T)
x_ref = configuration_to_state(q_ref)
x0 = deepcopy(x_ref)
u0 = [1.0e-4 * rand(model.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
include_snopt()
@time z̄ , info = solve(prob, copy(z0),
    nlp = :SNOPT7,
    tol = 1.0e-5, c_tol = 1.0e-5)

check_slack(z̄, prob)

x̄, ū = unpack(z̄, prob)

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
open(vis)
visualize!(vis, model,
    state_to_configuration([[x̄[1] for i = 1:10]...,x̄..., [x̄[end] for i = 1:10]...]),
    Δt = h)