text:
module Exterior

using NLsolve
#using DifferentialEquations
using OrdinaryDiffEq # using this module to reduce slow precompile of full DifferentialEquations package

using Compat
using Compat: reverse, repeat, findmin
using Compat.LinearAlgebra
using Compat.Statistics

using ..MapTypes
using ..Polygons

include("Integration.jl")
using .Integration

include("Reindex.jl")
using .Reindex

export PowerSeries,PowerSeriesDerivatives,PowerMap,ExteriorMap,
        KarmanTrefftzMap,JoukowskiMap,
        summary,parameters,coefficients,
        moments,area,centroid,Jmoment,addedmass

include("exterior/powerseries.jl")
include("exterior/polygon.jl")

struct PowerMap <: ConformalMap
    "Power series"
    ps::PowerSeries

    "Derivatives of power series"
    dps::PowerSeriesDerivatives

    "number of plotting control points"
    N::Int

    "control point coordinates in circle space"
    ζ::Vector{ComplexF64}

    "control point coordinates in body-fixed space"
    z::Vector{ComplexF64}

    "map Jacobian in body-fixed coordinates"
    dzdζ::Vector{ComplexF64}

    "Area enclosed by the mapped shape"
    area      :: Float64

    "Centroid of the mapped shape"
    Zc        :: ComplexF64

    "2nd area moment of the mapped shape"
    J         :: Float64

    "Added mass tensor"
    Ma        :: Array{Float64,2}

end

circle(N) = [exp(im*2π*(i-1)/N) for i in 1:N]

@doc raw"""
    PowerMap(c::Vector{ComplexF64}[;N = 200]) <: ConformalMap

Create a power series map from the exterior of the unit
circle to the exterior of a shape defined by the power series coefficients `c`.

The form of the mapping is

```math
z(\zeta) = c_{1}\zeta + c_{0} + \sum_{j=1}^{N_{c}} \frac{c_{-j}}{\zeta^{j}}
```

The entries in `c` correspond as follows: `c[1]` ``\rightarrow c_{1}``,
`c[2]` ``\rightarrow c_{0}``, `c[3]` ``\rightarrow c_{-1}``, etc.

The resulting map `m` can be evaluated at a single or a vector of points `ζ`
with `m(ζ)`.

# Example

```jldoctest
julia> c = ComplexF64[1,0,1/4];

julia> m = PowerMap(c)
Power series map

julia> ζ = [1.0+3.0im,-2.0-2.0im,0.0+1.1im];

julia> m(ζ)
3-element Array{Complex{Float64},1}:
   1.025+2.925im
 -2.0625-1.9375im
     0.0+0.872727im
```
"""
function PowerMap(ccoeff::Vector{ComplexF64}; N::Int = 200)
  # Must have at least two entries. If no entries, make leading entry
  # a 1 (an identity map)
  if isempty(ccoeff)
    push!(ccoeff,1)
  end
  # if only one entry, add a zero for c₀
  if length(ccoeff) == 1
    push!(ccoeff,0)
  end
  ncoeff = length(ccoeff)-2

  # Coefficients of |z(ζ)|²
  # dcoeff[1] = d₀, dcoeff[2] = d₋₁, etc.
  # Note that d₁ = conj(d₋₁) = conj(dcoeff[2]), d₂ = conj(d₋₂) = conj(dcoeff[3])
  dcoeff = [ccoeff'ccoeff]
  for k = 1:ncoeff+1
      push!(dcoeff,ccoeff[1:end-k]'ccoeff[k+1:end])
  end
  ps = PowerSeries(ccoeff,dcoeff)
  dps = PowerSeriesDerivatives(ps)

  ζ = circle(N)
  z = ps(ζ)
  dz, ddz = dps(ζ)

  area, Zc, J = shape_moments(ps)

  Ma = addedmass(ps,area)

  PowerMap(ps, dps, N, ζ, z, dz, area, Zc, J, Ma)
end

function Base.show(io::IO, m::PowerMap)
    println(io, "Power series map")
end

function Base.summary(m::PowerMap)
  m.ps
end

(m::PowerMap)(ζ) = m.ps(ζ)

(minv::InverseMap{PowerMap})(z::Vector{ComplexF64}) =
        evalinv_exterior(z,minv.m.ps,minv.m.dps)

(minv::InverseMap{PowerMap})(z::ComplexF64) = getindex(minv([z]),1)

(dm::DerivativeMap{PowerMap})(ζ) = dm.m.dps(ζ)

function shape_moments(ps::PowerSeries)

  ncoeff = length(ps.ccoeff)-2
  k = -1:ncoeff
  l = -1:ncoeff
  kml = repeat(transpose(k), length(ps.ccoeff))-repeat(l, 1, length(ps.ccoeff))

  c = Reflect(ShiftReindex(ps.ccoeff,-2))
  d = Reflect(OddReindex(ps.dcoeff))

  area = -π*sum(k.*abs.(c(-k)).^2)

  if area > 0
    Zc = -π/area*sum(k.*c(-k).*d(k))
  else
    Zc = ComplexF64(0)
  end
  J = real(-0.5π*c(-k)'*d(-kml)*(l.*c(-l)))

  return area, Zc, J
end

#=  Karman-Trefftz map from unit circle to airfoil
  This is actually a composition of maps, first from the unit circle in the ζ
  plane to a circle of radius a at center ϵCexp(iδ) in the ζ̃ plane. The mapping is
  (z-νC)/(z+νC) = ((ζ̃-C)/(ζ̃+C))^ν. This gives an interior angle at the
  trailing edge (ζ̃ = C) of angle (2-ν)π. Note that ζ̃ = ϵCexp(iδ) + aζ.
  The basic Joukowski airfoil is specified by ν = 2.

=#
struct KarmanTrefftzMap <: ConformalMap

  "K-T power ν, generally <∼ 2"
  nu::Real

  "Displacement of circle from origin, as fraction of C. Should be smaller than ~0.2."
  epsilon::Real

  "Angular orientation of circle, π/2 <= δ <= π"
  delta::Real

  "Scaling factor C. Airfoil chord is approximately 4C"
  C::Real

  "number of plotting control points"
  N::Int

  "radius of circle in ζ̃ plane"
  a::Real

  "center of circle in ζ̃ plane"
  mu::Number

  "Power series"
  ps::PowerSeries

  "control point coordinates in circle space"
  ζ::Vector{ComplexF64}

  "control point coordinates in body-fixed space"
  z::Vector{ComplexF64}

  "map Jacobian in body-fixed coordinates"
  dzdζ::Vector{ComplexF64}

  "Area enclosed by the mapped shape"
  area      :: Float64

  "Centroid of the mapped shape"
  Zc        :: ComplexF64

  "2nd area moment of the mapped shape"
  J         :: Float64

  "Added mass tensor"
  Ma        :: Array{Float64,2}

end

@doc raw"""
    KarmanTrefftzMap(ν,ϵ,δ,C[;N = 200]) <: ConformalMap

Create a map from the exterior of the unit
circle to the exterior of a Karman-Trefftz airfoil.

The form of the mapping is

```math
\frac{z-\nu C}{z+\nu C}  =
\left(\frac{\tilde{\zeta}-C}{\tilde{\zeta}+C}\right)^\nu
```
where $\tilde{\zeta}$ are the coordinates in an intermediate plane, in which
the circle is of radius ``a`` and centered at ``\epsilon C e^{i\delta}``:
```math
\tilde{\zeta} = \epsilon C e^{i\delta} + a \zeta
```
Note that ``a/C \geq 1`` and is determined by the choices for ``\epsilon`` and ``\delta``.

The trailing edge angle, ``(2-\nu)\pi`` is specified by ``\nu``. The thickness
is controlled by ``\epsilon C \cos\delta`` and the camber by ``\epsilon C \sin\delta``.
The airfoil chord length is approximately ``4C``. Generally, ``\epsilon`` should be
much smaller than 1 and ``\delta`` between ``\pi/2`` and ``\pi``.

The resulting map `m` can be evaluated at a single or a vector of points `ζ`
with `m(ζ)`.

# Example

```jldoctest
julia> ν = 1.9; ϵ = 0.1; δ = π; C = 0.25;

julia> m = KarmanTrefftzMap(ν,ϵ,δ,C)
Karman-Trefftz map

julia> ζ = [1.0+3.0im,-2.0-2.0im,0.0+1.1im];

julia> m(ζ)
3-element Array{Complex{Float64},1}:
   0.268188+0.764722im
  -0.624265-0.502634im
 -0.0390996+0.126737im
```
""" KarmanTrefftzMap
function KarmanTrefftzMap(ν,ϵ,δ,C; N::Int = 200)

  β = atan(ϵ*sin(δ)/(1-ϵ*cos(δ)))
  a_C = cos(β) + sqrt(ϵ^2-sin(β)^2)
  a = a_C*C
  μ = ϵ*C*exp(im*δ)

  # Compute power series coefficients
  ncoeff = 4
  ccoeff = zeros(ComplexF64,ncoeff)
  ccoeff[1] = a
  ccoeff[2] = μ
  ccoeff[3] = (ν^2-1)*C/(3*a_C)
  ccoeff[4] = -ccoeff[3]*μ/a

  # Coefficients of |z(ζ)|²
  # dcoeff[1] = d₀, dcoeff[2] = d₋₁, etc.
  # Note that d₁ = conj(d₋₁) = conj(dcoeff[2]), d₂ = conj(d₋₂) = conj(dcoeff[3])
  dcoeff = [dot(ccoeff,ccoeff)]
  for k = 1:ncoeff+1
      push!(dcoeff,dot(ccoeff[1:end-k],ccoeff[k+1:end]))
  end
  ps = PowerSeries(ccoeff,dcoeff)

  ζ = circle(N)
  ζ̃ = similar(ζ)
  z = similar(ζ)
  dzdζ = similar(ζ)

  @. ζ̃ = μ + a*ζ
  @. z = ν*C*((ζ̃+C)^ν+(ζ̃-C)^ν)/((ζ̃+C)^ν-(ζ̃-C)^ν)
  @. dzdζ = 4ν^2*C^2*a*(ζ̃^2-C^2)^(ν-1)/((ζ̃+C)^ν-(ζ̃-C)^ν)^2

  area, Zc, J = shape_moments(ps)

  Ma = addedmass(ps,area)

  KarmanTrefftzMap(ν,ϵ,δ,C,N,a,μ,ps,ζ,z,dzdζ,area, Zc, J, Ma)
end

JoukowskiMap(y...;x...) = KarmanTrefftzMap(2,y...;x...)

function (m::KarmanTrefftzMap)(ζ::Number)
  ζ̃ = m.mu + m.a*ζ
  return m.nu*m.C*((ζ̃+m.C)^m.nu+(ζ̃-m.C)^m.nu)/((ζ̃+m.C)^m.nu-(ζ̃-m.C)^m.nu)
end

function (minv::InverseMap{KarmanTrefftzMap})(z::Number)
  ζ̃ = minv.m.C*((z+minv.m.nu*minv.m.C)^(1/minv.m.nu)+(z-minv.m.nu*minv.m.C)^(1/minv.m.nu))/
           ((z+minv.m.nu*minv.m.C)^(1/minv.m.nu)-(z-minv.m.nu*minv.m.C)^(1/minv.m.nu))
  return  (ζ̃ - minv.m.mu)/minv.m.a
end

function (dm::DerivativeMap{KarmanTrefftzMap})(ζ::Number)
  ζ̃ = dm.m.mu + dm.m.a*ζ
  return 4dm.m.nu^2*dm.m.C^2*dm.m.a*(ζ̃^2-dm.m.C^2)^(dm.m.nu-1)/
                      ((ζ̃+dm.m.C)^dm.m.nu-(ζ̃-dm.m.C)^dm.m.nu)^2
end

(m::KarmanTrefftzMap)(ζ::Vector{T}) where T<:Number = m.(ζ)

(minv::InverseMap{KarmanTrefftzMap})(z::Vector{T}) where T<:Number = minv.(z)

(dm::DerivativeMap{KarmanTrefftzMap})(ζ::Vector{T}) where T<:Number = dm.(ζ)

function Base.show(io::IO, m::KarmanTrefftzMap)
    if m.nu == 2
      println(io, "Joukowski map")
    else
      println(io, "Karman-Trefftz map")
    end
end

function Base.summary(m::KarmanTrefftzMap)
  println("Karman-Trefftz map of unit circle to exterior of airfoil")
  print("   ")
  print("circle center: ")
  println("($(round(m.epsilon*cos(m.delta),4)),$(round(m.epsilon*sin(m.delta),4)))")
  print("   ")
  print("chord length: ")
  println("$(4*m.C)")
  print("   ")
  print("trailing edge angle (deg): ")
  println("$(round((2-m.nu)*180,4))")
end

#=   Exterior map from circle to polygon  =#

struct ExteriorMap <: ConformalMap

  "Number of vertices on polygon"
  N :: Int

  "Coordinates of vertices on the polygon, defined ccw"
  z :: Vector{ComplexF64}

  "Interior angles in the polygon"
  angle  :: Vector{Float64}

  "Quadrature nodes and weights for the Gauss-Jacobi"
  qdata  :: Tuple{Array{Float64,2},Array{Float64,2}}

  "Accuracy of the quadrature"
  accuracy  :: Float64

  "Coordinates of the prevertices on the unit circle (interior)"
  ζ :: Vector{ComplexF64}

  "Constant factor of the mapping"
  constant  :: ComplexF64

  "Coordinates of the pre-prevertices on the unit circle (exterior)"
  preprev   :: Vector{ComplexF64}

  "Angular locations of the pre-prevertices on the unit circle (exterior)"
  prevangle :: Vector{Float64}

  "Number of multipole coefficients"
  ncoeff :: Int

  "Multipole coefficients for z(ζ) and |z(ζ)|²"
  ps    :: PowerSeries

  "Moments of prevertices"
  mom       :: Vector{ComplexF64}

  "Area enclosed by the mapped polygon"
  area      :: Float64

  "Centroid of the mapped polygon"
  Zc        :: ComplexF64

  "2nd area moment of the mapped polygon"
  J         :: Float64

  "Added mass tensor"
  Ma        :: Array{Float64,2}

end

"""
    ExteriorMap(p::Polygon[;tol::Float64][,ncoeff::Int]) <: ConformalMap

Create a Schwarz-Christoffel map from the interior or exterior of
the unit circle to the exterior of polygon `p`.

# Example

```jldoctest
julia> p = Polygon([-1.0,0.2,1.0,-1.0],[-1.0,-1.0,0.5,1.0]);

julia> m = ExteriorMap(p)
Schwarz-Christoffel map of unit circle to exterior of polygon with 4 vertices
```

`ExteriorMap(p;tol=1e-12)` manually sets the tolerance to `1e-12`
(the default is 1e-8).

`ExteriorMap(p;ncoeff=200)` manually sets the number of coefficients of
negative powers of the multipole expansion of the mapping to `200`
(the default is 100).

The resulting map `m` can be evaluated at a single or vector of points `ζ`
with `m(ζ[;inside::Bool])`. The points are assumed to lie outside the unit
circle, unless the optional argument `inside=true`, in which case they are
assumed to lie inside the circle.

# Example

```jldoctest
julia> p = Polygon([-1.0,0.2,1.0,-1.0],[-1.0,-1.0,0.5,1.0]);

julia> m = ExteriorMap(p);

julia> ζ = [0.1,0.5-0.75im,-0.25-0.3im];

julia> m(ζ;inside=true)
3-element Array{Complex{Float64},1}:
   -6.9344-7.68965im
 0.0439774-1.11249im
   2.41181-0.044779im

julia> ζ = [1.0+3.0im,-2.0-2.0im,0.0+1.1im];

julia> m(ζ)
3-element Array{Complex{Float64},1}:
   0.81614+3.02956im
  -2.25237-2.08523im
 -0.333104+0.975837im
```
"""
function ExteriorMap(p::Polygon;tol::Float64 = 1e-8,ncoeff::Int = 100)

  n = length(p.vert)

  zeta0 = ComplexF64[]

  w = reverse(vertex(p), dims = 1)
  beta = 1 .- reverse(interiorangle(p), dims = 1)

  # do some fixing
  n = length(p)
  renum = 1:n
  shift = [2:n;1]
  w = reverse(vertex(p), dims = 1)
  beta = 1 .- reverse(interiorangle(p), dims = 1)
  while any(abs.(beta[n] .- [0;1]) .< eps()) & (n > 2)
    renum = renum[shift]
    w = w[shift]
    beta = beta[shift]
    if renum[1]==1
      deg = abs.(beta) .< eps()
      w[deg] = []
      beta[deg] = []
      renum = 1:2
      n = 2
      break
    end
  end

  p = Polygon(reverse(w, dims = 1),1 .- reverse(beta, dims = 1))

  nqpts = max(ceil(Int,-log10(tol)),2)
  qdat = qdata(beta,nqpts)

  zeta, c = param(w,beta,zeta0,qdat)
  preprev = -c/abs(c)./reverse(zeta, dims = 1)
  prevangle = angle.(preprev)

  # first two entries are for c₁ and c₀.
  betaflip = reverse(beta, dims = 1)
  mom = [sum(betaflip.*preprev)]
  for k = 1:ncoeff
    push!(mom,sum(betaflip.*preprev.^(k+1)))
  end

  # use numerical quadrature to compute the first two coefficients
  # the first coefficient should be abs(c), but we don't set this
  # explicitly
  nθ = 1024
  dθ = 2π/nθ
  ζs = exp.(im*collect(0:nθ-1)*dθ)
  σs = -c/abs(c)./ζs
  dc = evaluate_exterior(σs,w,beta,zeta,c,qdat).*sin(dθ)
  dc ./= ζs
  ccoeff = ComplexF64[]
  for k = -1:0
    push!(ccoeff,sum(dc))
    dc .*= ζs
  end
  ccoeff ./= 2π

  # now, we use Gauss-Jacobi quadrature for the remaining coefficients.
  ζ = [zeta;zeta[1]]
  N = n+1
  argz1 = angle.(ζ[1:N-1])
  argz2 = angle.(ζ[2:N])
  ang21 = angle.(ζ[2:N]./ζ[1:N-1])
  discont = (argz2-argz1).*ang21 .< 0
  argz2[discont] += 2π*sign.(ang21[discont])
  mid = exp.(im*0.5*(argz1+argz2))
  qdat2 = qdata(beta,max(ceil(Int,ncoeff/2),nqpts))
  dabsquad = DabsQuad(beta,qdat2)
  for k = 1:ncoeff
     push!(ccoeff,c*(-c/abs(c))^k*sum(dabsquad(ζ[1:N-1],mid,collect(1:N-1),zeta,-k)+
                  dabsquad(ζ[2:N],mid,[collect(2:N-1);1],zeta,-k))/(2π*k))
  end

  # Now compute the coefficients of the square mapping
  # first entry is d₀
  dcoeff = [ccoeff'ccoeff]
  for k = 1:ncoeff+1
    push!(dcoeff,ccoeff[1:end-k]'ccoeff[k+1:end])
  end
  ps = PowerSeries(ccoeff,dcoeff)

  area, Zc, J = shape_moments(vertex(p))

  Ma = addedmass(ps,area)

  ExteriorMap(n,p.vert,p.angle,qdat,tol,zeta,c,preprev,prevangle,
              ncoeff,ps,mom,area,Zc,J,Ma)
end

function shape_moments(z::Vector{ComplexF64})

  zmid = 0.5*(z+circshift(z,-1))
  dz = circshift(z,-1)-z
  area = 0.5*imag(sum(conj.(zmid).*dz))
  if area > 0
    Zc = -0.5im/area*sum(dz.*(abs.(zmid).^2+abs.(dz).^2/12))
  else
    Zc = mean(z)
  end
  J = 0.25*imag(sum(conj.(zmid).*dz.*(abs.(zmid).^2+abs.(dz).^2/12)))

  return area, Zc, J
end

function Base.show(io::IO, m::ExteriorMap)
    println(io, "Schwarz-Christoffel map of unit circle to exterior of"*
                " polygon with $(m.N) vertices")
end

function Base.summary(m::ExteriorMap)

  println("Schwarz-Christoffel map of unit circle to exterior"*
          " of polygon with $(m.N) vertices")
  print("   ")
  print("vertices: ")
  for i = 1:length(m.z)
    print("($(round(real(m.z[i]), digits=4)),$(round(imag(m.z[i]), digits=4))), ")
  end
  println()
  print("   ")
  print("interior angles/π: ")
  for i = 1:length(m.angle)
    print("$(round(m.angle[i], digits=4)), ")
  end
  println()
  print("   ")
  print("prevertices on circle: ")
  for i = length(m.ζ):-1:1
    print("($(round(real(m.ζ[i]), digits=4)),$(round(imag(m.ζ[i]), digits=4))), ")
  end
  println()
  print("   ")
  print("prevertex angles/π: ")
  for i = 1:length(m.prevangle)
    print("$(round(m.prevangle[i]/π, digits=4)), ")
  end
  println()
  print("   ")
  print("constant = $(round(m.constant, digits=4)), ")
  print("accuracy = $(m.accuracy), ")
  println()
  print("   ")
  print("number of multipole coefficients = $(m.ncoeff)")
  println()

end

function (m::ExteriorMap)(ζ::Vector{ComplexF64};inside::Bool=false)
  if inside
    return evaluate_exterior(ζ,reverse(m.z, dims = 1),1 .- reverse(m.angle, dims = 1),
            m.ζ,m.constant,m.qdata)
  else
    b = -m.constant/abs(m.constant)
    ζ[ζ.==0] .= eps();
    ζ[abs.(ζ).<1] = ζ[abs.(ζ).<1]./abs.(ζ[abs.(ζ).<1])

    σ = b./ζ
    return evaluate_exterior(σ,reverse(m.z, dims = 1),1 .- reverse(m.angle, dims = 1),
            m.ζ,m.constant,m.qdata)
  end

end

(m::ExteriorMap)(ζ::ComplexF64;inside::Bool=false) = getindex(m([ζ];inside=inside),1)

"""
    InverseMap(m::ConformalMap)

Constructs the inverse conformal map of the conformal map `m`.

This inverse conformal map can be evaluated at a single or vector of points.
Points should be outside the body. Whether the resulting point in the circle
plane is interpreted inside or outside the circle is determined by the optional
argument `inside`, which defaults to `false`.

# Example

```jldoctest
julia> p = Polygon([-1.0,0.2,1.0,-1.0],[-1.0,-1.0,0.5,1.0]);

julia> m = ExteriorMap(p);

julia> m⁻¹ = InverseMap(m);

julia> ζ = [1.0+3.0im,-2.0-2.0im,0.1+1.1im];

julia> m⁻¹(m(ζ))
3-element Array{Complex{Float64},1}:
  1.0+3.0im
 -2.0-2.0im
  0.1+1.1im
```
""" InverseMap

function (minv::InverseMap{ExteriorMap})(z::Vector{ComplexF64};inside::Bool=false)
  if inside
    return evalinv_exterior(z,reverse(minv.m.z, dims = 1),1 .- reverse(minv.m.angle, dims = 1),
            minv.m.ζ,minv.m.constant,minv.m.qdata)
  else
    σ = evalinv_exterior(z,reverse(minv.m.z, dims = 1),1 .- reverse(minv.m.angle, dims = 1),
            minv.m.ζ,minv.m.constant,minv.m.qdata)
    b = -minv.m.constant/abs(minv.m.constant)

    return b./σ
  end

end

(minv::InverseMap{ExteriorMap})(z::ComplexF64;inside::Bool=false) =
                getindex(minv([z];inside=inside),1)

"""
    DerivativeMap(m::ConformalMap)

Constructs new conformal maps from the first and second derivatives of the
conformal map `m`.

These new conformal maps can be evaluated at a single or vector of points just as
 `m` is. The first entry in the tuple returned is the first derivative,
the second entry is the second derivative.

# Example

```jldoctest
julia> p = Polygon([-1.0,0.2,1.0,-1.0],[-1.0,-1.0,0.5,1.0]);

julia> m = ExteriorMap(p);

julia> dm = DerivativeMap(m);

julia> ζ = [0.1,0.5-0.75im,-0.25-0.3im];

julia> dz, ddz = dm(ζ;inside=true);

julia> dz
3-element Array{Complex{Float64},1}:
  67.2068+76.6284im
 -1.11666+0.544576im
  3.99129-5.30641im
```
"""
function DerivativeMap() end

function (dm::DerivativeMap{ExteriorMap})(ζ::Vector{ComplexF64};inside::Bool=false)
  if inside
    return evalderiv_exterior(ζ,1 .- reverse(dm.m.angle, dims = 1),dm.m.ζ,dm.m.constant)
  else
    b = -dm.m.constant/abs(dm.m.constant)
    ζ[ζ.==0] .= eps();
    ζ[abs.(ζ).<1] = ζ[abs.(ζ).<1]./abs.(ζ[abs.(ζ).<1])

    σ = b./ζ
    dσ = -σ./ζ
    ddσ = -2.0*dσ./ζ
    dz, ddz = evalderiv_exterior(σ,1 .- reverse(dm.m.angle, dims = 1),dm.m.ζ,dm.m.constant)
    ddz = ddz.*dσ.^2 + dz.*ddσ
    dz .*= dσ
    return dz, ddz
  end
end

function (dm::DerivativeMap{ExteriorMap})(ζ::ComplexF64;inside::Bool=false)
  dz, ddz = dm([ζ];inside=inside)
  return dz[1],ddz[1]
end

#= get various data about the map =#
"""
    summary(m::ConformalMap)

Returns a summary of data for a conformal map

# Example

```jldoctest
julia> p = Polygon([-1.0,0.2,1.0,-1.0],[-1.0,-1.0,0.5,1.0]);

julia> m = ExteriorMap(p);

julia> summary(m)
Schwarz-Christoffel map of unit circle to exterior of polygon with 4 vertices
   vertices: (-1.0,-1.0), (0.2,-1.0), (1.0,0.5), (-1.0,1.0),
   interior angles/π: 0.5, 0.656, 0.422, 0.422,
   prevertices on circle: (1.0,0.0), (0.3764,-0.9265), (-0.9024,-0.4309), (-0.1868,0.9824),
   prevertex angles/π: -0.7291, -0.3519, 0.1291, 0.7111,
   constant = 0.6722 + 0.7669im, accuracy = 1.0e-8,
   number of multipole coefficients = 100
```

"""
Base.summary

"""
    length(m::ConformalMap) -> Integer

Returns the number of control points/vertices of the map `m`.

# Example

```jldoctest
julia> p = Polygon([-1.0,0.2,1.0,-1.0],[-1.0,-1.0,0.5,1.0]);

julia> m = ExteriorMap(p);

julia> length(m)
4
```
"""
Base.length(m::ConformalMap) = m.N

"""
    parameters(m::ExteriorMap) -> Tuple{Vector{ComplexF64},ComplexF64}

Returns a tuple of a vector of the prevertices and the complex factor of
the exterior polygon mapping `m`.

# Example

```jldoctest paramtest
julia> p = Polygon([-1.0,0.2,1.0,-1.0],[-1.0,-1.0,0.5,1.0]);

julia> m = ExteriorMap(p);

julia> prev, C = parameters(m);

julia> prev
4-element Array{Complex{Float64},1}:
       1.0+0.0im
  0.376406-0.926455im
 -0.902383-0.430935im
 -0.186756+0.982406im
```
"""
parameters(m::ExteriorMap) = reverse(m.ζ, dims = 1), m.constant

"""
    coefficients(m::ConformalMap) -> Tuple{Vector{ComplexF64},Vector{ComplexF64}}

Returns a tuple of vectors of the complex coefficients of the multipole
expansion of the mapping \$z(\\zeta)\$ described by `m` as well as the
coefficients of the square magnitude of the mapping \$|z(\\zeta)|^2\$.

# Example

```jldoctest
julia> p = Polygon([-1.0,0.2,1.0,-1.0],[-1.0,-1.0,0.5,1.0]);

julia> m = ExteriorMap(p);

julia> ccoeff, dcoeff = coefficients(m);
```
"""
coefficients(m::ConformalMap) = m.ps.ccoeff, m.ps.dcoeff

"""
    moments(m::ExteriorMap) -> Vector{ComplexF64}

Return the moments of the prevertices for exterior polygon mapping `m`.

# Example

```jldoctest
julia> p = Polygon([-1.0,0.2,1.0,-1.0],[-1.0,-1.0,0.5,1.0]);

julia> m = ExteriorMap(p);

julia> mom = moments(m);
```
"""
moments(m::ExteriorMap) = m.mom

"""
    area(m::ConformalMap) -> Float64

Returns the area of the shape described by the mapping `m`.

# Example

```jldoctest
julia> p = Polygon([-1.0,0.2,1.0,-1.0],[-1.0,-1.0,0.5,1.0]);

julia> m = ExteriorMap(p);

julia> area(m)
2.9
```

```jldoctest
julia> c = ComplexF64[1];

julia> m = PowerMap(c);

julia> area(m)
3.141592653589793
```

"""
area(m::ConformalMap) = m.area

"""
    centroid(m::ConformalMap) -> ComplexF64

Returns the complex centroid position of the shape described by the
mapping `m`.

# Example

```jldoctest
julia> p = Polygon([-1.0,0.2,1.0,-1.0],[-1.0,-1.0,0.5,1.0]);

julia> m = ExteriorMap(p);

julia> centroid(m)
-0.20919540229885059 - 0.04022988505747128im
```
"""
centroid(m::ConformalMap) = m.Zc

"""
    Jmoment(m::ConformalMap) -> Float64

Returns the second area moment of the shape described by the
mapping `m`.

# Example

```jldoctest
julia> p = Polygon([-1.0,0.2,1.0,-1.0],[-1.0,-1.0,0.5,1.0]);

julia> m = ExteriorMap(p);

julia> Jmoment(m)
1.5768333333333333
```
"""
Jmoment(m::ConformalMap) = m.J

"""
    addedmass(m::ConformalMap) -> Array{Float64,2}

Returns the added mass matrix of the shape described by the
conformal mapping `m`.

# Example

```jldoctest
julia> p = Polygon([-1.0,0.2,1.0,-1.0],[-1.0,-1.0,0.5,1.0]);

julia> m = ExteriorMap(p);

julia> addedmass(m)
3×3 Array{Float64,2}:
  0.725129    0.0944902  -1.37387
  0.0944902   3.67634    -0.255119
 -1.37387    -0.255119    3.59231
```
"""
addedmass(m::ConformalMap) = m.Ma

function addedmass(ps::PowerSeries,area::Float64)
  M = zeros(3,3)
  c = Reflect(ShiftReindex(ps.ccoeff,-2))
  d = Reflect(OddReindex(ps.dcoeff))

  k = 1:length(ps.ccoeff)+1
  M[1,1] = 0.5sum(k.*abs.(d(-k)).^2)
  M[1,2] = M[2,1] = -imag(c(1)*d(-1))
  M[1,3] = M[3,1] =  real(c(1)*d(-1))
  M[2,3] = M[3,2] = -imag(c(1)*c(-1))
  M[2,2] = abs(c(1))^2-real(c(1)*c(-1))-area/(2π)
  M[3,3] = abs(c(1))^2+real(c(1)*c(-1))-area/(2π)
  return 2π*M

end

end


text:
function part1(input)
    house = (0, 0)
    lucky_children = Set((house,))
    directions = Dict('^' => (0, -1), 'v' => (0, 1), '>' => (1, 0), '<' => (-1, 0))
    for c in collect(readchomp(input))
        house = house .+ directions[c]
        push!(lucky_children, house)
    end
    return length(lucky_children)
end

function part2(input)
    house = (0, 0)
    robo_house = (0, 0)
    lucky_children = Set((house,))
    directions = Dict('^' => (0, -1), 'v' => (0, 1), '>' => (1, 0), '<' => (-1, 0))
    for (i, c) in enumerate(collect(readchomp(input)))
        if i & 1 == 1
            house = house .+ directions[c]
            push!(lucky_children, house)
        else
            robo_house = robo_house .+ directions[c]
            push!(lucky_children, robo_house)
        end
    end
    return length(lucky_children)
end

ext: jl


text:
module Day16

using ..ReTest, ..DocStringExtensions

import ..solve, ..PuzzleInfo, ..TEST_DATA_DIRECTORY

abstract type Day16PuzzleInfo <: PuzzleInfo end

"""
Day 16 part 1 puzzle information.

$(FIELDS)
"""
struct Part1PuzzleInfo <: Day16PuzzleInfo
    """Path to input data file."""
    datapath::String
end

"""
Day 16 part 2 puzzle information.

$(FIELDS)
"""
struct Part2PuzzleInfo <: Day16PuzzleInfo
    """Path to input data file."""
    datapath::String
end

"""
    $(FUNCTIONNAME)(input)

Read Buoyancy Interchange Transmission System (BITS) packets from input stream `input`.
"""
function readbitspackets(input::IO)
    hexadecimalstring = readline(input)
    bytes = hex2bytes(hexadecimalstring)
    BitVector(reduce(vcat, reverse(b) for b in digits.(bytes, base=2, pad=8)))
end

"""
    $(FUNCTIONNAME)(bits)

Get the integer value correspond to a bit vector `bits` with most significant bit first.
"""
function bits2int(bits)
    powerof2 = 1
    intval = 0
    for b in Iterators.reverse(bits)
        intval += b * powerof2
        powerof2 <<= 1
    end
    intval
end

"""Buoyancy Interchange Transmission System (BITS) packet."""
abstract type Packet end

"""
Buoyancy Interchange Transmission System (BITS) literal packet.

$(FIELDS)
"""
struct LiteralPacket <: Packet
    """Version number for packet."""
    version::Int
    """Value of packet."""
    value::Int
end

"""
Buoyancy Interchange Transmission System (BITS) operator packet.

$(FIELDS)
"""
struct OperatorPacket <:Packet
    """Version number for packet."""
    version::Int
    """Type ID for packet."""
    typeid::Int
    """Subpackets of packet."""
    subpackets::Vector{Packet}
end

OperatorPacket(version, operator) = OperatorPacket(version, operator, [])

"""Map from operator packet type IDs to corresponding operators."""
const operatormap = Dict{Int, Function}(
    0 => sum,
    1 => prod,
    2 => minimum,
    3 => maximum,
    5 => x -> (x[1] > x[2]),
    6 => x -> (x[1] < x[2]),
    7 => x -> (x[1] == x[2])
)

"""
    $(FUNCTIONNAME)(bits)

Process packet represented by leading elements BITS transmission `bits`. Returns
tuple consisting of packet and remaining elements of `bits` after packet bits removed.
"""
function processpacket(bits)
    version, bits = bits2int(bits[1:3]), bits[4:end]
    typeid, bits = bits2int(bits[1:3]), bits[4:end]
    if typeid == 4
        literalbits = BitVector()
        while bits[1] != 0
            append!(literalbits, bits[2:5])
            bits = bits[6:end]
        end
        append!(literalbits, bits[2:5])
        packet = LiteralPacket(version, bits2int(literalbits))
        bits = bits[6:end]
    else
        packet = OperatorPacket(version, typeid)
        if bits[1] == 0
            subpacketslength, bits = bits2int(bits[2:16]), bits[17:end]
            subpacketbits, bits = bits[1:subpacketslength], bits[subpacketslength+1:end]
            while !isempty(subpacketbits)
                subpacket, subpacketbits = processpacket(subpacketbits)
                push!(packet.subpackets, subpacket)
            end
        else
            subpacketsnumber, bits = bits2int(bits[2:12]), bits[13:end]
            for p in 1:subpacketsnumber
                subpacket, bits = processpacket(bits)
                push!(packet.subpackets, subpacket)
            end
        end
    end
    packet, bits
end

"""
    $(FUNCTIONNAME)(packet)

Evaluate the value represented by `packet`.
"""
function getvalue(::Packet) end

getvalue(packet::LiteralPacket) = packet.value

function getvalue(packet::OperatorPacket)
    operatormap[packet.typeid]([getvalue(subpacket) for subpacket in packet.subpackets])
end

"""
    $(FUNCTIONNAME)(packet)

Recursively compute sum of version of `packet` and all subpackets.
"""
function getversionsum(packet::Packet)
    if isa(packet, OperatorPacket)
        packet.version + sum(getversionsum(subpacket) for subpacket in packet.subpackets)
    else
        packet.version
    end
end

function solve(info::Day16PuzzleInfo)
    bits = open(readbitspackets, info.datapath)
    packet, bits = processpacket(bits)
    @assert !any(bits) "Expected remaining bits to all be zero"
    isa(info, Part1PuzzleInfo) ? getversionsum(packet) : getvalue(packet)
end

@testset begin
    @test solve(Part1PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day16-1.txt"))) == 16
    @test solve(Part1PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day16-2.txt"))) == 12
    @test solve(Part1PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day16-3.txt"))) == 23
    @test solve(Part1PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day16-4.txt"))) == 31
    @test solve(Part2PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day16-5.txt"))) == 3
    @test solve(Part2PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day16-6.txt"))) == 54
    @test solve(Part2PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day16-7.txt"))) == 7
    @test solve(Part2PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day16-8.txt"))) == 9
    @test solve(Part2PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day16-9.txt"))) == 1
    @test solve(Part2PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day16-10.txt"))) == 0
    @test solve(Part2PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day16-11.txt"))) == 0
    @test solve(Part2PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day16-12.txt"))) == 1
end

end # module


text:
#
# This file is part of the SlowActors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

using Test, SafeTestsets

function redirect_devnull(f)
    open(@static(Sys.iswindows() ? "nul" : "/dev/null"), "w") do io
        redirect_stdout(io) do
            f()
        end
    end
end

@safetestset "Basics"       begin include("test_basics.jl") end
@safetestset "Interface"    begin include("test_interface.jl") end

println("running examples, output suppressed!")
redirect_devnull() do
    @safetestset "Factorial"     begin include("../examples/factorial.jl") end
    @safetestset "Fib"           begin include("../examples/fib.jl") end
    @safetestset "Simple"        begin include("../examples/simple.jl") end
    @safetestset "Stack"         begin include("../examples/stack.jl") end
end


text:
using NumericalMethodsforEngineers
using Test

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  function f3(x::Float64, y::Vector{Float64})
    [y[2], -2y[2]/y[1]]
  end

  (iters, nsteps, converged, res) = shootingmethod(f3, [0.0, 1.0, 1.0, 2.0], [1.0, 3.0],
    10, 0.000001, 80)

  if !converged
    println("No conversion achieved.")
  elseif iters > 50
    println("Iter limit reached.")
  else
    @test round.(res[:, 2], digits=4)' ==
      [1.0  1.1742  1.321  1.4471  1.5569  1.6534  1.7389  1.8151  1.8833  1.9446  2.0]
    #println()
    #res |> display
  end

end


text:
import Images.imread, Images.imwrite

function imread{S<:IO}(stream::S, ::Type{Images.Dummy})
    pixels = reshape(uint8(1:12), 3, 4)
    Image(pixels, Compat.@compat Dict("colorspace" => "Gray", "spatialorder" => "xy"))
end

function imwrite(img, filename::String, ::Type{Images.Dummy})
    open(filename, "w") do stream
        println(stream, "Dummy Image")
    end
end


text:
# Julia wrapper for header: /usr/include/libavutil/encryption_info.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_encryption_info_alloc,
    av_encryption_info_clone,
    av_encryption_info_free,
    av_encryption_info_get_side_data,
    av_encryption_info_add_side_data,
    av_encryption_init_info_alloc,
    av_encryption_init_info_free,
    av_encryption_init_info_get_side_data,
    av_encryption_init_info_add_side_data

function av_encryption_info_alloc(subsample_count::Integer, key_id_size::Integer, iv_size::Integer)
    ccall((:av_encryption_info_alloc, libavutil), Ptr{AVEncryptionInfo}, (UInt32, UInt32, UInt32), subsample_count, key_id_size, iv_size)
end

function av_encryption_info_clone(info)
    ccall((:av_encryption_info_clone, libavutil), Ptr{AVEncryptionInfo}, (Ptr{AVEncryptionInfo},), info)
end

function av_encryption_info_free(info)
    ccall((:av_encryption_info_free, libavutil), Cvoid, (Ptr{AVEncryptionInfo},), info)
end

function av_encryption_info_get_side_data(side_data, side_data_size::Csize_t)
    ccall((:av_encryption_info_get_side_data, libavutil), Ptr{AVEncryptionInfo}, (Ptr{UInt8}, Csize_t), side_data, side_data_size)
end

function av_encryption_info_add_side_data(info, side_data_size)
    ccall((:av_encryption_info_add_side_data, libavutil), Ptr{UInt8}, (Ptr{AVEncryptionInfo}, Ptr{Csize_t}), info, side_data_size)
end

function av_encryption_init_info_alloc(system_id_size::Integer, num_key_ids::Integer, key_id_size::Integer, data_size::Integer)
    ccall((:av_encryption_init_info_alloc, libavutil), Ptr{AVEncryptionInitInfo}, (UInt32, UInt32, UInt32, UInt32), system_id_size, num_key_ids, key_id_size, data_size)
end

function av_encryption_init_info_free(info)
    ccall((:av_encryption_init_info_free, libavutil), Cvoid, (Ptr{AVEncryptionInitInfo},), info)
end

function av_encryption_init_info_get_side_data(side_data, side_data_size::Csize_t)
    ccall((:av_encryption_init_info_get_side_data, libavutil), Ptr{AVEncryptionInitInfo}, (Ptr{UInt8}, Csize_t), side_data, side_data_size)
end

function av_encryption_init_info_add_side_data(info, side_data_size)
    ccall((:av_encryption_init_info_add_side_data, libavutil), Ptr{UInt8}, (Ptr{AVEncryptionInitInfo}, Ptr{Csize_t}), info, side_data_size)
end


text:
@testset "Adversarial problem tests" begin
    @test isapprox(adversarialProblem([1, 3], [3, 1], [2, 2], 2, getKnapsackConstraints([1, 2], 1), 0.5, KnapsackProblemDescriptor(2)), 5.0, atol = 1e-5)

    @test isapprox(adversarialProblem([1, 3], [3, 1], [2, 2], 2, getKnapsackConstraints([1, 2], 1), 1.0, KnapsackProblemDescriptor(2)), 3.99988, atol = 1e-5)
end


text:
export TFIDF

"""
    TFIDF(
        data::DataAccessor,
        tf::AbstractMatrix,
        idf::AbstractMatrix
    )

Content-based recommendation using TF-IDF scoring. TF and IDF matrix are respectively specified as `tf` and `idf`.

More concretely, each item is represented as a set of words, and the items are modeled by TF-IDF weighting of the words. The technique was initially used in information retrieval to model a document-word matrix. In case of our item-word matrices, for a given item ``i``, term frequency (TF) for a term ``t`` is defined as:

```math
\\mathrm{tf}(t, i) = \\frac{n_{t,i}}{N_i},
```

where ``n_{t,i}`` denotes an ``(i, t)`` element in ``I``, and ``N_i`` is the total number of words that an item ``i`` contains. So, ``\\mathrm{tf}(t, i)`` characterizes an item-term pair by normalized frequency of terms. Meanwhile, inverse document frequency (IDF) is computed over ``M`` items as:

```math
\\mathrm{idf}(t) = \\log \\frac{M}{\\mathrm{df}(t)} + 1,
```

where ``\\mathrm{df}(t)`` counts the number of items which associate with a term ``t``. What ``\\mathrm{idf}(t)`` does is to decrease the effect of commonly used terms, because general words cannot characterize a specific item.

Finally, each item-term pair is weighted by ``\\mathrm{tf}(t, i) \\cdot \\mathrm{idf}(t)`` in the TF-IDF scheme. For instance, if we like to recommend web pages to users, we first need to parse sentences on a page and then construct a vector based on the frequency of each term as follows. Similarly, in case that a recommender is running on a folksonomic (i.e. social tagging) service, the frequency of tags corresponds to each dimension of a vector, and ``\\mathrm{tf}(t, i)`` and ``\\mathrm{idf}(t)`` are finally calculated for each item-tag pair.

![tfidf](./assets/images/tfidf.png)
"""
struct TFIDF <: Recommender
    data::DataAccessor # document x attribute
    tf::AbstractMatrix # 1 x attribute
    idf::AbstractMatrix # 1 x attribute

    function TFIDF(data::DataAccessor, tf::AbstractMatrix, idf::AbstractMatrix)
        # instanciate with dummy status (i.e., always true)
        new(data, tf, idf)
    end
end

function predict(recommender::TFIDF, u::Int, i::Int)
    uv = recommender.data.user_attributes[u]
    profile = recommender.data.R' * uv # attribute x 1
    ((recommender.data.R .* recommender.idf) * profile)[i]
end


text:
function NamedTuple(::Missing)
    return (object = missing,)
end

function NamedTuple(d::Dict)
    return (dict_json_string = json(d),)
end
Base.convert(T::Type{NamedTuple}, x::Dict) = T(x)

function Dict(nt::NamedTuple)
    JSON.parse(nt.dict_json_string)
end
# Base.convert(T::Type{Dict}, x::NamedTuple) = T(x) # is now defined in NamedTupleTools-Package


text:
using Functors, Test

@testset "Functors.jl" begin

include("basics.jl")
include("base.jl")

end


text:
using BinDeps

@BinDeps.setup

using Compat

# check for cairo version
function validate_cairo_version(name,handle)
    f = Libdl.dlsym_e(handle, "cairo_version")
    f == C_NULL && return false
    v = ccall(f, Int32,())
    return v > 10800
end

group = library_group("cairo")

deps = [
    libpng = library_dependency("png", aliases = ["libpng","libpng-1.5.14","libpng15","libpng12.so.0","libpng12"], runtime = false, group = group)
    pixman = library_dependency("pixman", aliases = ["libpixman","libpixman-1","libpixman-1-0","libpixman-1.0"], depends = [libpng], runtime = false, group = group)
    libffi = library_dependency("ffi", aliases = ["libffi"], runtime = false, group = group)
    gettext = library_dependency("gettext", aliases = ["libintl", "preloadable_libintl", "libgettextpo"], os = :Unix, group = group)
    gobject = library_dependency("gobject", aliases = ["libgobject-2.0-0", "libgobject-2.0", "libgobject-2_0-0", "libgobject-2.0.so.0"], depends=[libffi, gettext], group = group)
    freetype = library_dependency("freetype", aliases = ["libfreetype"], runtime = false, group = group)
    fontconfig = library_dependency("fontconfig", aliases = ["libfontconfig-1", "libfontconfig", "libfontconfig.so.1"], depends = [freetype], runtime = false, group = group)
    cairo = library_dependency("cairo", aliases = ["libcairo-2", "libcairo","libcairo.so.2", "libcairo2"], depends = [gobject,fontconfig,libpng], group = group, validate = validate_cairo_version)
    pango = library_dependency("pango", aliases = ["libpango-1.0-0", "libpango-1.0","libpango-1.0.so.0", "libpango-1_0-0"], group = group)
    pangocairo = library_dependency("pangocairo", aliases = ["libpangocairo-1.0-0", "libpangocairo-1.0", "libpangocairo-1.0.so.0"], depends = [cairo], group = group)
    zlib = library_dependency("zlib", aliases = ["libzlib","zlib1"], os = :Windows, group = group)
]

if is_windows()
    using WinRPM
    provides(WinRPM.RPM,"libpango-1_0-0",[pango,pangocairo],os = :Windows)
    provides(WinRPM.RPM,["glib2", "libgobject-2_0-0"],gobject,os = :Windows)
    provides(WinRPM.RPM,"zlib-devel",zlib,os = :Windows)
    provides(WinRPM.RPM,["libcairo2","libharfbuzz0"],cairo,os = :Windows)
end

if is_apple()
    if Pkg.installed("Homebrew") === nothing
        error("Homebrew package not installed, please run Pkg.add(\"Homebrew\")")
    end
    using Homebrew
    provides( Homebrew.HB, "cairo", cairo, os = :Darwin )
    provides( Homebrew.HB, "pango", [pango, pangocairo], os = :Darwin, onload =
    """
    function __init__()
        ENV["PANGO_SYSCONFDIR"] = joinpath("$(Homebrew.prefix())", "etc")
    end
    """ )
    provides( Homebrew.HB, "fontconfig", fontconfig, os = :Darwin )
    provides( Homebrew.HB, "glib", gobject, os = :Darwin )
    provides( Homebrew.HB, "libpng", libpng, os = :Darwin )
    provides( Homebrew.HB, "gettext", gettext, os = :Darwin )
    provides( Homebrew.HB, "freetype", freetype, os = :Darwin )
    provides( Homebrew.HB, "libffi", libffi, os = :Darwin )
    provides( Homebrew.HB, "pixman", pixman, os = :Darwin )
end

# System Package Managers
provides(AptGet,
    @compat Dict(
        "libcairo2" => cairo,
        "libfontconfig1" => fontconfig,
        "libpango1.0-0" => [pango,pangocairo],
        "libglib2.0-0" => gobject,
        "libpng12-0" => libpng,
        "libpixman-1-0" => pixman,
        "gettext" => gettext
    ))

# TODO: check whether these are accurate
provides(Yum,
    @compat Dict(
        "cairo" => cairo,
        "fontconfig" => fontconfig,
        "pango" => [pango,pangocairo],
        "glib2" => gobject,
        "libpng" => libpng,
        "gettext-libs" => gettext
    ))

provides(Zypper,
    @compat Dict(
        "libcairo" => cairo,
        "libfontconfig" => fontconfig,
        "libpango-1.0" => [pango,pangocairo],
        "libglib-2.0" => gobject,
        "libpng12" => libpng,
        "libpixman-1" => pixman,
        "gettext" => gettext
    ))

const png_version = "1.5.14"

provides(Sources,
    @compat Dict(
        URI("http://www.cairographics.org/releases/pixman-0.28.2.tar.gz") => pixman,
        URI("http://www.cairographics.org/releases/cairo-1.12.16.tar.xz") => cairo,
        URI("http://download.savannah.gnu.org/releases/freetype/freetype-2.4.11.tar.gz") => freetype,
        URI("http://www.freedesktop.org/software/fontconfig/release/fontconfig-2.10.2.tar.gz") => fontconfig,
        URI("http://ftp.gnu.org/pub/gnu/gettext/gettext-0.18.2.tar.gz") => gettext,
        URI("ftp://ftp.simplesystems.org/pub/libpng/png/src/history/libpng15/libpng-$(png_version).tar.gz") => libpng,
        URI("ftp://sourceware.org/pub/libffi/libffi-3.0.11.tar.gz") => libffi,
        URI("http://ftp.gnome.org/pub/gnome/sources/glib/2.34/glib-2.34.3.tar.xz") => gobject,
        URI("http://ftp.gnome.org/pub/GNOME/sources/pango/1.32/pango-1.32.6.tar.xz") => [pango,pangocairo],
        URI("http://zlib.net/zlib-1.2.7.tar.gz") => zlib
    ))

xx(t...) = (is_windows() ? t[1] : (is_linux() || length(t) == 2) ? t[2] : t[3])

provides(BuildProcess,
    @compat Dict(
        Autotools(libtarget = "pixman/libpixman-1.la", installed_libname = xx("libpixman-1-0.","libpixman-1.","libpixman-1.0.")*Libdl.dlext) => pixman,
        Autotools(libtarget = xx("objs/.libs/libfreetype.la","libfreetype.la")) => freetype,
        Autotools(libtarget = "src/libfontconfig.la") => fontconfig,
        Autotools(libtarget = "src/libcairo.la", configure_options = append!(append!(
                AbstractString[],
                !is_linux() ? AbstractString["--without-x","--disable-xlib","--disable-xcb"] : AbstractString[]),
                is_apple() ? AbstractString["--enable-quartz","--enable-quartz-font","--enable-quartz-image","--disable-gl"] : AbstractString[])) => cairo,
        Autotools(libtarget = "gettext-tools/gnulib-lib/.libs/libgettextlib.la") => gettext,
        Autotools(libtarget = "libffi.la") => libffi,
        Autotools(libtarget = "gobject/libgobject-2.0.la") => gobject,
        Autotools(libtarget = "pango/libpango-1.0.la") => [pango,pangocairo]
    ))

provides(BuildProcess,Autotools(libtarget = "libpng15.la"),libpng,os = :Unix)

provides(SimpleBuild,
    (@build_steps begin
        GetSources(zlib)
        @build_steps begin
            ChangeDirectory(joinpath(BinDeps.depsdir(zlib),"src","zlib-1.2.7"))
            MakeTargets(["-fwin32/Makefile.gcc"])
            #MakeTargets(["-fwin32/Makefile.gcc","DESTDIR=../../usr/","INCLUDE_PATH=include","LIBRARY_PATH=lib","SHARED_MODE=1","install"])
        end
    end),zlib, os = :Windows)

prefix=joinpath(BinDeps.depsdir(libpng),"usr")
uprefix = replace(replace(prefix,"\\","/"),"C:/","/c/")
pngsrcdir = joinpath(BinDeps.depsdir(libpng),"src","libpng-$png_version")
pngbuilddir = joinpath(BinDeps.depsdir(libpng),"builds","libpng-$png_version")
provides(BuildProcess,
    (@build_steps begin
        GetSources(libpng)
        CreateDirectory(pngbuilddir)
        @build_steps begin
            ChangeDirectory(pngbuilddir)
            FileRule(joinpath(prefix,"lib","libpng15.dll"),@build_steps begin
                `cmake -DCMAKE_INSTALL_PREFIX="$prefix" -G"MSYS Makefiles" $pngsrcdir`
                `make`
                `cp 'libpng*.dll' $prefix/lib`
                `cp 'libpng*.a' $prefix/lib`
                `cp 'libpng*.pc' $prefix/lib/pkgconfig`
                `cp pnglibconf.h $prefix/include`
                `cp $pngsrcdir/png.h $prefix/include`
                `cp $pngsrcdir/pngconf.h $prefix/include`
            end)
        end
    end),libpng, os = :Windows)

@BinDeps.install Dict([(:gobject, :_jl_libgobject),
                       (:cairo, :_jl_libcairo),
                       (:pango, :_jl_libpango),
                       (:pangocairo, :_jl_libpangocairo)])


text:
using ConstrainedDynamics

# Parameters
joint_axis = [.0;1.0;1.0]
v1 = -[0.5;0;0.]
v2 = -[0.;0.5;0.]
v3 = -[0.;0;0.5]

length1 = 1.
width, depth = 0.1, 0.1

# Links
origin = Origin{Float64}()
link1 = Box(length1, width, depth, length1, color = RGBA(1., 1., 0.))
link2 = Box(width, length1, depth, length1, color = RGBA(1., 0., 0.))
link3 = Box(width, depth, length1, length1, color = RGBA(1., 0., 1.))

# Constraints
# joint1 = EqualityConstraint(Spherical(origin, link1; p2=v1))
# joint2 = EqualityConstraint(Spherical(origin, link2; p2=v2))
# joint3 = EqualityConstraint(Spherical(origin, link3; p2=v3))
joint1 = EqualityConstraint(Revolute(origin, link1, joint_axis; p2=v1))
joint2 = EqualityConstraint(Revolute(origin, link2, joint_axis; p2=v2))
joint3 = EqualityConstraint(Revolute(origin, link3, joint_axis; p2=v3))

links = [link1;link2;link3]
constraints = [joint1;joint2;joint3]

mech = Mechanism(origin, links, constraints, g = 0.)
setPosition!(origin,link1,Δq = UnitQuaternion(RotZ(0.)),p2=v1)
setPosition!(origin,link2,Δq = UnitQuaternion(RotX(0.)),p2=v2)
setPosition!(origin,link3,Δq = UnitQuaternion(RotX(0.)),p2=v3)


text:
using InfluxDB
using HTTPClient.HTTPC
using Base.Test

# write your own tests here
influx = InfluxDBClient("http://localhost", 8086, "root", "root", "telegraf")
@test influx.baseurl == "http://localhost:8086"
@test influx.database == "telegraf"


text:
module Layout

using Revise
import Genie
using Stipple

export layout

function layout(output::Union{String,Vector}; partial::Bool = false, title::String = "", class::String = "", style::String = "") :: String
  isa(output, Vector) && (output = join(output, '\n'))

  content = string(
    theme(),
    output,
    Stipple.deps()
  )

  partial && return content

  Genie.Renderer.Html.doc(
    Genie.Renderer.Html.html([
      Genie.Renderer.Html.head([
        Genie.Renderer.Html.title(title)
        Genie.Renderer.Html.meta(name="viewport", content="width=device-width, initial-scale=1, maximum-scale=1, user-scalable=no, minimal-ui")
      ])
      Genie.Renderer.Html.body(content, class=class, style=style)
    ])
  )
end

include(joinpath("layout", "page.jl"))
include(joinpath("layout", "theme.jl"))

end


text:
struct ConstantManifoldSelector{Tm} <: AbstractUpdateSelector
    manifold::Tm
end

function select_update!(constantsel::ConstantManifoldSelector, state, optimizer, pb)
    state.previous_update = state.selected_update
    state.selected_update = optimizer.manifold_update

    state.M = constantsel.manifold
    return
end

Base.summary(sel::ConstantManifoldSelector) = string("Cst ", sel.manifold)


text:
using SharedMATLABEngine
using Test

@testset "SharedMATLABEngine.jl" begin
    # Write your tests here.
end


text:
abstract type EoS end
abstract type SAFT <: EoS end

abstract type PCSAFTFamily <: SAFT end
abstract type sPCSAFTFamily <: SAFT end
abstract type ogSAFTFamily <: SAFT end
abstract type SAFTVRMieFamily <: SAFT end

struct SAFTVRMie <: SAFTVRMieFamily; components; params::SAFTVRMieParams end
struct PCSAFT <: PCSAFTFamily; components; params::PCSAFTParams end
struct sPCSAFT <: sPCSAFTFamily; components; params::sPCSAFTParams end
struct ogSAFT <: ogSAFTFamily; components; params::ogSAFTParams end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: finspace
using AWS.Compat
using AWS.UUIDs

"""
    create_environment(name)
    create_environment(name, params::Dict{String,<:Any})

Create a new FinSpace environment.

# Arguments
- `name`: The name of the FinSpace environment to be created.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"description"`: The description of the FinSpace environment to be created.
- `"federationMode"`: Authentication mode for the environment.    FEDERATED - Users access
  FinSpace through Single Sign On (SSO) via your Identity provider.    LOCAL - Users access
  FinSpace via email and password managed within the FinSpace environment.
- `"federationParameters"`: Configuration information when authentication mode is FEDERATED.
- `"kmsKeyId"`: The KMS key id to encrypt your data in the FinSpace environment.
- `"tags"`: Add tags to your FinSpace environment.
"""
create_environment(name; aws_config::AbstractAWSConfig=global_aws_config()) = finspace("POST", "/environment", Dict{String, Any}("name"=>name); aws_config=aws_config)
create_environment(name, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = finspace("POST", "/environment", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("name"=>name), params)); aws_config=aws_config)

"""
    delete_environment(environment_id)
    delete_environment(environment_id, params::Dict{String,<:Any})

Delete an FinSpace environment.

# Arguments
- `environment_id`: The identifier for the FinSpace environment.

"""
delete_environment(environmentId; aws_config::AbstractAWSConfig=global_aws_config()) = finspace("DELETE", "/environment/$(environmentId)"; aws_config=aws_config)
delete_environment(environmentId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = finspace("DELETE", "/environment/$(environmentId)", params; aws_config=aws_config)

"""
    get_environment(environment_id)
    get_environment(environment_id, params::Dict{String,<:Any})

Returns the FinSpace environment object.

# Arguments
- `environment_id`: The identifier of the FinSpace environment.

"""
get_environment(environmentId; aws_config::AbstractAWSConfig=global_aws_config()) = finspace("GET", "/environment/$(environmentId)"; aws_config=aws_config)
get_environment(environmentId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = finspace("GET", "/environment/$(environmentId)", params; aws_config=aws_config)

"""
    list_environments()
    list_environments(params::Dict{String,<:Any})

A list of all of your FinSpace environments.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to return in this request.
- `"nextToken"`: A token generated by FinSpace that specifies where to continue pagination
  if a previous request was truncated. To get the next set of pages, pass in the nextToken
  value from the response object of the previous page call.
"""
list_environments(; aws_config::AbstractAWSConfig=global_aws_config()) = finspace("GET", "/environment"; aws_config=aws_config)
list_environments(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = finspace("GET", "/environment", params; aws_config=aws_config)

"""
    list_tags_for_resource(resource_arn)
    list_tags_for_resource(resource_arn, params::Dict{String,<:Any})

A list of all tags for a resource.

# Arguments
- `resource_arn`: The Amazon Resource Name of the resource.

"""
list_tags_for_resource(resourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = finspace("GET", "/tags/$(resourceArn)"; aws_config=aws_config)
list_tags_for_resource(resourceArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = finspace("GET", "/tags/$(resourceArn)", params; aws_config=aws_config)

"""
    tag_resource(resource_arn, tags)
    tag_resource(resource_arn, tags, params::Dict{String,<:Any})

Adds metadata tags to a FinSpace resource.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) for the resource.
- `tags`: One or more tags to be assigned to the resource.

"""
tag_resource(resourceArn, tags; aws_config::AbstractAWSConfig=global_aws_config()) = finspace("POST", "/tags/$(resourceArn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resourceArn, tags, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = finspace("POST", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), params)); aws_config=aws_config)

"""
    untag_resource(resource_arn, tag_keys)
    untag_resource(resource_arn, tag_keys, params::Dict{String,<:Any})

Removes metadata tags from a FinSpace resource.

# Arguments
- `resource_arn`: A FinSpace resource from which you want to remove a tag or tags. The
  value for this parameter is an Amazon Resource Name (ARN).
- `tag_keys`: The tag keys (names) of one or more tags to be removed.

"""
untag_resource(resourceArn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()) = finspace("DELETE", "/tags/$(resourceArn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceArn, tagKeys, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = finspace("DELETE", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), params)); aws_config=aws_config)

"""
    update_environment(environment_id)
    update_environment(environment_id, params::Dict{String,<:Any})

Update your FinSpace environment.

# Arguments
- `environment_id`: The identifier of the FinSpace environment.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"description"`: The description of the environment.
- `"federationMode"`: Authentication mode for the environment.    FEDERATED - Users access
  FinSpace through Single Sign On (SSO) via your Identity provider.    LOCAL - Users access
  FinSpace via email and password managed within the FinSpace environment.
- `"federationParameters"`:
- `"name"`: The name of the environment.
"""
update_environment(environmentId; aws_config::AbstractAWSConfig=global_aws_config()) = finspace("PUT", "/environment/$(environmentId)"; aws_config=aws_config)
update_environment(environmentId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = finspace("PUT", "/environment/$(environmentId)", params; aws_config=aws_config)


text:
# This file was generated, do not modify it. # hide
y11 = predict(m11, restrict(X, f, 1));
y12 = predict(m12, restrict(X, f, 2));
y13 = predict(m13, restrict(X, f, 3));


text:
type ProdMean <: Mean
    means::Vector{Mean}
    function ProdMean(args...)
        means = Array(Mean, length(args))
        for (i,m) in enumerate(args)
            isa(m, Mean) || throw(ArgumentError("All arguments of ProdMean must be Mean objects"))
            means[i] = m
        end
        return new(means)
    end
end

function show(io::IO, pm::ProdMean, depth::Int = 0)
    pad = repeat(" ", 2 * depth)
    println(io, "$(pad)Type: $(typeof(pm))")
    for m in pm.means
        show(io, m, depth+1)
    end
end

function mean(prodmean::ProdMean, x::Matrix{Float64})
    p = 1.0
    for m in prodmean.means
        p = p.*mean(m, x)
    end
    return p
end

function get_params(prodmean::ProdMean)
    p = Array(Float64, 0)
    for m in prodmean.means
        append!(p, get_params(m))
    end
    p
end

get_param_names(prodmean::ProdMean) = composite_param_names(prodmean.means, :pm)

function num_params(prodmean::ProdMean)
    n = 0
    for m in prodmean.means
        n += num_params(m)
    end
    n
end

function set_params!(prodmean::ProdMean, hyp::Vector{Float64})
    i, n = 1, num_params(prodmean)
    length(hyp) == num_params(prodmean) || throw(ArgumentError("ProdMean object requires $(n) hyperparameters"))
    for m in prodmean.means
        np = num_params(m)
        set_params!(m, hyp[i:(i+np-1)])
        i += np
    end
end

function grad_mean(prodmean::ProdMean, x::Vector{Float64})
     dm = Array(Float64, 0)
      for m in prodmean.means
          p = 1.0
          for j in prodkern.means[find(k.!=prodkern.means)]
              p = p.*mean(j, x)
          end
        append!(dm,grad_mean(m, x).*p)
      end
    dm
end

# Multiplication operators
function *(m1::ProdMean, m2::Mean)
    means = [m1.means, m2]
    ProdMean(means...)
end
function *(m1::ProdMean, m2::ProdMean)
    means = [m1.means, m2.means]
    ProdMean(means...)
end
*(m1::Mean, m2::Mean) = ProdMean(m1,m2)
*(m1::Mean, m2::ProdMean) = *(m2,m1)


text:
function KS_solve_ChebySCF( pw::PWGrid,
                            V_ionic, Focc, Nstates::Int64;
                            cheby_degree=5,
                            β = 0.5, E_NN = nothing,
                            rho0 = nothing,
                            Potentials0 = nothing,
                            NiterMax=100,
                            verbose=false )

    Ngwx = pw.gvecw.Ngwx
    Ns = pw.Ns
    Npoints = prod(Ns)
    ΔV = pw.Ω/Npoints

    # starting rho
    if rho0 == nothing
        srand(1234)
        v = rand( ComplexF64, Ngwx, Nstates )
        v = ortho_gram_schmidt( v )
        rho = calc_rho( pw, Focc, v )
    else
        srand(1234)
        v = rand( ComplexF64, Ngwx, Nstates ) # What if v0 is provided ?
        rho = copy( rho )
    end

    # starting Potentials
    if Potentials0 == nothing
        # initialize potential: Ionic, Hartree, XC
        V_Hartree = real( G_to_R( Ns, Poisson_solve(pw, rho) ) )
        V_xc = excVWN( rho ) + rho .* excpVWN( rho )
        Potentials = PotentialsT( V_ionic, V_Hartree, V_xc )
    else
        Potentials = PotentialsT( Potentials0.Ionic,
                                  Potentials0.Hartree,
                                  Potentials0.XC )
    end

    Etot_old = 0.0
    rho_old = copy(rho)

    if( E_NN != nothing )
        Energies = EnergiesT( 0.0, 0.0, 0.0, 0.0, 0.0, E_NN )
    else
        Energies = EnergiesT( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 )
    end

    # eigenvalues
    λ = zeros(Nstates)
    resnrm = zeros(Nstates)

    for iter = 1:NiterMax

        ub, lb = get_ub_lb_lanczos(pw, Potentials, Nstates*2)
        #@printf("lb = %f, ub = %f\n", lb, ub)

        v = chebyfilt(pw, Potentials, v, cheby_degree, lb, ub)
        v = ortho_gram_schmidt(v)

        #
        rho_new = calc_rho( pw, Focc, v )
        diffRho = norm(rho_new - rho)
        #
        rho = β*rho_new[:] + (1-β)*rho[:]
        #
        #integRho = sum(rho)*ΔV
        #@printf("integRho = %18.10f\n", integRho)
        #
        V_Hartree = real( G_to_R( Ns, Poisson_solve(pw, rho) ) )
        V_xc = excVWN( rho ) + rho .* excpVWN( rho )
        Potentials = PotentialsT( V_ionic, V_Hartree, V_xc )

        # Calculate energies and update potetials
        Energies = calc_Energies( pw, Potentials, Focc, v, Energies.NN )

        Etot = Energies.Total
        diffE = abs( Etot - Etot_old )

        #
        @printf("ChebySCF: %8d %18.10f %18.10e %18.10e\n", iter, Etot, diffE, diffRho )

        if diffE < 1e-7
            @printf("ChebySCF is is converged: iter: %d , diffE = %10.7e\n", iter, diffE)
            break
        end
        #
        Etot_old = Etot
    end

    # Calculate eigenvalues
    Hv = op_H(pw,Potentials,v)
    G = v' * Hv
    R = Hv - v*G
    λ = real(eigvals(G))
    #for j = 1:Nstates
    #   resnrm[j] = norm(R[:,j]);
    #   @printf("λ[%2d] = %11.3e, resnrm = %11.3e\n", j, λ[j], resnrm[j]);
    #end

    return Energies, Potentials, v, λ

end


text:
module Declares

    # Импортируем бибилиотеку работы со временем
    include("Chronometer.jl")
    import .Chronometer

    export queues, exchangers, in_channel

    """ Функция регистрации очередей по именам.

    channel: канал связи с брокером
    queues_names: список имён очередей
    Adapter: адаптер брокера очередей
    """
    function queues(channel, queues_names, Adapter)
        @info Chronometer.message_with_time("Выполняется функция регистрации очередей...")

        @debug "Параметры вызова функции:"
        @debug "> канал связи с брокером" channel
        @debug "> список имён очередей" queues_names
        @debug "> адаптер брокера очередей" Adapter

        queues = nothing
        @debug "Инициализируем пустой список очередей" queues

        @debug "Перебираем имена очередей и регистрируем каждую"
        for name in queues_names
            @debug "> Имя для регистрации очереди" name

            queue = Adapter.Queues.declare(channel, name)
            @debug "> Зарегистрированная очередь" queue

            if queue !== nothing
                if queues === nothing
                    queues = Dict{String, Any}(name => queue)
                    @debug ">> Создали словарь зарегистрированных очередей" queues
                else
                    push!(queues, name => queue)
                    @debug ">> Обновили словарь зарегистрированных очередей" queues
                end
            end
        end

        @debug "Список зарегистрированных очередей" queues
        return queues
    end

    """ Функция регистрации обменников.

    channel: канал связи с брокером
    exchanger_config: конфигурация обменников
    Adapter: адаптер брокера очередей
    """
    function exchangers(channel, exchangers_configs, Adapter)
        @info Chronometer.message_with_time("Выполняется функция регистрации обменников...")

        @debug "Параметры вызова функции:"
        @debug "- канал связи с брокером" channel
        @debug "- конфигурации обменников канала" exchangers_configs
        @debug "- адаптер брокера очередей" Adapter

        exchangers = nothing
        @debug "Инициализируем пустое значение словаря обменников" exchangers

        name_field::String = "name"
        @debug "Имя поля наименования обменника" name_field

        type_field::String = "type"
        @debug "Имя поля типа обменника" type_field

        @debug "Перебираем список настроек объетов конфигурации обменников и регистрируем их"
        for exchanger_config in exchangers_configs
            @debug "> Конфигурация обмненика" exchanger_config

            name = get(exchanger_config, name_field, nothing)
            @debug "> Имя обменника из конфигурации" name

            type = get(exchanger_config, type_field, "direct")
            @debug "> Тим обменника из конфигурации" type

            if name !== nothing
                exchanger::Dict{String, Any} = Adapter.Exchanges.declare(channel, name, type)
                @debug ">> Зарегистрированный обменник" exchanger

                if exchangers === nothing
                    exchangers = Dict{String, Any}(name => exchanger)
                    @debug ">>> Создали словарь зарегистрированных обменников" exchangers
                else
                    push!(exchangers, name => exchanger)
                    @debug ">>> Обновоили словарь зарегистрированных обменников" exchangers
                end
            end
        end

        @debug "Зарегистрированные обменники" exchangers
        return exchangers
    end

    """ Функця регистрации очередей и обмеников в канале

    channel: канал связи с брокером
    channel_config: конфигурация канала
    Adapter: адаптер брокера
    """
    function in_channel(channel, channel_config, Adapter)
        @info Chronometer.message_with_time("Выполняется функция регистрации очередей, обменников и т.п. в канале...")

        @debug "Параметры вызова функции:"
        @debug "- канал связи с брокером" channel
        @debug "- конфигурация канала" channel_config
        @debug "- адаптер брокера" Adapter

        exchangers_dict = nothing
        @debug "Задаём пустое значение для зарегистрированных обменников" exchangers_dict

        queues_dict = nothing
        @debug "Задаём пустое значене для зарегистрированных очередей" queues_dict

        queues_field::String = "queues"
        exchangers_field::String = "exchangers"
        @debug "Имя поля конфигурации очередей" queues_field
        @debug "Имя поля конфигурации обменников" exchangers_field

        @debug "Проверяем наличие конфигурации канала и самого канала, если все есть - регистрируем очереди и обменники"
        if channel_config !== nothing && channel !== nothing
            exchangers_list = get(channel_config, exchangers_field, nothing)
            @debug "> Список настроек для обменников из конфигурации" exchangers_list

            queues_list = get(channel_config, queues_field, nothing)
            @debug "> Список настроек для обнменников из конфигурации" queues_list

            exchangers_dict = exchangers(channel, exchangers_list, Adapter)
            @debug "> Словарь зарегистрированных обменников" exchangers_dict

            queues_dict = queues(channel, queues_list, Adapter)
            @debug "> Словарь зарегистрированных очередей" queues_dict
        end

        @debug "Зарегистрированные обменники" exchangers_dict
        @debug "Зарегистрированные очереди" queues_dict
        return exchangers_dict, queues_dict
    end
end


text:
"""
find_kbest_assignments(C, k, maximize=false) -> kbest_sols

Find the k lowest [or highest] cost 2D assignments for the two-dimensional
assignment problem with a rectangular cost matrix C.

# Example
```julia
julia> M=rand(1:100,3,4)
3×4 Matrix{Int64}:
 77  51  42  67
 72  53  47   4
 24  50  77  96

 julia> sols = find_kbest_assigments(M, 5)
 5-element Vector{Assignment.AssignmentSolution{Int64, Int64}}:
  AssignmentSolution(CartesianIndex.(1:3, [3, 4, 1]), 70)
  AssignmentSolution(CartesianIndex.(1:3, [2, 4, 1]), 79)
  AssignmentSolution(CartesianIndex.(1:3, [3, 4, 2]), 96)
  AssignmentSolution(CartesianIndex.(1:3, [3, 2, 1]), 119)
  AssignmentSolution(CartesianIndex.(1:3, [2, 3, 1]), 122)

 julia> max_sols = find_kbest_assignments(M, 5, true)
 5-element Vector{Assignment.AssignmentSolution{Int64, Int64}}:
  AssignmentSolution(CartesianIndex.(1:3, [1, 2, 4]), 226)
  AssignmentSolution(CartesianIndex.(1:3, [1, 3, 4]), 220)
  AssignmentSolution(CartesianIndex.(1:3, [2, 1, 4]), 219)
  AssignmentSolution(CartesianIndex.(1:3, [4, 1, 3]), 216)
  AssignmentSolution(CartesianIndex.(1:3, [3, 1, 4]), 210)
```

This is an implementation of Murty's method, which is described in [1].
The algorithm relies on the existence of a 2D assignment algorithm. The 2D
assignment algorithm of [2] is used. Additionally, the dual variable
inheritance methods described in [3] is used to reduce the computational
complexity of the technique.

Murty's algorithm runs 2D assignment algorithms a number of times with an
increasing number of constraints. Instances of MurtyData are stored in an
ordered list implemented using the BinaryMinHeap class.

This code is a port of Matlab code released in the public domain by the
US Naval Research Laboratory [4].

This work is not affliated with or endorsed by the US Naval Research Laboratory.

# References
[1] K. G. Murty, "An algorithm for ranking all the assignments in order of
   increasing cost;" Operations Research; vol. 16; no. 3; pp. 682-687
   May-Jun. 1968.
[2] D. F. Crouse, "On Implementing 2D Rectangular Assignment Algorithms,"
   IEEE Transactions on Aerospace & Electronic Systems; vol. 52; no. 4
   pp. 1679-1696; Aug. 2016.
[3] M. L. Miller, H. S. Stone, & J. Cox, Ingemar, "Optimizing Murty's
   ranked assignment method;" IEEE Transactions on Aerospace &
   Electronic Systems; vol. 33; no. 3; pp. 851-862; Jul. 1997.
[4] D. F. Crouse, "The Tracker Component Library: Free Routines for Rapid
   Prototyping," IEEE Aerospace and Electronic Systems Magazine, vol. 32,
   no. 5, pp. 18-27, May. 2017

October 2013 David F. Crouse; Naval Research Laboratory; Washington D.C.
(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.
"""
function find_kbest_assignments(C::AbstractMatrix{T}, k::Integer, maximize=false) where T
    numRow, numCol = size(C)

	if numRow > numCol
		sols = find_kbest_assignments(C', k, maximize)
		for (i, sol) in enumerate(sols)
			sols[i] = sol'
		end
		return sols
	end

	#The cost matrix must have all non-negative elements for the assignment
	#algorithm to work. This forces all of the elements to be positive. The
	#delta is added back in when computing the gain in the end.
    if maximize
        CDelta = maximum(C)
        C = -C.+ CDelta
	else
        CDelta = minimum(C)
        C =  C .- CDelta
    end

	#Step 1: Find the best assignment. It will become the root from which()
	#everything else splits.
    #MUST PAD IT TO MAKE IT SQUARE; because we are using suboptimal
    #dual variables.
	numPad = numCol-numRow
	numPad > 0 && (C = [C; zeros(T, numPad, numCol)])

    LCHyp = MurtyData(C, numRow)

	#Check for feasibility.
    if isnothing(LCHyp.gainFull)
        return AssignmentSolution{T, _dual_eltype(T)}[]
    end

	k_best = [AssignmentSolution(LCHyp, CDelta, maximize)]
	sizehint!(k_best, k)

	hyp_list = BinaryMinHeap([LCHyp])

	ScannedRows = falses(numCol)
	ScannedCol = falses(numCol)
	pred = zeros(Int, numCol)
	shortestPathCost = fill(typemax(float(T)), numCol)

    for curSweep in 2:k
        #We have to successively split the LC hypothesis either k times |
        #until splitting is no longer possible & get rid of the
        #hypothesis from the list.
		data = pop!(hyp_list)

		#Basically; it is all columns that have not been assigned in
		#one of the subproblems.
		col2Scan = @view data.col4rowLCFull[data.activeRow:end]

		#The loop ends at numVarRow; because we are not going to allow
		#variable hypotheses to be applied to dummy rows.
		for curRow=data.activeRow:data.numVarRow
			#forbiddenColumns contains a 1 for every column that one is
			#not allowed to assign to curRow.
			if curRow==data.activeRow
				forbiddenColumns = data.forbiddenActiveCol
			else
				forbiddenColumns = falses(numCol)
				forbiddenColumns[data.col4rowLCFull[curRow]] = true
			end

			#Get rid of the current assignment.
			row4colInit = data.row4colLCFull[:]
			col4rowInit = data.col4rowLCFull[:]
			row4colInit[col4rowInit[curRow]] = 0
			col4rowInit[curRow] = 0

			splitHyp = ConstrainedShortestPath!(
				data.A, data.numVarRow, curRow,
				forbiddenColumns, col4rowInit, row4colInit, col2Scan[:], data.u[:], data.v[:],
				ScannedRows, ScannedCol, pred, shortestPathCost
			)

			#If a valid assignment was produced.
			!isnothing(splitHyp.gainFull) && push!(hyp_list, splitHyp)

			#Remove the current assigned column from the list of
			#columns that can be scanned.
			col2Scan = @view col2Scan[col2Scan .!= data.col4rowLCFull[curRow]]
		end

        if !isempty(hyp_list)
			#Save the ordered best solutions.
			push!(k_best, AssignmentSolution(first(hyp_list), CDelta, maximize))
        else
            break #All possible hypotheses have been enumerated.
        end
    end

	return k_best
end

struct MurtyData{T, M<:AbstractMatrix{T}, S}
	A::M
	numVarRow::Int
	activeRow::Int
	forbiddenActiveCol::BitVector
	col4rowLCFull::Vector{Int}
	row4colLCFull::Vector{Int}
	gainFull::Union{T, Nothing}
	u::Vector{S}
	v::Vector{S}
end

function MurtyData(A::AbstractMatrix, numVarRow::Integer)
	numRow, numCol = size(A)

	sol = find_best_assignment(A')
	forbiddenActiveCol = falses(numCol)
	if !isempty(sol.col4row)
		forbiddenActiveCol[sol.col4row[1]] = 1
	end

	return MurtyData(
		A, numVarRow, 1, forbiddenActiveCol,
		sol.col4row, sol.row4col, sol.cost, sol.v, sol.u)
end

Base.isless(a::MurtyData, b::MurtyData) = Base.isless(a.gainFull, b.gainFull)

for op in (:<, :>, :(<=), :(>=), :(==), :(!=))
	@eval begin
		$op(a::MurtyData, B::MurtyData) = $op(a.gainFull, b.gainFull)
		$op(a::MurtyData, B) = $op(a.gainFull, b)
		$op(a, B::MurtyData) = $op(a, b.gainFull)
	end
end

function AssignmentSolution(data::MurtyData{T, M, S}, CDelta=T(0), maximize=false) where {T, M, S}
	numCol = size(data.A, 2)
	numRow = data.numVarRow

	col4row = data.col4rowLCFull[1:numRow]
	row4col = data.row4colLCFull[:]
	if numRow < numCol
		for (i, j) in enumerate(row4col)
			if j > numRow
				row4col[i] = 0
			end
		end
	end
	gain =  maximize ? CDelta * numRow  - data.gainFull : CDelta * numRow + data.gainFull

	return AssignmentSolution{T, S}(col4row, row4col, T(gain), data.u[:], data.v[1:numRow])
end

function ConstrainedShortestPath!(
		C::AbstractMatrix{T},
		numVarRow::Integer,
		activeRow::Integer,
		forbiddenActiveCols::AbstractVector{Bool},
		col4row::AbstractVector{Int},
		row4col::AbstractVector{Int},
		col2Scan::AbstractVector{Int},
		u::AbstractVector{S},
		v::AbstractVector{S},
		ScannedRows::AbstractVector{Bool} = falses(size(C,1)),
		ScannedCol::AbstractVector{Bool} = falses(size(C,2)),
    	pred::AbstractVector{Int} = zeros(Int, size(C,2)),
		shortestPathCost::AbstractVector = fill(Inf, size(C,2))
	) where {T, S}

	numRow = size(C,1)
	numCol = size(C,2)
    numCol2Scan=length(col2Scan)
	ScannedRows[:] .= false
	ScannedCol[:] .= false
	pred[:] .= 0
	shortestPathCost[:] .= Inf

	sink = 0
	delta = S(0)
    curRow = activeRow
	closestColScan=0

	while sink==0
		ScannedRows[curRow]=1;

        # Scan all of the columns that have not already been scanned.
        minVal=Inf;
        @inbounds for curColScan=1:numCol2Scan
            curCol=col2Scan[curColScan];

			# Do not scan forbidden columns.
            curRow==activeRow && forbiddenActiveCols[curCol] && continue

            reducedCost=delta+C[curRow,curCol]-u[curRow]-v[curCol]
            if reducedCost < shortestPathCost[curCol]
                pred[curCol] = curRow
                shortestPathCost[curCol] = reducedCost
            end

            # Find the minimum unassigned column that was
            # scanned.
            if(shortestPathCost[curCol]<minVal)
                minVal=shortestPathCost[curCol];
                closestColScan=curColScan;
            end
        end

		# If the minimum cost column is not finite, then the problem is
		# not feasible.
		!isfinite(minVal) && return MurtyData(
			C, numVarRow, activeRow, forbiddenActiveCols,
			Int[], Int[], nothing, u, v)

        closestCol=col2Scan[closestColScan];

        # Add the column to the list of scanned columns and delete it from
        # the list of columns to scan.
        ScannedCol[closestCol] = 1
        numCol2Scan -= 1
		@inbounds for i in closestColScan:numCol2Scan
			col2Scan[i] = col2Scan[i + 1]
		end

        delta = S(shortestPathCost[closestCol])

        # If we have reached an unassigned row.
        if row4col[closestCol]==0
            sink=closestCol
        else
            curRow=row4col[closestCol]
        end
    end

	# Dual Update Step
    # Update the first row in the augmenting path.
    u[activeRow] = u[activeRow] + delta;
    # Update the rest of the rows in the agumenting path.
	@inbounds for (i, col) in enumerate(ScannedRows)
		if col != 0 && i != activeRow
			u[i] += delta - S(shortestPathCost[col4row[i]])
		end
	end
    # Update the scanned columns in the augmenting path.
	@inbounds for (i, row) in enumerate(ScannedCol)
		if row != 0
			v[i] -= delta - S(shortestPathCost[i])
		end
	end

	# Augmentation
    # We have to remove node k from those that must be assigned.
    j = sink
    while true
        i = row4col[j] = pred[j]
		j, col4row[i] = col4row[i], j
		i==activeRow && break
    end

    # Calculate the gain that should be returned. Because the duality gap is
    # zero, the gain can be most easily obtained from the dual cost
    # function.
    gain=T(0)
	@inbounds for (i, j) in enumerate(col4row)
        gain=gain+C[i, j];
    end

	forbiddenActiveCols[col4row[activeRow]] = true

	return MurtyData(
		C, numVarRow, activeRow, forbiddenActiveCols,
		col4row, row4col, gain, u, v)
end


text:
function grammar_koza(; erc=false)
    Base.eval(Main, :(protectedDiv(x, y) = Benchmarks.protectedDiv(x, y)))
    Base.eval(Main, :(protectedSin(x) = Benchmarks.protectedSin(x)))
    Base.eval(Main, :(protectedCos(x) = Benchmarks.protectedCos(x)))
    grammar = erc ?
        (@grammar begin
            R = x
            R = R + R
            R = R - R
            R = R * R
            R = protectedDiv(R,R)
            R = protectedSin(R)
            R = protectedCos(R)
            R = exp(R)
            R = log(abs(R))
            R = _(2*rand()-1)
        end) :
        (@grammar begin
            R = x
            R = R + R
            R = R - R
            R = R * R
            R = protectedDiv(R,R)
            R = protectedSin(R)
            R = protectedCos(R)
            R = exp(R)
            R = log(abs(R))
            R = |(-1.0:0.25:1.0)
        end)
    grammar
end

gt_koza_1(x) = x^4 + x^3 + x^2 + x
gt_koza_2(x) = x^5 - 2x^3 + x
gt_koza_3(x) = x^6 - 2x^4 + x^2

function loss_koza(gt::Function, S::SymbolTable, tree::RuleNode, grammar::Grammar)
    ex = get_executable(tree, grammar)
    los = 0.0
    n = 0
    for x in range(-1.0,stop=1.0,length=20)
        S[:x] = x
        los += abs2(Core.eval(S,ex) - gt(x))
        n += 1
    end
    los / n
end

function run_koza_1(seed::Int=0, n_pop::Int=300)
    Random.seed!(seed)

    grammar = grammar_koza(erc=true)
    S = SymbolTable(grammar, ExprOptimization.Benchmarks)
    loss(tree::RuleNode, grammar::Grammar) = loss_koza(gt_koza_1, S, tree, grammar)

    init_method = RandomInit()
    select_method = TournamentSelection(3)
    p = GeneticProgram(n_pop, #pop_size
                       40,  #iterations
                       10,  #max_depth
                       0.4, #p_reproduction
                       0.5, #p_crossover
                       0.1, #p_mutation
                       init_method=init_method,
                       select_method=select_method)
    tstart = CPUtime_us()
    result = optimize(p, grammar, :R, loss)
    dt = CPUtime_us() - tstart
    result, dt
end


text:
# Note this file is outdates
using Test, Printf, Random

# TODO What is this?
#@everywhere include(joinpath(dirname(@__FILE__),"..","src","jl"))

@everywhere using Polynomials, IterTools, Plots

function plotall(Λ,ℓ,t,g, both=true)
    M = size(Λ,1)
    N = length(t)
    γ = both ? 2 : 1
    p = plot(layout=(M*γ,N-1), size=(1400,1000))
    k = 1
    for m = M:-1:1
        for n = 1:N-1
            idx = (n-1)*M+m #Linear index in Λ
            if isassigned(Λ,idx)
                x, y, x_all, y_all = get_vals(Λ[m,n])
                for i in eachindex(x_all)
                    plot!(p, y_all[i][[1,end]], x_all[i][[1,end]], l = 0, m=(1,:cross,:orange), subplot=k, lab="")
                    plot!(p, y_all[i], x_all[i], l=(1,:blue), lab="", subplot=k)
                end
                I, yI, fI = recover_solution(Λ[m,n], n, N)
                ystr = @sprintf("%3.2f", yI)
                fstr = @sprintf("%3.2f", fI)
                plot!(p, y, x, l=:red, subplot=k, lab="$I, $ystr, $fstr", ylims=(-4,4), xlims=(0,6))
                if both && length(I) > 1 && fI != NaN
                    Y, _ = find_optimal_y_values(ℓ, I)
                    tcont = range(t[1], stop=t[end], length=100*length(t))
                    plot!(p, tcont/t[end]*(N-1)+1, g.(tcont), l=:blue, lab="", subplot=k+N-1)
                    plot!(p, I, Y, m=:circle, l=:red, lab="", subplot=k+N-1)
                end
            end
            k += 1
        end
        k += both ? N-1 : 0
    end
    return p
end

##
Random.seed!(9)
N = 20
ran = cumsum(0.1*randn(10*N))

t = range(0, stop=1, length=N)

g6(t) = ran[floor(Int64, t*10(N-1)+1)]
ℓ = compute_transition_costs(g6, t);

Λ_0 = [create_new_pwq(minimize_wrt_x2(ℓ[i, N], QuadraticPolynomial{Float64}(0.,0.,0.))) for i in 1:N-1];
Λ, t2, _, _, _ = @timed pwq_dp_constrained(Λ_0, ℓ, 19);
##
#[(isassigned(Λ,i) ? length(Λ[i]) : 0) for i in eachindex(Λ)]
lengths = reshape([(isassigned(Λ,i) ? length(Λ[i]) : 0) for i in eachindex(Λ)], size(Λ))
plotly()
heatmap(lengths, show=true)
# p = plotall(Λ,ℓ,t,g, true);
# @time I2, y2, f2 = recover_solution(Λ[19, 1], 1, N)
# Y2, _ = find_optimal_y_values(ℓ, I2)
# plot(range(1, stop=N, length=N), g5.(range(1, stop=N, length=N)))
# plot!(I2, Y2)
# plot!(p, show=true)


text:
module Wflow

using Dates
using TOML
using LightGraphs
using NCDatasets
using StaticArrays
using Statistics
using UnPack
using CSV
using Random
using DataFrames

mutable struct Clock
    time::DateTime
    iteration::Int
    Δt::Second
end

include("io.jl")

"""
    Model{N,L,V,R,W}

Composite type that represents all different aspects of a Wflow Model, such as the
network, parameters, clock, configuration and input and output.
"""
struct Model{N,L,V,R,W}
    config::Config  # all configuration options
    network::N  # connectivity information, directed graph
    lateral::L  # lateral model that holds lateral state, moves along network
    vertical::V  # vertical model that holds vertical state, independent of each other
    clock::Clock  # to keep track of simulation time
    reader::R  # provides the model with dynamic input
    writer::W  # writes model output
end

# prevent a large printout of model components and arrays
Base.show(io::IO, m::Model) = print(io, "model of type ", typeof(m))

include("horizontal_process.jl")
include("hbv.jl")
include("sbm.jl")
include("reservoir_lake.jl")
include("hbv_model.jl")
include("sbm_model.jl")
include("flow.jl")
include("vertical_process.jl")
include("groundwater/connectivity.jl")
include("groundwater/aquifer.jl")
include("groundwater/boundary_conditions.jl")
include("sbm_gwf_model.jl")
include("utils.jl")

"""
    run_simulation(tomlpath::String)
    run_simulation(config::Config)
    run_simulation(model::Model)

Run an entire simulation starting either from a path to a TOML settings file,
a prepared `Config` object, or an initialized `Model` object. This allows more flexibility
if you want to for example modify a `Config` before initializing the `Model`.
"""
function run_simulation(tomlpath)
    config = Config(tomlpath)
    run_simulation(config)
end

function run_simulation(config::Config)
    modeltype = config.model.type

    model = if modeltype == "sbm"
        initialize_sbm_model(config)
    elseif modeltype == "sbm_gwf"
        initialize_sbm_gwf_model(config)
    elseif modeltype == "hbv"
        initialize_hbv_model(config)
    else
        error("unknown model type")
    end

    run_simulation(model)
end

function run_simulation(model::Model; close_files = true)
    @unpack network, config, writer, clock = model

    # in the case of sbm_gwf it's currently a bit hard to use dispatch
    update_func = config.model.type == "sbm_gwf" ? update_sbm_gwf : update

    while true
        model = update_func(model)
        is_finished(clock, config) && break
    end

    # write output state NetCDF
    # undo the clock advance at the end of the last iteration, since there won't
    # be a next step, and then the output state falls on the correct time
    rewind!(clock)
    write_netcdf_timestep(model, writer.state_dataset, writer.state_parameters)

    reset_clock!(model.clock, config)

    # option to support running function twice without re-initializing
    # and thus opening the NetCDF files
    if close_files
        Wflow.close_files(model, delete_output = false)
    end
    return model
end

end # module


text:
function sv_next!(Car::CarSt, a::Float64, Δt::Float64)
    Car.s += Car.v*Δt + Δt^2*a/2
    Car.v += Δt*a
    return nothing
end

function sv_next(Car::CarSt, a::Float64, Δt::Float64)
    Car_next = deepcopy(Car)
    sv_next!(Car_next, a, Δt)
    return Car_next
end

########## transition ########
# generate the probability of ending in state x' when executing action a in state x
function POMDPs.transition(DP::DrivePOMDP, Ss::Sts, Aego::Symbol)
    acc_k = 0.0 # we can also use Aref
    akV = [e for e in DP.Aset.min:DP.Aset.max]
    acc_distribution = DiscreteND1D(acc_k, DP.Aset.comfort/2, akV, 1.0)
    PD = zeros(POMDPs.n_states(DP))
    for (i, ak) in enumerate(akV)
        if acc_distribution[i] > 1.0e-2
            result = IDMtransit(DP, Ss, Aego, ak)
            SIndex = 0
            try
                SIndex = POMDPs.stateindex(DP, Sts(result.Ego, result.Other)) # round is made in stateindex()
            catch
                error("Transition Result error: $result")
            end
            PD[SIndex] += acc_distribution[i]
        end
    end
    normalize!(PD, 1)
    #return DiscreteBelief(DP, PD) # toooooo slow
    return SparseCat(DP.SSpace, PD)
end

function IDMtransit(DP::DrivePOMDP, Ss::Sts, Aego::Symbol, acc_k::Float64)

    if Ss.Ego.s <= 0.0 && Ss.Ego.v <= 0.0 && Aego == :giveup # avoid transition of edge states going out of state range
        Other = sv_next(Ss.Other, acc_k, DP.Δt)
        sv_boundry!(DP, Other)
        Ego = deepcopy(Ss.Ego)
        return (Ego=Ego, Other=Other, accs=zeros(Int64(DP.Δt/0.1)))
    end

    Ego = deepcopy(Ss.Ego)
    Other = deepcopy(Ss.Other)
    accVec = Vector{Float64}()
    overlap = false
    dist1 = 0.0
    distk = 0.0

    if haskey(DP.Routes[Ego.r].intersect_Infos, Other.r)
        PtDistVec1 = DP.Routes[Ego.r].intersect_Infos[Other.r]
        if length(PtDistVec1) > 1 && PtDistVec1[2][1] == "overlap"
            dist1 = PtDistVec1[2][3]
            distk = DP.Routes[Other.r].intersect_Infos[Ego.r][2][3]
            overlap = true
        end
    end

    step = Int64(DP.Δt/0.1)
    for i in 1:step+1 # update the state of EgoCar and OtherCar
        acc_ego = 0.0
        if overlap
            if Ego.s <= DP.Stopline
                acc_ego = AccCalculate(DP, Ego, Aego) # only based on stopline, the state of other car is not directly used.
            else # Ego.s >= DP.Stopline
                if Aego == :giveup
                    Δs = Other.s - Ego.s + dist1 - distk
                    if Δs < DP.Smin # ego vehicle can not follow other vehicle because too close/ego car is in front of other car.
                        if Ego.v > 0 # vehicle should not go back
                            acc_ego = DP.Aset.min
                        else
                            acc_ego = 0.0
                        end
                    else # Δs >= DP.Smin; ego should follow other
                        acc_ego = AccCalculate(DP, Ego, Other, Aego, Δs)
                    end
                else # Aego == :takeover -> drives freely
                    Vref = DP.Routes[Ego.r].Vref[min(UInt16(floor(Ego.s/DP.Δs)+1), 21)]
                    acc_ego = IDM(Vego=Ego.v, Vfront=Vref, Vref=Vref, Snet=Inf, T=1.0, Amax=DP.Aset.max, Bdec=DP.Aset.comfort, Smin=DP.Smin)
                end
            end
        else # no overlap
            acc_ego = AccCalculate(DP, Ego, Aego)
        end
        acc_ego = Acclimit!(DP, acc_ego)
        push!(accVec, acc_ego)
        #@show acc_ego
        if i <= step
            sv_next!(Ego, acc_ego, 0.1)
            Ego.s = round(Ego.s, 1e-13)
            sv_boundry!(DP, Ego)
            sv_next!(Other, acc_k, 0.1)
            Other.s = round(Other.s, 1e-13)
            sv_boundry!(DP, Other)
        end
        #@show Ego
        #@show Other
    end

    #CSRound!(Ego, DP.Δs, DP.Δv)
    #CSRound!(Other, DP.Δs, DP.Δv)
    return (Ego=Ego, Other=Other, accs=accVec)
end


text:
"""
A simple implementation of calculating edit distance of two strings
"""
function edit_distance(s1::ASCIIString, s2::ASCIIString)
    m=length(s1)+1
    n=length(s2)+1
    i=0
    j=0
    tbl = Array{Int64}(m,n)
    for i=1:m
        tbl[i,1]=i-1
    end
    for j=1:n
        tbl[1,j]=j-1
    end
    for i=2:m
        for j=2:n
            cost = s1[i-1] == s2[j-1]?0:1
            tbl[i,j] = min(tbl[i, j-1]+1, tbl[i-1, j]+1, tbl[i-1, j-1]+cost)
        end
    end
    return tbl[i,j]
end

function hamming_distance(s1::ASCIIString, s2::ASCIIString)
    d = 0
    len = min(length(s1), length(s2))
    for i in 1:len
        if s1[i] != s2[i]
            d += 1
        end
    end
    return d
end

function edit_distance(s1::Sequence, s2::Sequence)
    return edit_distance(s1.seq, s2.seq)
end

function hamming_distance(s1::Sequence, s2::Sequence)
    return hamming_distance(s1.seq, s2.seq)
end


text:
using DashBootstrapComponents, DashHtmlComponents

cards = dbc_cardgroup([
    dbc_card(
        dbc_cardbody([
            html_h5("Card 1", className="card-title"),
            html_p(
                "This card has some text content, which is a little " *
                "bit longer than the second card.",
                className="card-text",
            ),
            dbc_button("Click here", color="success", className="mt-auto"),
        ]),
    ),
    dbc_card(
        dbc_cardbody([
            html_h5("Card 2", className="card-title"),
            html_p("This card has some text content.", className="card-text"),
            dbc_button("Click here", color="warning", className="mt-auto"),
        ]),
    ),
    dbc_card(
        dbc_cardbody([
            html_h5("Card 3", className="card-title"),
            html_p(
                "This card has some text content, which is longer " *
                "than both of the other two cards, in order to " *
                "demonstrate the equal height property of cards in a " *
                "card group.",
                className="card-text",
            ),
            dbc_button("Click here", color="danger", className="mt-auto"),
        ]),
    ),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
pyfraction(x::Rational) = pyfraction(numerator(x), denominator(x))
pyfraction(x, y) = pyfractiontype(x, y)
pyfraction(x) = pyfractiontype(x)
pyfraction() = pyfractiontype()
export pyfraction

# works for any collections.abc.Rational
function pyconvert_rule_fraction(::Type{R}, x::Py, ::Type{Rational{T0}}=Utils._type_lb(R), ::Type{Rational{T1}}=Utils._type_ub(R)) where {R<:Rational,T0,T1}
    a = @pyconvert_and_del(Utils._typeintersect(Integer, T1), x.numerator)
    b = @pyconvert_and_del(Utils._typeintersect(Integer, T1), x.denominator)
    a, b = promote(a, b)
    T2 = Utils._promote_type_bounded(T0, typeof(a), typeof(b), T1)
    pyconvert_return(Rational{T2}(a, b))
end

# works for any collections.abc.Rational
function pyconvert_rule_fraction(::Type{T}, x::Py) where {T<:Number}
    pyconvert_tryconvert(T, @pyconvert(Rational{<:Integer}, x))
end


text:
include("./imports.jl")

@testset ExtendedTestSet "Optimizer and OptimParams" begin
    opt = ADAM(0.01)
    @test getoptimparam(opt, LR) == 0.01
    setoptimparam!(opt, LR, 0.001)
    @test getoptimparam(opt, LR) == 0.001
end


text:
function validate_uci_move(uci_move::String)
    if  ~occursin(r"[a-h][1-8] [a-h][1-8] [queen, rook, night, bishop, none]",
        uci_move)
        println("wrong move format, try again")
        return false
    end
    return true
end

function ask_for_move()
    while true
        println("Enter move:")
        user_move = readline(stdin)

        if user_move == "quit"
            return "xx", "xx"
        end

        if validate_uci_move(user_move)
            return split(user_move, ' ')
        end
    end
    return "xx xx"
end

function validate_promotion_type(promotion_type::String)
    if ~occursin(r"[queen, rook, night, bishop]", promotion_type)
        println("invalid promotion type, try again")
        return false
    end
    return true
end

function ask_for_promotion()
    while true
        println("Promove to [queen], [rook], k[night], or [bishop]?")
        promotion_type = readline(stdin)
        if validate_promotion_type(promotion_type)
            return promotion_type
        end
    end
end

function play(human_color::String="white")
    b = set_board()
    check = false

    while true
        pretty_print(b)

        moves = get_all_valid_moves(b, b.player_color)
        if length(moves) == 0
            if check
                println("check mate!")
                return
            else
                println("stalemate...")
                return
            end
        end
        println(b.player_color," to move")

        if b.player_color == human_color
            while true
                s, t, p = ask_for_move()
                if s == "xx"
                    return
                end

                for move in moves
                    if move.source == PGN2UINT[s] && move.target == PGN2UINT[t]
                        if move.promotion_type != "none"
                            promotion_type = ask_for_promotion()
                        end

                        b = move_piece(b, move, b.player_color)
                        b = update_attacked(b)
                        b = update_castling_rights(b)
                        @goto next_move
                    end
                end
                println("Move not available, try again")
            end
        else
            b = move_piece(b, moves[rand(1:length(moves))], b.player_color)
            b = update_attacked(b)
            b = update_castling_rights(b)
        end

        @label next_move
        b.player_color = change_color(b.player_color)
        if king_in_check(b, b.player_color)
            check = true
            println("check!")
        end
    end
end