text:
type Gaussian1DNonstandardSwaptionEngine{G <: Gaussian1DModel, Y <: YieldTermStructure, P <: GaussianProbabilities} <: PricingEngine
  model::G
  integrationPoints::Int
  stddevs::Float64
  extrapolatePayoff::Bool
  flatPayoffExtrapolation::Bool
  oas::Quote
  discountCurve::Y
  probabilities::P
end

Gaussian1DNonstandardSwaptionEngine{G <: Gaussian1DModel, Y <: YieldTermStructure, P <: GaussianProbabilities}(model::G, integrationPoints::Int = 64, stddevs::Float64 = 7.0,
                        extrapolatePayoff::Bool = true, flatPayoffExtrapolation::Bool = false, oas::Quote = Quote(-1.0), discountCurve::Y = NullYieldTermStructure(),
                        probabilities::P = NoneProbabilities()) =
                        Gaussian1DNonstandardSwaptionEngine{G, Y, P}(model, integrationPoints, stddevs, extrapolatePayoff, flatPayoffExtrapolation, oas, discountCurve, probabilities)

# methods #
function _calculate!(pe::Gaussian1DNonstandardSwaptionEngine, swaption::NonstandardSwaption)
  isa(swaption.delivery, SettlementPhysical) || error("cash-settled swaptions not yet implemented")
  settlement = reference_date(pe.model.ts)

  if swaption.exercise.dates[end] <= settlement
    swaption.results.value = 0.0 # swaption is expired
    return swaption
  end

  # rebatedExercise = RebatedExercise(swaption.exercise)
  idx = length(swaption.exercise.dates)

  minIdxAlive = upper_bound(swaption.exercise.dates, settlement)

  swap = swaption.swap
  optType = isa(swap.swapT, Payer) ? Call() : Put()
  fixedSchedule = swap.fixedSchedule
  floatingSchedule = swap.floatSchedule

  npv0 = zeros(2 * pe.integrationPoints + 1)
  npv1 = zeros(2 * pe.integrationPoints + 1)
  z = y_grid(pe.model, pe.stddevs, pe.integrationPoints)
  p = zeros(length(z))

  # for probability computation
  npvp0 = Vector{Vector{Float64}}()
  npvp1 = Vector{Vector{Float64}}()
  if !isa(pe.probabilities, NoneProbabilities)
    for i = 1:(idx - (minIdxAlive - 1) + 1)
      npvTmp0 = zeros(2 * pe.integrationPoints + 1)
      npvTmp1 = zeros(2 * integrationPoints + 1)
      push!(npvp0, npvTmp0)
      push!(npvp1, npvTmp1)
    end
  end
  # end probability computation

  expiry0 = expiry1 = Date()
  expiry0Time = expiry1Time = -1.0

  while true
    if idx == minIdxAlive - 1
      expiry0 = settlement
    else
      expiry0 = swaption.exercise.dates[idx]
    end

    expiry0Time = max(time_from_reference(pe.model.ts, expiry0), 0.0)
    j1 = upper_bound(fixedSchedule.dates, expiry0 - Dates.Day(1))
    k1 = upper_bound(floatingSchedule.dates, expiry0 - Dates.Day(1))

    for k = 1:(expiry0 > settlement ? length(npv0) : 1)
      price = 0.0

      if expiry1Time != -1.0
        zSpreadDf = pe.oas.value == -1.0 ? 1.0 : exp(-pe.oas.value * (expiry1Time - expiry0Time))

        yg = y_grid(pe.model, pe.stddevs, pe.integrationPoints, expiry1Time, expiry0Time, expiry0 > settlement ? z[k] : 0.0)
        payoff0 = CubicInterpolation(Spline(), Lagrange(), Lagrange(), z, npv1)
        for i in eachindex(yg)
          p[i] = payoff0(yg[i])
        end

        payoff1 = CubicInterpolation(Spline(), Lagrange(), Lagrange(), z, p)
        for i = 1:length(z) - 1
          price += gaussian_shifted_polynomial_integral(pe.model, 0.0, payoff1.c[i], payoff1.b[i], payoff1.a[i], p[i], z[i], z[i], z[i + 1]) * zSpreadDf
        end
        if pe.extrapolatePayoff
          if pe.flatPayoffExtrapolation
            price += gaussian_shifted_polynomial_integral(pe.model, 0.0, 0.0, 0.0, 0.0, p[length(z) - 1], z[end - 1], z[end], 100.0) * zSpreadDf
            price += gaussian_shifted_polynomial_integral(pe.model, 0.0, 0.0, 0.0, 0.0, p[1], z[1], -100.0, z[1]) * zSpreadDf
          else
            if isa(optType, Call)
              price += gaussian_shifted_polynomial_integral(pe.model, 0.0, payoff1.c[length(z) - 1], payoff1.b[length(z) - 1], payoff1.a[length(z) - 1],
                        p[length(z) - 1], z[end - 1], z[end], 100.0) * zSpreadDf
            else
              price += gaussian_shifted_polynomial_integral(pe.model, 0.0, payoff1.c[1], payoff1.b[1], payoff1.a[1], p[1], z[1], -100.0, z[1]) * zSpreadDf
            end
          end
        end
      end
      npv0[k] = price
      # for probability computation
      if !isa(pe.probabilities, NoneProbabilities)
        for m = 1:length(npvp0)
          price = 0.0
          if expiry1Time != -1.0
            zSpreadDf = pe.oas.value == -1.0 ? 1.0 : exp(-pe.oas.value * (expiry1Time - expiry0Time))
            yg = y_grid(pe.model, pe.stddevs, pe.integrationPoints, expiry1Time, expiry0Time, expiry0 > settlement ? z[k] : 0.0)
            payoff0 = CubicInterpolation(Spline(), Lagrange(), Lagrange(), z, npvp1[m])
            for i in eachindex(yg)
              p[i] = payoff0(yg[i])
            end

            payoff1 = CubicInterpolation(Spline(), Lagrange(), Lagrange(), z, p)
            for i = 1:length(z) - 1
              price += gaussian_shifted_polynomial_integral(pe.model, 0.0, payoff1.c[i], payoff1.b[i], payoff1.a[i], p[i], z[i], z[i], z[i + 1]) * zSpreadDf
            end

            if pe.extrapolatePayoff
              if pe.flatPayoffExtrapolation
                price += gaussian_shifted_polynomial_integral(pe.model, 0.0, 0.0, 0.0, 0.0, p[length(z) - 1], z[end - 1], z[end], 100.0) * zSpreadDf
                price += gaussian_shifted_polynomial_integral(pe.model, 0.0, 0.0, 0.0, 0.0, p[1], z[1], -100.0, z[1]) * zSpreadDf
              else
                if isa(optType, Call)
                  price += gaussian_shifted_polynomial_integral(pe.model, 0.0, payoff1.c[length(z) - 1], payoff1.b[length(z) - 1], payoff1.a[length(z) - 1],
                            p[length(z) - 1], z[end - 1], z[end], 100.0) * zSpreadDf
                else
                  price += gaussian_shifted_polynomial_integral(pe.model, 0.0, payoff1.c[1], payoff1.b[1], payoff1.a[1], p[1], z[1], -100.0, z[1]) * zSpreadDf
                end
              end
            end
          end
          npvp0[m][k] = price
        end
      end
      # end probability computation
      if expiry0 > settlement
        floatingLegNPV = 0.0
        for l = k1:length(get_floating_coupons(swaption.swap))
          zSpreadDf = pe.oas.value == -1.0 ? 1.0 : exp(-pe.oas.value * (year_fraction(pe.model.ts.dc, expiry0, get_floating_pay_dates(swaption.swap)[l])))
          if swaption.swap.args.floatingIsRedemptionFlow[l]
            amount = get_floating_coupons(swaption.swap)[l]
          else
            amount = swaption.swap.floatingNominal[l] * get_floating_accrual_times(swaption.swap)[l] * (swaption.swap.args.floatingGearings[l] *
                      forward_rate(pe.model, swaption.swap.legs[2].coupons[l].fixingDate, expiry0, z[k], swaption.swap.iborIndex) +
                      get_floating_spreads(swaption.swap)[l])
          end
          floatingLegNPV += amount * zerobond(pe.model, get_floating_pay_dates(swaption.swap)[l], expiry0, z[k], pe.discountCurve) * zSpreadDf
        end
        fixedLegNPV = 0.0
        for l = j1:length(get_fixed_coupons(swaption.swap))
          zSpreadDf = pe.oas.value == -1.0 ? 1.0 : exp(-pe.oas.value * (year_fraction(pe.model.ts.dc, expiry0, get_fixed_pay_dates(swaption.swap)[l])))
          fixedLegNPV += get_fixed_coupons(swaption.swap)[l] * zerobond(pe.model, get_fixed_pay_dates(swaption.swap)[l], expiry0, z[k], pe.discountCurve) *
                        zSpreadDf
        end
        rebate = 0.0
        zSpreadDf = 1.0
        rebateDate = expiry0
        # TODO rebate Exercise code
        exerciseValue = ((isa(optType, Call) ? 1.0 : -1.0) * (floatingLegNPV - fixedLegNPV) + rebate *
                          zerobond(pe.model, rebateDate, expiry0, z[k], pe.discountCurve) * zSpreadDf) /
                        numeraire(pe.model, expiry0Time, z[k], pe.discountCurve)

        # probability computation
        if !isa(pe.probabilities, NoneProbabilities)
          if idx == length(swaption.exercise.dates) # if true we are at the latest date so we init the no call prob
            npvp0[end][k] = isa(pe.probabilities, NaiveProbabilities) ? 1.0 : 1.0 / (zerobond(pe.model, expiry0Time, 0.0, 0.0, pe.discountCurve) * numeraire(pe.model, expiry0, z[k], pe.discountCurve))
          end

          if exerciseValue >= npv0[k]
            npvp0[idx - (minIdxAlive - 1)][k] = isa(pe.probabilities, NaiveProbabilities) ? 1.0 : (zerobond(pe.model, expiry0Time, 0.0, 0.0, pe.discountCurve) * numeraire(pe.model, expiry0Time, z[k], pe.discountCurve))
            for ii = idx - (minIdxAlive - 1) + 1:length(npvp0)
              npvp0[ii][k] = 0.0
            end
          end
        end
        # end probability computation
        npv0[k] = max(npv0[k], exerciseValue)
      end
    end

    npv1, npv0 = npv0, npv1

    # for probability computation
    if !isa(pe.probabilities, NoneProbabilities)
      for i in eachindex(npvp0)
        npvp1[i], npvp0[i] = npvp0[i], npvp1[i]
      end
    end
    # end probability computation

    expiry1 = expiry0
    expiry1Time = expiry0Time
    idx -= 1
    if idx < minIdxAlive - 1
      break
    end
  end

  swaption.results.value = npv1[1] * numeraire(pe.model, 0.0, 0.0, pe.discountCurve)

  # for probability computation
  if !isa(pe.probabilities, NoneProbabilities)
    prob = zeros(length(npvp0))
    for i in eachindex(npvp0)
      prob[i] = npvp[i][1] * (isa(pe.probabilities, NaiveProbabilities) ? 1.0 : numeraire(pe.model, 0.0, 0.0, pe.discountCurve))
    end
    swaption.results.additionalResults["probabilities"] = prob
  end

  return swaption
end

underlying_last_date(::Gaussian1DNonstandardSwaptionEngine, swaption::NonstandardSwaption) = get_fixed_pay_dates(swaption.swap)[end]

function calibration_basket(swaptionEngine::Gaussian1DNonstandardSwaptionEngine, swaption::NonstandardSwaption, swapIndex::SwapIndex,
                            swaptionVol::SwaptionVolatilityStructure, basketType::NaiveBasketType)

  n = length(swaption.exercise.dates)
  minIdxAlive = upper_bound(swaption.exercise.dates, settings.evaluation_date)
  minIdxAlive = minIdxAlive == 0 ? 1 : minIdxAlive
  result = Vector{SwaptionHelper}(n - (minIdxAlive - 1))

  # rebEx = RebatedExercise(exercise)

  for i = minIdxAlive:n
    expiry = swaption.exercise.dates[i]
    rebate = 0.0
    rebateDate = expiry

    swapLength = year_fraction(swaptionVol.dc, value_date(swapIndex, expiry), underlying_last_date(swaptionEngine, swaption))
    sec = smile_section(swaptionVol, expiry, Dates.Month(round(Int, floor(swapLength * 12.0 + 0.5))))
    atmStrike = sec.atmLevel
    if atmStrike == -1.0
      atmVol = volatility(sec, 0.03)
    else
      atmVol = volatility(sec, atmStrike)
    end

    shift = sec.shift
    ts = swapIndex.exogenousDiscount ? swapIndex.discount : swapIndex.iborIndex.ts
    helper = SwaptionHelper(expiry, underlying_last_date(swaptionEngine, swaption), Quote(atmVol), swapIndex.iborIndex, TenorPeriod(swapIndex.fixedLegTenor),
                            swapIndex.fixedLegDayCount, swapIndex.iborIndex.dc,  ts, NullSwaptionEngine(), -1.0, 1.0, shift)

    result[i] = helper
  end

  return result
end


text:
""" Generic data structures for C-sets (presheaves).
"""
module CSets
export AbstractCSet, AbstractCSetType, CSet, CSetType,
  nparts, has_part, subpart, has_subpart, incident,
  add_part!, add_parts!, copy_parts!, set_subpart!, set_subparts!

using Compat
using LabelledArrays, StaticArrays

using ...Present
using ...Theories: Category, dom, codom

# C-set data types
##################

const NamedSVector{Names,T,N} = SLArray{Tuple{N},T,1,N,Names}

""" Abstract type for C-sets (presheaves).

The type parameters are:
- `Ob`: tuple of symbols naming objects
- `Hom`: tuple of symbols naming morphims
- `Dom`: tuple of integers giving domain of each morphism
- `Codom`: tuple of integers giving codomain of each morphism
- `Data`: tuple of symbols naming extra morphisms for data/attributes
- `DataDom`: tuple of integers giving domain of each data morphism

For example, a graph with one vertex attribute could be represented as
`AbstractCSet{(:V,:E),(:src,:tgt),(2,2),(1,1),(:vattr,),(1,)}`.

See also: [`CSet`](@ref) and [`CSetType`](@ref).
"""
abstract type AbstractCSet{Ob,Hom,Dom,Codom,Data,DataDom} end

""" Data type for C-sets (presheaves).

Instead of filling out the type parameters manually, we recommend using the
function [`CSetType`](@ref) to generate a `CSet` type from a presentation of a
category. Nevertheless, the first six type parameters are documented at
[`AbstractCSet`](@ref). The remaining type parameters are an implementation
detail and should be ignored.

Following LightGraphs.jl, the incidence vectors, stored in the `incidence`
field, are kept in sorted order. To ensure consistency, no field of the struct
should ever be mutated directly.
"""
mutable struct CSet{Ob,Hom,Dom,Codom,Data,DataDom,Index,DataIndex,IsUnique} <:
       AbstractCSet{Ob,Hom,Dom,Codom,Data,DataDom}
  nparts::NamedSVector{Ob,Int}
  subparts::NamedSVector{Hom,Vector{Int}}
  indices::NamedSVector{Index,Vector{Vector{Int}}}
  data::NamedTuple{Data}
  data_indices::NamedTuple{DataIndex}
end

function CSet{Ob,Hom,Dom,Codom,Data,DataDom,Index,DataIndex,IsUnique}(;
    kw...) where {Ob,Hom,Dom,Codom,Data,DataDom,Index,DataIndex,IsUnique}
  CSet{Ob,Hom,Dom,Codom,Data,DataDom,Index,DataIndex,IsUnique}((; kw...))
end
function CSet{Ob,Hom,Dom,Codom,Data,DataDom,Index,DataIndex,IsUnique}(
    datatypes::NamedTuple{Data}) where
    {Ob,Hom,Dom,Codom,Data,DataDom,Index,DataIndex,IsUnique}
  NOb, NHom, NIndex = length(Ob), length(Hom), length(Index)
  CSet{Ob,Hom,Dom,Codom,Data,DataDom,Index,DataIndex,IsUnique}(
    NamedSVector{Ob,Int,NOb}(zeros(SVector{NOb,Int})),
    NamedSVector{Hom,Vector{Int},NHom}(Tuple(Int[] for i in 1:NHom)),
    NamedSVector{Index,Vector{Vector{Int}},NIndex}(
      Tuple(Vector{Int}[] for i in 1:NIndex)),
    NamedTuple{Data}(Union{T,Missing}[] for T in datatypes),
    NamedTuple{DataIndex}(
      Dict{datatypes[name], is_unique ? Int : Vector{Int}}()
      for (name, is_unique) in zip(DataIndex, IsUnique)))
end

""" Generate an abstract C-set type from a presentation of a category.

See also: [`CSetType`](@ref).
"""
function AbstractCSetType(pres::Presentation{Category}; data=())
  # Only include `Data` and `DataDom` params if there are data morphisms.
  Ts = CSetTypeParams(pres, data=data)
  AbstractCSet{(isempty(Ts[5]) ? Ts[1:4] : Ts[1:6])...}
end

""" Generate a C-set data type from a presentation of a category.

See also: [`AbstractCSetType`](@ref).
"""
function CSetType(pres::Presentation{Category}; kw...)
  CSet{CSetTypeParams(pres; kw...)...}
end

function CSetTypeParams(pres::Presentation{Category};
                        data=(), index=(), unique_index=())
  obs, homs = generators(pres, :Ob), generators(pres, :Hom)
  data_obs, obs = separate(ob -> nameof(ob) ∈ data, obs)
  data_homs, homs = separate(hom -> codom(hom) ∈ data_obs, homs)
  data_index, index = separate(name -> name ∈ nameof.(data_homs),
                               sort!(index ∪ unique_index))
  ob_num = ob -> findfirst(obs .== ob)::Int
  (Tuple(nameof.(obs)), Tuple(nameof.(homs)),
   Tuple(@. ob_num(dom(homs))), Tuple(@. ob_num(codom(homs))),
   Tuple(nameof.(data_homs)), Tuple(@. ob_num(dom(data_homs))),
   Tuple(index), Tuple(data_index), Tuple(k ∈ unique_index for k in data_index))
end
separate(f, a::AbstractArray) = (i = f.(a); (a[i], a[.!i]))

function Base.:(==)(x1::T, x2::T) where T <: CSet
  # The subpart and data indices are redundant, so need not be compared.
  x1.nparts == x2.nparts && x1.subparts == x2.subparts && x1.data == x2.data
end

function Base.copy(cset::T) where T <: CSet
  T(cset.nparts, map(copy, cset.subparts), map(copy, cset.indices),
    map(copy, cset.data), map(copy, cset.data_indices))
end

Base.empty(cset::T) where T <: CSet = T(map(eltype, cset.data))

# C-set interface
#################

""" Number of parts of given type in a C-set.
"""
nparts(cset::CSet, type::Symbol) = cset.nparts[type]

""" Whether a C-set has a part with the given name.
"""
has_part(cset::CSet, type::Symbol) = _has_part(cset, Val(type))

@generated _has_part(cset::CSet{obs}, ::Val{type}) where {obs,type} =
  type ∈ obs

has_part(cset::CSet, type::Symbol, part::Int) = 1 <= part <= nparts(cset, type)
has_part(cset::CSet, type::Symbol, part::AbstractVector{Int}) =
  let n=nparts(cset, type); [ 1 <= x <= n for x in part ] end

""" Get subpart of part in C-set.

Both single and vectorized access are supported.
"""
subpart(cset::CSet, part, name::Symbol) = subpart(cset, name)[part]
subpart(cset::CSet, name::Symbol) = _subpart(cset, Val(name))

@generated function _subpart(cset::T, ::Val{name}) where
    {name, obs,homs,doms,codoms,data, T <: CSet{obs,homs,doms,codoms,data}}
  if name ∈ homs
    :(cset.subparts.$name)
  elseif name ∈ data
    :(cset.data.$name)
  else
    throw(KeyError(name))
  end
end

""" Whether a C-set has a subpart with the given name.
"""
has_subpart(cset::CSet, name::Symbol) = _has_subpart(cset, Val(name))

@generated function _has_subpart(cset::T, ::Val{name}) where
    {name, obs,homs,doms,codoms,data, T <: CSet{obs,homs,doms,codoms,data}}
  name ∈ homs || name ∈ data
end

""" Get superparts incident to part in C-set.
"""
incident(cset::CSet, part, name::Symbol) = _incident(cset, part, Val(name))

@generated function _incident(cset::T, part, ::Val{name}) where
    {name, obs,homs,doms,codoms,data,data_doms,indexed,data_indexed,
     T <: CSet{obs,homs,doms,codoms,data,data_doms,indexed,data_indexed}}
  if name ∈ indexed
    :(cset.indices.$name[part])
  elseif name ∈ data_indexed
    :(get_data_index.(Ref(cset.data_indices.$name), part))
  else
    throw(KeyError(name))
  end
end

""" Add part of given type to C-set, optionally setting its subparts.

Returns the ID of the added part.

See also: [`add_parts!`](@ref).
"""
function add_part!(cset::CSet, type::Symbol, args...; kw...)
  part = only(_add_parts!(cset, Val(type), 1))
  set_subparts!(cset, part, args...; kw...)
  part
end

""" Add parts of given type to C-set, optionally setting their subparts.

Returns the range of IDs for the added parts.

See also: [`add_part!`](@ref).
"""
function add_parts!(cset::CSet, type::Symbol, n::Int, args...; kw...)
  parts = _add_parts!(cset, Val(type), n)
  set_subparts!(cset, parts, args...; kw...)
  parts
end

@generated function _add_parts!(cset::T, ::Val{type}, n::Int) where
    {type, obs,homs,doms,codoms,data,data_doms,indexed,
     T <: CSet{obs,homs,doms,codoms,data,data_doms,indexed}}
  ob_num = findfirst(obs .== type)::Int
  in_homs = [ homs[i] for (i, codom) in enumerate(codoms) if codom == ob_num ]
  out_homs = [ homs[i] for (i, dom) in enumerate(doms) if dom == ob_num ]
  data_homs = [ data[i] for (i, dom) in enumerate(data_doms) if dom == ob_num ]
  indexed_homs = filter(hom -> hom ∈ indexed, in_homs)
  # Question: The three loops could (should?) be unrolled. Or is the Julia
  # compiler smart enough to do this on its own?
  quote
    if n == 0; return 1:0 end
    nparts = cset.nparts.$type + n
    cset.nparts = SLVector(cset.nparts; $type=nparts)
    start = nparts - n + 1
    for name in $(Tuple(out_homs))
      resize!(cset.subparts[name], nparts)
      @inbounds cset.subparts[name][start:nparts] .= 0
    end
    for name in $(Tuple(indexed_homs))
      index = cset.indices[name]
      resize!(index, nparts)
      @inbounds for part in start:nparts; index[part] = Int[] end
    end
    for name in $(Tuple(data_homs))
      resize!(cset.data[name], nparts)
      @inbounds cset.data[name][start:nparts] .= missing
    end
    start:nparts
  end
end

""" Copy parts from one C-set to another.

All subparts between the selected parts, including any attached data, are
preserved. Thus, if the selected parts form a sub-C-set, then the whole
sub-C-set is preserved. On the other hand, if the selected parts do *not* form a
sub-C-set, then some copied parts will have undefined subparts.
"""
copy_parts!(cset::CSet, from::CSet; kw...) = copy_parts!(cset, from, (; kw...))

copy_parts!(cset::CSet, from::CSet, type::Symbol, parts) =
  copy_parts!(cset, from, (; type => parts))

@generated function copy_parts!(cset::T, from::T, parts::NamedTuple{types}) where
    {types, obs,homs,doms,codoms, T <: CSet{obs,homs,doms,codoms}}
  obnums = [ findfirst(obs .== type)::Int for type in types ]
  in_obs, out_homs = Symbol[], Tuple{Symbol,Symbol,Symbol}[]
  for (hom, dom, codom) in zip(homs, doms, codoms)
    if dom ∈ obnums && codom ∈ obnums
      push!(in_obs, obs[codom])
      push!(out_homs, (hom, obs[dom], obs[codom]))
    end
  end
  in_obs = Tuple(unique!(in_obs))
  quote
    newparts = NamedTuple{$types}(tuple($(map(types) do type
      :(_copy_parts_data!(cset, from, Val($(QuoteNode(type))), parts.$type))
    end...)))
    partmaps = NamedTuple{$in_obs}(tuple($(map(in_obs) do type
      :(Dict{Int,Int}(zip(parts.$type, newparts.$type)))
    end...)))
    for (name, dom, codom) in $(Tuple(out_homs))
      for (p, newp) in zip(parts[dom], newparts[dom])
        q = subpart(from, p, name)
        newq = get(partmaps[codom], q, nothing)
        if !isnothing(newq)
          set_subpart!(cset, newp, name, newq)
        end
      end
    end
    newparts
  end
end

""" Copy parts of a single type from one C-set to another.

Any data attached to the parts is also copied.
"""
@generated function _copy_parts_data!(cset::T, from::T, ::Val{type}, parts) where
    {type, obs,homs,doms,codoms,data,data_doms,
     T <: CSet{obs,homs,doms,codoms,data,data_doms}}
  obnum = findfirst(obs .== type)::Int
  data_homs = [ data[i] for (i, dom) in enumerate(data_doms) if dom == obnum ]
  quote
    newparts = add_parts!(cset, $(QuoteNode(type)), length(parts))
    for name in $(Tuple(data_homs))
      set_subpart!(cset, newparts, name, from.data[name][parts])
    end
    newparts
  end
end

""" Mutate subpart of a part in a C-set.

Both single and vectorized assignment are supported.

See also: [`set_subparts!`](@ref).
"""
set_subpart!(cset::CSet, part::Int, name::Symbol, subpart) =
  _set_subpart!(cset, part, Val(name), subpart)

function set_subpart!(cset::CSet, part::AbstractVector{Int},
                      name::Symbol, subpart)
  broadcast(part, subpart) do part, subpart
    _set_subpart!(cset, part, Val(name), subpart)
  end
end

set_subpart!(cset::CSet, ::Colon, name::Symbol, subpart) =
  set_subpart!(cset, name, subpart)
set_subpart!(cset::CSet, name::Symbol, new_subpart) =
  set_subpart!(cset, 1:length(subpart(cset, name)), name, new_subpart)

@generated function _set_subpart!(cset::T, part::Int, ::Val{name}, subpart) where
    {name, obs,homs,doms,codoms,data,data_doms,indexed,data_indexed,
     T <: CSet{obs,homs,doms,codoms,data,data_doms,indexed,data_indexed}}
  if name ∈ homs
    codom = obs[codoms[findfirst(name .== homs)]]
  end
  if name ∈ indexed
    quote
      @assert 0 <= subpart <= cset.nparts.$codom
      old = cset.subparts.$name[part]
      cset.subparts.$name[part] = subpart
      if old > 0
        deletesorted!(cset.indices.$name[old], part)
      end
      if subpart > 0
        insertsorted!(cset.indices.$name[subpart], part)
      end
    end
  elseif name ∈ homs
    quote
      @assert 0 <= subpart <= cset.nparts.$codom
      cset.subparts.$name[part] = subpart
    end
  elseif name ∈ data_indexed
    quote
      old = cset.data.$name[part]
      cset.data.$name[part] = subpart
      if !ismissing(old)
        unset_data_index!(cset.data_indices.$name, old, part)
      end
      if !ismissing(subpart)
        set_data_index!(cset.data_indices.$name, subpart, part)
      end
    end
  elseif name ∈ data
    :(cset.data.$name[part] = subpart)
  else
    throw(KeyError(name))
  end
end

""" Mutate subparts of a part in a C-set.

Both single and vectorized assignment are supported.

See also: [`set_subpart!`](@ref).
"""
set_subparts!(cset::CSet, part; kw...) = set_subparts!(cset, part, (; kw...))

function set_subparts!(cset::CSet, part, subparts)
  for (name, subpart) in pairs(subparts)
    set_subpart!(cset, part, name, subpart)
  end
end

""" Look up key in C-set data index.
"""
get_data_index(d::AbstractDict{K,Int}, k::K) where K =
  get(d, k, nothing)
get_data_index(d::AbstractDict{K,<:AbstractVector{Int}}, k::K) where K =
  get(d, k, 1:0)

""" Set key and value for C-set data index.
"""
function set_data_index!(d::AbstractDict{K,Int}, k::K, v::Int) where K
  if haskey(d, k)
    error("Key $k already defined in unique index")
  end
  d[k] = v
end
function set_data_index!(d::AbstractDict{K,<:AbstractVector{Int}},
                         k::K, v::Int) where K
  insertsorted!(get!(d, k) do; Int[] end, v)
end

""" Unset key and value from C-set data index.
"""
function unset_data_index!(d::AbstractDict{K,Int}, k::K, v::Int) where K
  @assert pop!(d, k) == v
end
function unset_data_index!(d::AbstractDict{K,<:AbstractVector{Int}},
                           k::K, v::Int) where K
  vs = d[k]
  deletesorted!(vs, v)
  if isempty(vs)
    delete!(d, k)
  end
end

""" Insert into sorted vector, preserving the sorting.
"""
function insertsorted!(a::AbstractVector, x)
  insert!(a, searchsortedfirst(a, x), x)
end

""" Delete one occurrence of value from sorted vector.
"""
function deletesorted!(a::AbstractVector, x)
  i = searchsortedfirst(a, x)
  @assert i <= length(a) && a[i] == x "Value $x is not contained in $a"
  deleteat!(a, i)
end

end

ext: jl


text:
"""
    delete!(tree::RTree, pt::Point, [id])
    delete!(tree::RTree, br::Rect, [id])

Deletes the value identified by `br` bounding box (or point `pt`) and
the `id` (if tree elements support `HasID` trait) from the `tree`.
"""
function Base.delete!(tree::RTree{T,N}, br::Rect{T,N}, id::Any = nothing) where {T,N}
    leafx = findfirst(tree, br, id)
    leafx === nothing && __spatial_keyerror(eltype(tree), br, id)
    _detach!(leafx[1], leafx[2], tree)
    tmpdetached = Vector{nodetype(tree)}() # FIXME use pool, Union{Leaf,Branch} ?
    _condense!(leafx[1], tree, tmpdetached)
    _reinsert!(tree, tmpdetached)
    tree.nelems -= 1
    tree.nelem_deletions += 1
    return tree
end

Base.delete!(tree::RTree{T,N}, pt::Point{T,N}, id::Any = nothing) where {T,N} =
    delete!(tree, Rect(pt), id)

# deletes the subtree with the `node` root from the `tree`
# (does not update the parent MBR)
function delete_subtree!(tree::RTree, node::Node)
    #@debug "delete_subtree(): lev=$(level(node))"
    _release_descendants!(tree, node)

    if hasparent(node) # remove the node from the parent
        tree.nnodes_perlevel[level(node)] -= 1 # don't count this node anymore
        _detach!(parent(node), pos_in_parent(node), tree)
        release(tree, node)
        # FIXME propagate parent mbr updates (if required)
    else # root node just stays empty
        node.mbr = empty(mbrtype(node))
        empty!(node.children)
    end
end

# recursively release all `node` descendants to the pool
function _release_descendants!(tree::RTree, node::Node)
    if node isa Branch
        for child in children(node)
            _release_descendants!(tree, child)
            release(tree, child)
        end
        # update nodes counts
        tree.nnodes_perlevel[level(node)-1] -= length(node)
    elseif node isa Leaf
        tree.nelems -= length(node)
        tree.nelem_deletions += length(node)
    end
    #@debug "_release_descendants(): done node lev=$(level(node)) parent=$(hasparent(node))"
end

# FIXME replace `region` with a more generic `filter`
"""
    subtract!(tree::RTree, reg::Region)

Subtracts the `region` from the `tree`, i.e. removes all elements within
`region`.
"""
function subtract!(tree::RTree{T,N}, reg::Region{T,N}) where {T,N}
    #@debug "subtract!(): region=$(reg)"
    isempty(tree) && return tree

    tmpdetached = Vector{nodetype(tree)}()
    status = _subtract!(tree.root, 0, reg, tree, tmpdetached)
    #@debug "subtract!(): status=$(status) tmpdetached=$(length(tmpdetached))"
    if status > 0 # tree changed
        _condense!(tree.root, tree, tmpdetached) # try to condense the root
    end
    _reinsert!(tree, tmpdetached)
    return tree
end

function _subtract!(node::Node, node_ix::Int, reg::Region, tree::RTree,
                    tmpdetached::AbstractVector{<:Node})
    nodembr = mbr(node)
    #@debug "_subtract!(): lev=$(level(node)) i=$node_ix len=$(length(node))"
    # TODO juxtaposition() method that combines in() and intersects()?
    if in(nodembr, reg)
        #@debug "_subtract!(): delete subtree lev=$(level(node)) i=$node_ix len=$(length(node))"
        delete_subtree!(tree, node)
        return 2 # node removed
    elseif intersects(nodembr, reg)
        mbr_dirty = false
        i = 1
        while i <= length(node)
            #@debug "1: lev=$(level(node)) i=$i len=$(length(node))"
            child = node[i]
            oldmbr = mbr(child)
            if node isa Branch
                child_res = _subtract!(child, i, reg, tree, tmpdetached)
            else
                @assert node isa Leaf
                if in(oldmbr, reg)
                    #@debug "_subtract!(): detach elem lev=$(level(node)) i=$i len=$(length(node))"
                    _detach!(node, i, tree, updatembr = false) # don't update MBR, we do it later
                    tree.nelems -= 1
                    tree.nelem_deletions += 1
                    child_res = 2
                    # FIXME intersect is not considered (would be important if elem mbr is not a point), should be handled by `filter`
                else
                    child_res = 0
                end
            end
            if !mbr_dirty && (child_res != 0) && tree.tight_mbrs
                mbr_dirty = touches(nodembr, oldmbr)
            end
            if child_res != 2 # don't increment if i-th node removed
                i += 1
            end
        end
        if hasparent(node) && length(node) < floor(Int, tree.reinsert_factor * capacity(node, tree))
            _detach!(parent(node), node_ix, tree)
            tree.nnodes_perlevel[level(node)] -= 1
            if isempty(node)
                #@debug "Releasing empty node (lv=$(level(node)))"
                release(tree, node)
            else
                push!(tmpdetached, node)
            end
            return 2 # node removed
        elseif mbr_dirty
            syncmbr!(node)
            return 1 # mbr changed
        else
            return 0 # no mbr change
        end
    else
        return 0 # no change
    end
end


text:
push!(LOAD_PATH,"../src/")

using Documenter, EcologicalNetworksPlots

const pages = [
    "Index" => "index.md",
    "Layouts" => [
        "Introduction" => "layouts/initial.md",
        "Circular" => "layouts/circular.md",
        "Bipartite" => "layouts/bipartite.md",
        "Force-directed" => "layouts/forcedirected.md",
        "Unravel" => "layouts/unravelled.md"
    ],
    "Advanced topics" => [
        "Nodes attributes" => "advanced/attributes.md",
        "Networks subsets" => "advanced/subsets.md"
    ]
    # TODO add plotting, heatmap, advanced uses
]

makedocs(
    sitename = "EcologicalNetworksPlots",
    authors = "Timothée Poisot",
    pages = pages
    )

deploydocs(
    repo = "github.com/EcoJulia/EcologicalNetworksPlots.jl.git",
    push_preview = true
)


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

# deep copying

# Note: deepcopy_internal(::Any, ::ObjectIdDict) is
#       only exposed for specialization by libraries

deepcopy(x) = deepcopy_internal(x, ObjectIdDict())::typeof(x)

deepcopy_internal(x::Union{Symbol,Core.MethodInstance,Method,GlobalRef,DataType,Union,Task},
                  stackdict::ObjectIdDict) = x
deepcopy_internal(x::Tuple, stackdict::ObjectIdDict) =
    ntuple(i->deepcopy_internal(x[i], stackdict), length(x))
deepcopy_internal(x::Module, stackdict::ObjectIdDict) = error("deepcopy of Modules not supported")

function deepcopy_internal(x::SimpleVector, stackdict::ObjectIdDict)
    if haskey(stackdict, x)
        return stackdict[x]
    end
    y = Core.svec(Any[deepcopy_internal(x[i], stackdict) for i = 1:length(x)]...)
    stackdict[x] = y
    return y
end

function deepcopy_internal(x::String, stackdict::ObjectIdDict)
    if haskey(stackdict, x)
        return stackdict[x]
    end
    y = unsafe_string(pointer(x), sizeof(x))
    stackdict[x] = y
    return y
end

function deepcopy_internal(x::ANY, stackdict::ObjectIdDict)
    T = typeof(x)::DataType
    nf = nfields(T)
    (isbits(T) || nf == 0) && return x
    if haskey(stackdict, x)
        return stackdict[x]
    end
    y = ccall(:jl_new_struct_uninit, Any, (Any,), T)
    if T.mutable
        stackdict[x] = y
    end
    for i in 1:nf
        if isdefined(x,i)
            ccall(:jl_set_nth_field, Void, (Any, Csize_t, Any), y, i-1,
                  deepcopy_internal(getfield(x,i), stackdict))
        end
    end
    return y::T
end

function deepcopy_internal(x::Array, stackdict::ObjectIdDict)
    if haskey(stackdict, x)
        return stackdict[x]
    end
    _deepcopy_array_t(x, eltype(x), stackdict)
end

function _deepcopy_array_t(x::ANY, T, stackdict::ObjectIdDict)
    if isbits(T)
        return (stackdict[x]=copy(x))
    end
    dest = similar(x)
    stackdict[x] = dest
    for i = 1:(length(x)::Int)
        if ccall(:jl_array_isassigned, Cint, (Any, Csize_t), x, i-1) != 0
            xi = ccall(:jl_arrayref, Any, (Any, Csize_t), x, i-1)
            if !isbits(typeof(xi))
                xi = deepcopy_internal(xi, stackdict)
            end
            ccall(:jl_arrayset, Void, (Any, Any, Csize_t), dest, xi, i-1)
        end
    end
    return dest
end

function deepcopy_internal(x::Dict, stackdict::ObjectIdDict)
    if haskey(stackdict, x)
        return stackdict[x]::typeof(x)
    end

    if isbits(eltype(x))
        return (stackdict[x] = copy(x))
    end

    dest = similar(x)
    stackdict[x] = dest
    for (k, v) in x
        dest[deepcopy_internal(k, stackdict)] = deepcopy_internal(v, stackdict)
    end
    dest
end


text:
"""
    GroupedDataFrame

The result of a `groupby` operation on an AbstractDataFrame; a
view into the AbstractDataFrame grouped by rows.

Not meant to be constructed directly, see `groupby`.
"""
struct GroupedDataFrame{T<:AbstractDataFrame}
    parent::T
    cols::Vector{Int}    # columns used for grouping
    idx::Vector{Int}     # indexing vector when grouped by the given columns
    starts::Vector{Int}  # starts of groups
    ends::Vector{Int}    # ends of groups
end

#
# Split
#

"""
A view of an AbstractDataFrame split into row groups

```julia
groupby(d::AbstractDataFrame, cols; sort = false, skipmissing = false)
groupby(cols; sort = false, skipmissing = false)
```

### Arguments

* `d` : an AbstractDataFrame to split (optional, see [Returns](#returns))
* `cols` : data table columns to group by
* `sort`: whether to sort rows according to the values of the grouping columns `cols`
* `skipmissing`: whether to skip rows with `missing` values in one of the grouping columns `cols`

### Returns

A `GroupedDataFrame` : a grouped view into `d`

### Details

An iterator over a `GroupedDataFrame` returns a `SubDataFrame` view
for each grouping into `d`. A `GroupedDataFrame` also supports
indexing by groups, `map` (which applies a function to each group)
and `combine` (which applies a function to each group
and combines the result into a data frame).

See the following for additional split-apply-combine operations:

* `by` : split-apply-combine using functions
* `aggregate` : split-apply-combine; applies functions in the form of a cross product
* `colwise` : apply a function to each column in an `AbstractDataFrame` or `GroupedDataFrame`
* `map` : apply a function to each group of a `GroupedDataFrame` (without combining)
* `combine` : combine a `GroupedDataFrame`, optionally applying a function to each group

### Examples

```julia
df = DataFrame(a = repeat([1, 2, 3, 4], outer=[2]),
               b = repeat([2, 1], outer=[4]),
               c = randn(8))
gd = groupby(df, :a)
gd[1]
last(gd)
vcat([g[:b] for g in gd]...)
for g in gd
    println(g)
end
```

"""
function groupby(df::AbstractDataFrame, cols::AbstractVector;
                 sort::Bool = false, skipmissing::Bool = false)
    intcols = index(df)[cols]
    sdf = df[intcols]
    df_groups = group_rows(sdf, false, sort, skipmissing)
    GroupedDataFrame(df, intcols, df_groups.rperm,
                     df_groups.starts, df_groups.stops)
end
groupby(d::AbstractDataFrame, cols;
        sort::Bool = false, skipmissing::Bool = false) =
    groupby(d, [cols], sort = sort, skipmissing = skipmissing)

function Base.iterate(gd::GroupedDataFrame, i=1)
    if i > length(gd.starts)
        nothing
    else
        (view(gd.parent, gd.idx[gd.starts[i]:gd.ends[i]], :), i+1)
    end
end

Base.length(gd::GroupedDataFrame) = length(gd.starts)
Compat.lastindex(gd::GroupedDataFrame) = length(gd.starts)
Base.first(gd::GroupedDataFrame) = gd[1]
Base.last(gd::GroupedDataFrame) = gd[end]

Base.getindex(gd::GroupedDataFrame, idx::Integer) =
    view(gd.parent, gd.idx[gd.starts[idx]:gd.ends[idx]], :)
Base.getindex(gd::GroupedDataFrame, idxs::AbstractArray) =
    GroupedDataFrame(gd.parent, gd.cols, gd.idx, gd.starts[idxs], gd.ends[idxs])
Base.getindex(gd::GroupedDataFrame, idxs::Colon) =
    GroupedDataFrame(gd.parent, gd.cols, gd.idx, gd.starts, gd.ends)

function Base.:(==)(gd1::GroupedDataFrame, gd2::GroupedDataFrame)
    gd1.cols == gd2.cols &&
        length(gd1) == length(gd2) &&
        all(x -> ==(x...), zip(gd1, gd2))
end

function Base.isequal(gd1::GroupedDataFrame, gd2::GroupedDataFrame)
    isequal(gd1.cols, gd2.cols) &&
        isequal(length(gd1), length(gd2)) &&
        all(x -> isequal(x...), zip(gd1, gd2))
end

Base.names(gd::GroupedDataFrame) = names(gd.parent)
_names(gd::GroupedDataFrame) = _names(gd.parent)

"""
    map(cols => f, gd::GroupedDataFrame)
    map(f, gd::GroupedDataFrame)

Apply a function to each group of rows and return a `GroupedDataFrame`.

If the first argument is a `cols => f` pair, `cols` must be a column name or index, or
a vector or tuple thereof, and `f` must be a callable. If `cols` is a single column index,
`f` is called with a `SubArray` view into that column for each group; else, `f` is called
with a named tuple holding `SubArray` views into these columns.

If the first argument is a vector, tuple or named tuple of such pairs, each pair is
handled as described above. If a named tuple, field names are used to name
each generated column.

If the first argument is a callable, it is passed a `SubDataFrame` view for each group,
and the returned `DataFrame` then consists of the returned rows plus the grouping columns.
Note that this second form is much slower than the first one due to type instability.

`f` can return a single value, a row or multiple rows. The type of the returned value
determines the shape of the resulting data frame:
- A single value gives a data frame with a single column and one row per group.
- A named tuple of single values or a `DataFrameRow` gives a data frame with one column
  for each field and one row per group.
- A vector gives a data frame with a single column and as many rows
  for each group as the length of the returned vector for that group.
- A data frame, a named tuple of vectors or a matrix gives a data frame
  with the same columns and as many rows for each group as the rows returned for that group.

As a special case, if a tuple or vector of pairs is passed as the first argument, each function
is required to return a single value or vector, which will produce each a separate column.

In all cases, the resulting `GroupedDataFrame` contains all the grouping columns in addition
to those listed above. Column names are automatically generated when necessary: for functions
operating on a single column and returning a single value or vector, the function name is
appended to the input column name; for other functions, columns are called `x1`, `x2`
and so on.

Note that `f` must always return the same type of object for
all groups, and (if a named tuple or data frame) with the same fields or columns.
Due to type instability, returning a single value or a named tuple is dramatically
faster than returning a data frame.

### Examples

```jldoctest
julia> df = DataFrame(a = repeat([1, 2, 3, 4], outer=[2]),
                      b = repeat([2, 1], outer=[4]),
                      c = 1:8);

julia> gd = groupby(df, :a);

julia> map(:c => sum, gd)
GroupedDataFrame{DataFrame} with 4 groups based on key: :a
First Group: 1 row
│ Row │ a     │ c_sum │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 6     │
⋮
Last Group: 1 row
│ Row │ a     │ c_sum │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 4     │ 12    │

julia> map(df -> sum(df.c), gd) # Slower variant
GroupedDataFrame{DataFrame} with 4 groups based on key: :a
First Group: 1 row
│ Row │ a     │ x1    │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 6     │
⋮
Last Group: 1 row
│ Row │ a     │ x1    │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 4     │ 12    │
```

See [`by`](@ref) for more examples.

### See also

`combine(f, gd)` returns a `DataFrame` rather than a `GroupedDataFrame`

"""
function Base.map(f::Any, gd::GroupedDataFrame)
    if length(gd) > 0
        idx, valscat = _combine(f, gd)
        parent = hcat!(gd.parent[idx, gd.cols], valscat, makeunique=true)
        starts = Vector{Int}(undef, length(gd))
        ends = Vector{Int}(undef, length(gd))
        starts[1] = 1
        j = 2
        @inbounds for i in 2:length(idx)
            if idx[i] != idx[i-1]
                starts[j] = i
                ends[j-1] = i - 1
                j += 1
            end
        end
        # In case some groups have to be dropped
        resize!(starts, j-1)
        resize!(ends, j-1)
        ends[end] = length(idx)
        return GroupedDataFrame(parent, gd.cols, collect(1:length(idx)), starts, ends)
    else
        return GroupedDataFrame(parent, gd.cols, Int[], Int[], Int[])
    end
end

"""
    combine(gd::GroupedDataFrame)
    combine(gd::GroupedDataFrame, cols => f...)
    combine(gd::GroupedDataFrame; (colname = cols => f)...)
    combine(gd::GroupedDataFrame, f)
    combine(f, gd::GroupedDataFrame)

Transform a `GroupedDataFrame` into a `DataFrame`.

If the last argument(s) consist(s) in one or more `cols => f` pair(s), or if
`colname = cols => f` keyword arguments are provided, `cols` must be
a column name or index, or a vector or tuple thereof, and `f` must be a callable.
A pair or a (named) tuple of pairs can also be provided as the first or last argument.
If `cols` is a single column index, `f` is called with a `SubArray` view into that
column for each group; else, `f` is called with a named tuple holding `SubArray`
views into these columns.

If the last argument is a callable `f`, it is passed a `SubDataFrame` view for each group,
and the returned `DataFrame` then consists of the returned rows plus the grouping columns.
Note that this second form is much slower than the first one due to type instability.
A method is defined with `f` as the first argument, so do-block
notation can be used.

`f` can return a single value, a row or multiple rows. The type of the returned value
determines the shape of the resulting data frame:
- A single value gives a data frame with a single column and one row per group.
- A named tuple of single values or a `DataFrameRow` gives a data frame with one column
  for each field and one row per group.
- A vector gives a data frame with a single column and as many rows
  for each group as the length of the returned vector for that group.
- A data frame, a named tuple of vectors or a matrix gives a data frame
  with the same columns and as many rows for each group as the rows returned for that group.

As a special case, if a tuple or vector of pairs is passed as the first argument, each function
is required to return a single value or vector, which will produce each a separate column.

In all cases, the resulting data frame contains all the grouping columns in addition
to those listed above. Column names are automatically generated when necessary: for functions
operating on a single column and returning a single value or vector, the function name is
appended to the input column name; for other functions, columns are called `x1`, `x2`
and so on.

Note that `f` must always return the same type of object for
all groups, and (if a named tuple or data frame) with the same fields or columns.
Due to type instability, returning a single value or a named tuple is dramatically
faster than returning a data frame.

The resulting data frame will be sorted if `sort=true` was passed to the [`groupby`](@ref)
call from which `gd` was constructed. Otherwise, ordering of rows is undefined.

### Examples

```jldoctest
julia> df = DataFrame(a = repeat([1, 2, 3, 4], outer=[2]),
                      b = repeat([2, 1], outer=[4]),
                      c = 1:8);

julia> gd = groupby(df, :a);

julia> combine(gd, :c => sum)
4×2 DataFrame
│ Row │ a     │ c_sum │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 6     │
│ 2   │ 2     │ 8     │
│ 3   │ 3     │ 10    │
│ 4   │ 4     │ 12    │

julia> combine(:c => sum, gd)
4×2 DataFrame
│ Row │ a     │ c_sum │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 6     │
│ 2   │ 2     │ 8     │
│ 3   │ 3     │ 10    │
│ 4   │ 4     │ 12    │

julia> combine(df -> sum(df.c), gd) # Slower variant
4×2 DataFrame
│ Row │ a     │ x1    │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 6     │
│ 2   │ 2     │ 8     │
│ 3   │ 3     │ 10    │
│ 4   │ 4     │ 12    │
```

See [`by`](@ref) for more examples.

### See also

[`by(f, df, cols)`](@ref) is a shorthand for `combine(f, groupby(df, cols))`.

[`map`](@ref): `combine(f, groupby(df, cols))` is a more efficient equivalent
of `combine(map(f, groupby(df, cols)))`.

"""
function combine(f::Any, gd::GroupedDataFrame)
    if length(gd) > 0
        idx, valscat = _combine(f, gd)
        return hcat!(gd.parent[idx, gd.cols], valscat, makeunique=true)
    else
        return similar(gd.parent[gd.cols], 0)
    end
end
combine(gd::GroupedDataFrame, f::Any) = combine(f, gd)
combine(gd::GroupedDataFrame, f::Pair...) = combine(f, gd)
combine(gd::GroupedDataFrame; f...) =
    isempty(f) ? combine(identity, gd) : combine(values(f), gd)

# Wrapping automatically adds column names when the value returned
# by the user-provided function lacks them
wrap(x::Union{AbstractDataFrame, NamedTuple, DataFrameRow}) = x
wrap(x::AbstractMatrix) =
    NamedTuple{Tuple(gennames(size(x, 2)))}(Tuple(view(x, :, i) for i in 1:size(x, 2)))
wrap(x::Any) = (x1=x,)

function do_call(f::Any, gd::GroupedDataFrame, incols::AbstractVector, i::Integer)
    idx = gd.idx[gd.starts[i]:gd.ends[i]]
    f(view(incols, idx))
end
function do_call(f::Any, gd::GroupedDataFrame, incols::NamedTuple, i::Integer)
    idx = gd.idx[gd.starts[i]:gd.ends[i]]
    f(map(c -> view(c, idx), incols))
end
do_call(f::Any, gd::GroupedDataFrame, incols::Nothing, i::Integer) =
    f(gd[i])

_nrow(df::AbstractDataFrame) = nrow(df)
_nrow(x::NamedTuple{<:Any, <:Tuple{Vararg{AbstractVector}}}) =
    isempty(x) ? 0 : length(x[1])
_ncol(df::AbstractDataFrame) = ncol(df)
_ncol(x::Union{NamedTuple, DataFrameRow}) = length(x)

function gen_fun(::Type{NT}, f2) where NT
    get_val(::Val{T}) where {T} = T

    function(incols)
        tup = map(f2) do p
            nms = get_val(first(p))
            if nms isa Tuple
                res = last(p)(NamedTuple{nms}(map(c -> incols[c], nms)))
            else
                res = last(p)(incols[nms])
            end
            if res isa Union{AbstractDataFrame, NamedTuple, DataFrameRow, AbstractMatrix}
                throw(ArgumentError("a single value or vector result is required when passing " *
                                    "a vector or tuple of functions (got $(typeof(res)))"))
            end
            res
        end
        NT(tup)
    end
end

function _combine(f::Union{AbstractVector{<:Pair}, Tuple{Vararg{Pair}},
                           NamedTuple{<:Any, <:Tuple{Vararg{Pair}}}},
                  gd::GroupedDataFrame)
    idxs = (first(p) isa Union{Integer, Symbol} ?
            index(gd.parent)[first(p)] :
            index(gd.parent)[collect(first(p))] for p in f)
    incols = [names(gd.parent)[idx] for idx in idxs]
    allcols = collect(reduce(union, (nam isa Symbol ? (nam,) : nam for nam in incols)))
    f2 = ntuple(length(f)) do i
        nms = incols[i]
        (nms isa AbstractArray ? Val(Tuple(nms)) : Val(nms)) => last(f[i])
    end
    # Use temporary names for columns, and rename after the fact where appropriate
    NT = NamedTuple{Tuple(gennames(length(f2)))}
    fun = gen_fun(NT, f2)
    idx, valscat = _combine(allcols => fun, gd)
    if f isa NamedTuple
        nams = collect(propertynames(f))
    else
        nams = names(valscat)
        for i in 1:ncol(valscat)
            if f[i] isa Pair{<:Union{Symbol,Integer}}
                 nams[i] = Symbol(names(gd.parent)[index(gd.parent)[first(f[i])]],
                                  '_', funname(last(f[i])))
            end
        end
    end
    names!(valscat, nams)
    idx, valscat
end

function _combine(f::Any, gd::GroupedDataFrame)
    if f isa Pair{<:Union{Symbol,Integer}}
        incols = gd.parent[first(f)]
        fun = last(f)
    elseif f isa Pair
        df = gd.parent[collect(first(f))]
        incols = NamedTuple{Tuple(names(df))}(columns(df))
        fun = last(f)
    else
        incols = nothing
        fun = f
    end
    firstres = do_call(fun, gd, incols, 1)
    idx, valscat = _combine(wrap(firstres), fun, gd, incols)
    if f isa Pair{<:Union{Symbol,Integer}} &&
       !isa(firstres, Union{AbstractDataFrame, NamedTuple, DataFrameRow, AbstractMatrix})
        nam = Symbol(names(gd.parent)[index(gd.parent)[first(f)]],
                     '_', funname(fun))
        names!(valscat, [nam])
    end
    return idx, valscat
end

function _combine(first::Union{NamedTuple, DataFrameRow, AbstractDataFrame},
                  f::Any, gd::GroupedDataFrame,
                  incols::Union{Nothing, AbstractVector, NamedTuple})
    if first isa AbstractDataFrame
        n = 0
        eltys = eltypes(first)
    elseif first isa NamedTuple{<:Any, <:Tuple{Vararg{AbstractVector}}}
        n = 0
        eltys = map(eltype, first)
    elseif first isa DataFrameRow
        n = length(gd)
        eltys = eltypes(parent(first))
    else # NamedTuple giving a single row
        n = length(gd)
        eltys = map(typeof, first)
        if any(x -> x <: AbstractVector, eltys)
            throw(ArgumentError("mixing single values and vectors in a (named) tuple is not allowed"))
        end
    end
    idx = Vector{Int}(undef, n)
    initialcols = ntuple(i -> Tables.allocatecolumn(eltys[i], n), _ncol(first))
    outcols = _combine!(first, initialcols, idx, 1, 1, f, gd, incols,
                        tuple(propertynames(first)...))
    valscat = DataFrame(collect(outcols), collect(Symbol, propertynames(first)))
    idx, valscat
end

function fill_row!(row, outcols::NTuple{N, AbstractVector},
                   i::Integer, colstart::Integer,
                   colnames::NTuple{N, Symbol}) where N
    if !isa(row, Union{NamedTuple, DataFrameRow})
        throw(ArgumentError("return value must not change its kind " *
                            "(single row or variable number of rows) across groups"))
    elseif _ncol(row) != N
        throw(ArgumentError("return value must have the same number of columns " *
                            "for all groups (got $N and $(length(row)))"))
    end
    @inbounds for j in colstart:length(outcols)
        col = outcols[j]
        cn = colnames[j]
        local val
        try
            val = row[cn]
        catch
            throw(ArgumentError("return value must have the same column names " *
                                "for all groups (got $colnames and $(propertynames(row)))"))
        end
        S = typeof(val)
        T = eltype(col)
        if S <: T || promote_type(S, T) <: T
            col[i] = val
        else
            return j
        end
    end
    return nothing
end

function _combine!(first::Union{NamedTuple, DataFrameRow}, outcols::NTuple{N, AbstractVector},
                   idx::Vector{Int}, rowstart::Integer, colstart::Integer,
                   f::Any, gd::GroupedDataFrame,
                   incols::Union{Nothing, AbstractVector, NamedTuple},
                   colnames::NTuple{N, Symbol}) where N
    len = length(gd)
    # Handle first group
    j = fill_row!(first, outcols, rowstart, colstart, colnames)
    @assert j === nothing # eltype is guaranteed to match
    idx[rowstart] = gd.idx[gd.starts[rowstart]]
    # Handle remaining groups
    @inbounds for i in rowstart+1:len
        row = wrap(do_call(f, gd, incols, i))
        j = fill_row!(row, outcols, i, 1, colnames)
        if j !== nothing # Need to widen column type
            local newcols
            let i = i, j = j, outcols=outcols, row=row # Workaround for julia#15276
                newcols = ntuple(length(outcols)) do k
                    S = typeof(row[k])
                    T = eltype(outcols[k])
                    U = promote_type(S, T)
                    if S <: T || U <: T
                        outcols[k]
                    else
                        copyto!(Tables.allocatecolumn(U, length(outcols[k])),
                                1, outcols[k], 1, k >= j ? i-1 : i)
                    end
                end
            end
            return _combine!(row, newcols, idx, i, j, f, gd, incols, colnames)
        end
        idx[i] = gd.idx[gd.starts[i]]
    end
    outcols
end

# This needs to be in a separate function
# to work around a crash due to JuliaLang/julia#29430
if VERSION >= v"1.1.0-DEV.723"
    @inline function do_append!(do_it, col, vals)
        do_it && append!(col, vals)
        return do_it
    end
else
    @noinline function do_append!(do_it, col, vals)
        do_it && append!(col, vals)
        return do_it
    end
end

function append_rows!(rows, outcols::NTuple{N, AbstractVector},
                      colstart::Integer, colnames::NTuple{N, Symbol}) where N
    if !isa(rows, Union{AbstractDataFrame, NamedTuple{<:Any, <:Tuple{Vararg{AbstractVector}}}})
        throw(ArgumentError("return value must not change its kind " *
                            "(single row or variable number of rows) across groups"))
    elseif _ncol(rows) != N
        throw(ArgumentError("return value must have the same number of columns " *
                            "for all groups (got $N and $(_ncol(rows)))"))
    end
    @inbounds for j in colstart:length(outcols)
        col = outcols[j]
        cn = colnames[j]
        local vals
        try
            vals = rows[cn]
        catch
            throw(ArgumentError("return value must have the same column names " *
                                "for all groups (got $(Tuple(colnames)) and $(Tuple(names(rows))))"))
        end
        S = eltype(vals)
        T = eltype(col)
        if !do_append!(S <: T || promote_type(S, T) <: T, col, vals)
            return j
        end
    end
    return nothing
end

function _combine!(first::Union{AbstractDataFrame,
                                NamedTuple{<:Any, <:Tuple{Vararg{AbstractVector}}}},
                   outcols::NTuple{N, AbstractVector},
                   idx::Vector{Int}, rowstart::Integer, colstart::Integer,
                   f::Any, gd::GroupedDataFrame,
                   incols::Union{Nothing, AbstractVector, NamedTuple},
                   colnames::NTuple{N, Symbol}) where N
    len = length(gd)
    # Handle first group
    j = append_rows!(first, outcols, colstart, colnames)
    @assert j === nothing # eltype is guaranteed to match
    append!(idx, Iterators.repeated(gd.idx[gd.starts[rowstart]], _nrow(first)))
    # Handle remaining groups
    @inbounds for i in rowstart+1:len
        rows = wrap(do_call(f, gd, incols, i))
        j = append_rows!(rows, outcols, 1, colnames)
        if j !== nothing # Need to widen column type
            local newcols
            let i = i, j = j, outcols=outcols, rows=rows # Workaround for julia#15276
                newcols = ntuple(length(outcols)) do k
                    S = eltype(rows[k])
                    T = eltype(outcols[k])
                    U = promote_type(S, T)
                    if S <: T || U <: T
                        outcols[k]
                    else
                        copyto!(Tables.allocatecolumn(U, length(outcols[k])), outcols[k])
                    end
                end
            end
            return _combine!(rows, newcols, idx, i, j, f, gd, incols, colnames)
        end
        append!(idx, Iterators.repeated(gd.idx[gd.starts[i]], _nrow(rows)))
    end
    outcols
end

"""
Apply a function to each column in an AbstractDataFrame or
GroupedDataFrame

```julia
colwise(f, d)
```

### Arguments

* `f` : a function or vector of functions
* `d` : an AbstractDataFrame of GroupedDataFrame

### Returns

* various, depending on the call

### Examples

```julia
df = DataFrame(a = repeat([1, 2, 3, 4], outer=[2]),
               b = repeat([2, 1], outer=[4]),
               c = randn(8))
colwise(sum, df)
colwise([sum, length], df)
colwise((minimum, maximum), df)
colwise(sum, groupby(df, :a))
```

"""
colwise(f, d::AbstractDataFrame) = [f(d[i]) for i in 1:ncol(d)]

# apply several functions to each column in a data frame
colwise(fns::Union{AbstractVector, Tuple}, d::AbstractDataFrame) = [f(d[i]) for f in fns, i in 1:ncol(d)]
colwise(f, gd::GroupedDataFrame) = [colwise(f, g) for g in gd]

"""
    by(d::AbstractDataFrame, keys, cols => f...; sort::Bool = false)
    by(d::AbstractDataFrame, keys; (colname = cols => f)..., sort::Bool = false)
    by(d::AbstractDataFrame, keys, f; sort::Bool = false)
    by(f, d::AbstractDataFrame, keys; sort::Bool = false)

Split-apply-combine in one step: apply `f` to each grouping in `d`
based on grouping columns `keys`, and return a `DataFrame`.

`keys` can be either a single column index, or a vector thereof.

If the last argument(s) consist(s) in one or more `cols => f` pair(s), or if
`colname = cols => f` keyword arguments are provided, `cols` must be
a column name or index, or a vector or tuple thereof, and `f` must be a callable.
A pair or a (named) tuple of pairs can also be provided as the first or last argument.
If `cols` is a single column index, `f` is called with a `SubArray` view into that
column for each group; else, `f` is called with a named tuple holding `SubArray`
views into these columns.

If the last argument is a callable `f`, it is passed a `SubDataFrame` view for each group,
and the returned `DataFrame` then consists of the returned rows plus the grouping columns.
Note that this second form is much slower than the first one due to type instability.
A method is defined with `f` as the first argument, so do-block
notation can be used.

`f` can return a single value, a row or multiple rows. The type of the returned value
determines the shape of the resulting data frame:
- A single value gives a data frame with a single column and one row per group.
- A named tuple of single values or a `DataFrameRow` gives a data frame with one column
  for each field and one row per group.
- A vector gives a data frame with a single column and as many rows
  for each group as the length of the returned vector for that group.
- A data frame, a named tuple of vectors or a matrix gives a data frame
  with the same columns and as many rows for each group as the rows returned for that group.

As a special case, if multiple pairs are passed as last arguments, each function
is required to return a single value or vector, which will produce each a separate column.

In all cases, the resulting data frame contains all the grouping columns in addition
to those listed above. Column names are automatically generated when necessary: for functions
operating on a single column and returning a single value or vector, the function name is
appended to the input colummn name; for other functions, columns are called `x1`, `x2`
and so on.

Note that `f` must always return the same type of object for
all groups, and (if a named tuple or data frame) with the same fields or columns.
Due to type instability, returning a single value or a named tuple is dramatically
faster than returning a data frame.

The resulting data frame will be sorted on `keys` if `sort=true`.
Otherwise, ordering of rows is undefined.

`by(d, cols, f)` is equivalent to `combine(f, groupby(d, cols))` and to the
less efficient `combine(map(f, groupby(d, cols)))`.

### Examples

```jldoctest
julia> df = DataFrame(a = repeat([1, 2, 3, 4], outer=[2]),
                      b = repeat([2, 1], outer=[4]),
                      c = 1:8);

julia> by(df, :a, :c => sum)
4×2 DataFrame
│ Row │ a     │ c_sum │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 6     │
│ 2   │ 2     │ 8     │
│ 3   │ 3     │ 10    │
│ 4   │ 4     │ 12    │

julia> by(df, :a, d -> sum(d.c)) # Slower variant
4×2 DataFrame
│ Row │ a     │ x1    │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 6     │
│ 2   │ 2     │ 8     │
│ 3   │ 3     │ 10    │
│ 4   │ 4     │ 12    │

julia> by(df, :a) do d # do syntax for the slower variant
           sum(d.c)
       end
4×2 DataFrame
│ Row │ a     │ x1    │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 6     │
│ 2   │ 2     │ 8     │
│ 3   │ 3     │ 10    │
│ 4   │ 4     │ 12    │

julia> by(df, :a, :c => x -> 2 .* x)
8×2 DataFrame
│ Row │ a     │ c_function │
│     │ Int64 │ Int64      │
├─────┼───────┼────────────┤
│ 1   │ 1     │ 2          │
│ 2   │ 1     │ 10         │
│ 3   │ 2     │ 4          │
│ 4   │ 2     │ 12         │
│ 5   │ 3     │ 6          │
│ 6   │ 3     │ 14         │
│ 7   │ 4     │ 8          │
│ 8   │ 4     │ 16         │

julia> by(df, :a, c_sum = :c => sum, c_sum2 = :c => x -> sum(x.^2))
4×3 DataFrame
│ Row │ a     │ c_sum │ c_sum2 │
│     │ Int64 │ Int64 │ Int64  │
├─────┼───────┼───────┼────────┤
│ 1   │ 1     │ 6     │ 26     │
│ 2   │ 2     │ 8     │ 40     │
│ 3   │ 3     │ 10    │ 58     │
│ 4   │ 4     │ 12    │ 80     │

julia> by(df, :a, (:b, :c) => x -> (minb = minimum(x.b), sumc = sum(x.c)))
4×3 DataFrame
│ Row │ a     │ minb  │ sumc  │
│     │ Int64 │ Int64 │ Int64 │
├─────┼───────┼───────┼───────┤
│ 1   │ 1     │ 2     │ 6     │
│ 2   │ 2     │ 1     │ 8     │
│ 3   │ 3     │ 2     │ 10    │
│ 4   │ 4     │ 1     │ 12    │
```

"""
by(d::AbstractDataFrame, cols::Any, f::Any; sort::Bool = false) =
    combine(f, groupby(d, cols, sort = sort))
by(f::Any, d::AbstractDataFrame, cols::Any; sort::Bool = false) =
    by(d, cols, f, sort = sort)
by(d::AbstractDataFrame, cols::Any, f::Pair; sort::Bool = false) =
    combine(f, groupby(d, cols, sort = sort))
by(d::AbstractDataFrame, cols::Any, f::Pair...; sort::Bool = false) =
    combine(f, groupby(d, cols, sort = sort))
by(d::AbstractDataFrame, cols::Any; sort::Bool = false, f...) =
    combine(values(f), groupby(d, cols, sort = sort))

#
# Aggregate convenience functions
#

# Applies a set of functions over a DataFrame, in the from of a cross-product
"""
Split-apply-combine that applies a set of functions over columns of an
AbstractDataFrame or GroupedDataFrame

```julia
aggregate(d::AbstractDataFrame, cols, fs)
aggregate(gd::GroupedDataFrame, fs)
```

### Arguments

* `d` : an AbstractDataFrame
* `gd` : a GroupedDataFrame
* `cols` : a column indicator (Symbol, Int, Vector{Symbol}, etc.)
* `fs` : a function or vector of functions to be applied to vectors
  within groups; expects each argument to be a column vector

Each `fs` should return a value or vector. All returns must be the
same length.

### Returns

* `::DataFrame`

### Examples

```julia
using Statistics
df = DataFrame(a = repeat([1, 2, 3, 4], outer=[2]),
               b = repeat([2, 1], outer=[4]),
               c = randn(8))
aggregate(df, :a, sum)
aggregate(df, :a, [sum, x->mean(skipmissing(x))])
aggregate(groupby(df, :a), [sum, x->mean(skipmissing(x))])
```

"""
aggregate(d::AbstractDataFrame, fs::Any; sort::Bool=false) =
    aggregate(d, [fs], sort=sort)
function aggregate(d::AbstractDataFrame, fs::AbstractVector; sort::Bool=false)
    headers = _makeheaders(fs, _names(d))
    _aggregate(d, fs, headers, sort)
end

# Applies aggregate to non-key cols of each SubDataFrame of a GroupedDataFrame
aggregate(gd::GroupedDataFrame, f::Any; sort::Bool=false) = aggregate(gd, [f], sort=sort)
function aggregate(gd::GroupedDataFrame, fs::AbstractVector; sort::Bool=false)
    headers = _makeheaders(fs, setdiff(_names(gd), _names(gd.parent[gd.cols])))
    res = combine(x -> _aggregate(without(x, gd.cols), fs, headers), gd)
    sort && sort!(res, headers)
    res
end

# Groups DataFrame by cols before applying aggregate
function aggregate(d::AbstractDataFrame,
                   cols::Union{S, AbstractVector{S}},
                   fs::Any;
                   sort::Bool=false) where {S<:ColumnIndex}
    aggregate(groupby(d, cols, sort=sort), fs)
end

function funname(f)
    n = nameof(f)
    String(n)[1] == '#' ? :function : n
end

_makeheaders(fs::AbstractVector, cn::AbstractVector{Symbol}) =
    [Symbol(colname, '_', funname(f)) for f in fs for colname in cn]

function _aggregate(d::AbstractDataFrame, fs::AbstractVector,
                    headers::AbstractVector{Symbol}, sort::Bool=false)
    res = DataFrame(AbstractVector[vcat(f(d[i])) for f in fs for i in 1:size(d, 2)], headers, makeunique=true)
    sort && sort!(res, headers)
    res
end


text:
import Base.show, Base.zero
using Tsp.Util
type City
    id::Int
    x::Float32
    y::Float32
    tabu::Int
    chunk_x::Int
    chunk_y::Int
end

"""
A rectangular region of a 2D space

       |
    y1 |....+------+
       |    |      |
       |    |      |
    y0 |....+------+
       |    :      :
    (0,0)------------
            x0     x1

x, y - chunk coordinates on the grid
"""
type Chunk
    x0::Float32
    x1::Float32
    y0::Float32
    y1::Float32
    x::Int
    y::Int
    cities::Array{City}
end

type Grid
    cities::Array{City} # city data
    chunks::Array{Chunk, 2} # grid sectors
    distances::Array{Float32, 2} # doesn't make sense :O
end

type GridGeometry
    sw::Tuple{Float32,Float32}
    ne::Tuple{Float32,Float32}
    area::Float32
    width::Float32
    height::Float32
end

# =================================
#   CITY
# =================================

City(id::Int, x::Float32, y::Float32) = City(id, x, y, 0, -1, -1)

zero(::Type{City}) = City(0, NaN, NaN, -1, -1, -1)

function Base.show(stream::IO, city::City)
    show(stream, "City($(city.id) @ ($(city.x), $(city.y)), in chunk ($(city.chunk_x), $(city.chunk_y)))")
end

"""
    distance(City, City)
Calculate euclidean distance between 2 cities.
"""
@inline function distance(a::City, b::City)
    sqrt((a.x-b.x)^2 + (a.y-b.y)^2)
end

@inline function distance(a::City, b::City, grid::Grid)
    sqrt((a.x-b.x)^2 + (a.y-b.y)^2)
    # super slow for some reason
    # if a.id > b.id
    #     if grid.distances[a.id, b.id] == 0 && a.id != b.id
    #         grid.distances[a.id, b.id] = sqrt((a.x-b.x)^2 + (a.y-b.y)^2)
    #     end
    #     return grid.distances[a.id, b.id]
    # else
    #     if grid.distances[b.id, a.id] == 0 && a.id != b.id
    #         grid.distances[b.id, a.id] = sqrt((a.x-b.x)^2 + (a.y-b.y)^2)
    #     end
    #     return grid.distances[b.id, a.id]
    # end
end

"Make the city tabu until round n"
@inline function make_tabu(city::City, n::Int)
    city.tabu = n
end

"Load cities from file"
function load_cities(file::IOStream)
    size = file
    lines = eachline(file)
    size = parse(Int, first(lines))
    cities = Array(City, size)
    timeLimit = -1
    for line in lines
        data = split(line)
        if length(data) == 1
            timeLimit = parse(Int, data[1])
        else
            id = parse(Int, data[1])
            cities[id] = City(
                id,
                parse(Float32, data[2]),
                parse(Float32, data[3])
            )
        end
    end
    return cities, timeLimit
end

function load_cities(file::String)
    open(file) do f
        load_cities(f)
    end
end

function random_cities(n::Int, xs, ys)
    cities = Array(City, n)
    for i in 1:n
        cities[i] = City(i, Float32(rand(xs)), Float32(rand(ys)))
    end
    return cities
end

function random_cities(n::Int, filename::String, xs, ys)
    open(filename, "w") do f
        println(f, n)
        for i in 1:n
            c = City(i, Float32(rand(xs)), Float32(rand(ys)))
            println(f, "$(c.id) $(c.x) $(c.y)")
        end
    end
end

# =================================
#   CHUNK
# =================================

zero(::Type{Chunk}) = Chunk(NaN,NaN,NaN,NaN,-1,-1, zeros(City, 0))
Chunk(x0::Float32, x1::Float32, y0::Float32, y1::Float32, x::Int, y::Int) = Chunk(x0, x1, y0, y1, x, y, [])

Base.in(city::City, chunk::Chunk) = chunk.x0 <= city.x < chunk.x1 && chunk.y0 <= city.y < chunk.y1

function Base.show(io::IO, chunk::Chunk)
    show(io, "Chunk([$(chunk.x), $(chunk.y)] @ (($(chunk.x0), $(chunk.y0)),($(chunk.x1), $(chunk.y1))), $(length(chunk.cities)) cities)")
end

# =================================
#   GRID
# =================================

zero(::Type{Grid}) = Grid(zeros(City, 0), zeros(Chunk, (0,0)), zeros(Float32, (0,0)))

Grid(cities::Array{City}, chunks::Array{Chunk, 2}) = Grid(cities, chunks, zeros(Float32, (0,0)))

function Base.show(io::IO, grid::Grid)
    if length(grid.chunks) > 0
        x0 = grid.chunks[1].x0
        y0 = grid.chunks[1].y0
        x1 = grid.chunks[end].x1
        y1 = grid.chunks[end].y1
        show(io, "Grid((($x0, $y0), ($x1, $y1)), $(length(grid.cities)) cities, $(length(grid.chunks)) chunks)")
    else
        show(io, "Grid(<empty>)")
    end
end

"""
    make_grid(Array{City})
Place cities in chunks based on cooridnates
"""
function make_grid(cities::Array{City})
    grid_size = floor(Int, length(cities)^(1/5))
    xmin, xmax = map(v->v[1], find_extremes(map(c->c.x, cities)))
    ymin, ymax = map(v->v[1], find_extremes(map(c->c.y, cities)))
    xsteps = linspace(xmin-1, xmax+1, grid_size+1)
    ysteps = linspace(ymin-1, ymax+1, grid_size+1)
    chunks = Array(Chunk, (grid_size, grid_size))
    grid = Grid(cities, chunks)
    for yi in 1:length(ysteps)-1, xi in 1:length(xsteps)-1
        @inbounds grid.chunks[xi, yi] = Chunk(xsteps[xi], xsteps[xi+1], ysteps[yi], ysteps[yi+1], xi, yi)
    end
    # println("[make_grid] Created $(length(chunks)) chunks")
    _populate_chunks!(chunks, cities)
    # println("[make_grid] Chunks populated")
    if length(cities) < 100 # precalculate ditances for low number of cities.
        grid.distances = _calculate_distances(cities)
    end
    return grid, GridGeometry((xmin, ymin), (xmax, ymax))
end

function GridGeometry(sw::Tuple{Float32,Float32}, ne::Tuple{Float32,Float32})
    width = ne[1]-sw[1]
    height = ne[2]-sw[2]
    area = width * height
    GridGeometry(sw, ne, area, width, height)
end

function print_chunk_density(grid::Grid)
    for i in 1:length(grid.chunks[:, 1])
           for j in 1:length(grid.chunks[:,1])
               chunk = grid.chunks[i,j]
               area = (chunk.x1-chunk.x0)*(chunk.y1-chunk.y0)
               @printf("%8.5f | ",length(chunk.cities)/area)
           end
           println("\n",repeat("-",length(grid.chunks[:,1])*11-1))
       end
end

function _calculate_distances(cities)
    distances = spzeros(length(cities), length(cities))
    ctr = 0
    for i in 1:length(cities), j in 1:i
        # print("\r\33[2K$ctr distances done")
        @inbounds distances[i, j] = distance(cities[i], cities[j])
        ctr += 1
    end
    return distances
end

"""
Assign cities to sectors they are in
"""
function _populate_chunks!(chunks, cities)
    ctr = 1
    interval = 0.01 * length(cities)
    for city in cities, chunk in chunks
        if city in chunk
            push!(chunk.cities, city)
            city.chunk_x = chunk.x
            city.chunk_y = chunk.y
        end
        ctr += 1
    end
end


text:
using Random, Makie, MeshArrays, DataFrames, ColorSchemes, Statistics

"""
    PlotMakie(df::DataFrame,nn::Integer)

Plot random subset of size nn trajectories.
"""
function PlotMakie(df::DataFrame,nn::Integer,dMax::Float64=0.)
   IDs = randperm(maximum(df.ID))
   COs=[:gray76 :gold :limegreen :black]

   #scene=Scene(limits=FRect(0, 0, 40, 40),show_axis = false)
   #scene=Scene(limits=FRect(-185, -95, 370, 190),show_axis = false)
   scene=Scene(show_axis = false)
   df_by_ID = groupby(df, :ID)
   ID_list=[df_by_ID[i][1,:ID] for i in 1:length(df_by_ID)]
   for ii=1:nn
      jj=findall(ID_list.==IDs[ii])[1]
      tmp=df_by_ID[jj]
      if dMax > 0.
         d=abs.(diff(tmp[!,:lon]))
         jj=findall(d .> dMax)
         tmp[jj,:lon].=NaN; tmp[jj,:lat].=NaN
         d=abs.(diff(tmp[!,:lat]))
         jj=findall(d .> dMax)
         tmp[jj,:lon].=NaN; tmp[jj,:lat].=NaN
      end
      CO=COs[mod(ii,4)+1]
      Makie.lines!(scene,tmp[!,:lon],tmp[!,:lat],color=CO,linewidth=1.0)
   end

   return scene
end

"""
MakieScatterMovie(scene,df,tt,fil::String)

Animate positions, according to time vector tt, and save movie to mp4 file.

```
using IndividualDisplacements
p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/recipes_Makie.jl"))
include(joinpath(p,"../examples/helper_functions.jl"))
module ex3
    fil="../examples/worldwide/three_dimensional_ocean.jl"
    include(joinpath(Main.p,fil))
    export set_up_individuals, Γ, 𝐷
end

using .ex3
𝐼=set_up_individuals(ex3.𝐼,nf=10000);
𝑇=(0.0,2*𝐼.𝑃.𝑇[2])
∫!(𝐼,𝑇)

θ=0.5*(𝐷.θ0+𝐷.θ1)
scene = MakieBase(θ,2:2:28,Tiso=15)
MakieScatterMovie(scene,𝐼.🔴,0:0.02:4,"tmp.mp4")
```
"""
function MakieScatterMovie(scene::Scene,df,tt,fil::String)

   🔴_by_t = groupby(𝐼.🔴, :t)
   _, threeD, twoD = MakieScatter(scene,🔴_by_t[end])

   np,nmax=length(🔴_by_t[end][:lon]),100000
   (xs,ys,zs)=fill(NaN,nmax),fill(NaN,nmax),fill(NaN,nmax)
   cs=fill(:blue,nmax)
   zmul=1/5

   ye=[🔴_by_t[i][1,:year] for i in 1:length(🔴_by_t)]
   tt,dt=collect(tt),median(diff(tt))

   scene
   record(scene, fil, 1:length(tt); framerate=12) do i
       jj = findall( (ye.>tt[i]-dt/2).&(ye.<=tt[i]+dt/2) )
       [xs[collect((1:np).+(j-jj[1])*np)]=🔴_by_t[j][:,:lon] for j in jj]
       [ys[collect((1:np).+(j-jj[1])*np)]=🔴_by_t[j][:,:lat] for j in jj]
       [zs[collect((1:np).+(j-jj[1])*np)]=🔴_by_t[j][:,:z] for j in jj]
       [cs[collect((1:np).+(j-jj[1])*np)]=🔴_by_t[j][:,:col] for j in jj]

       #xs[xs.>180.0] .-= 360.0
       xs[xs.<20.0] .+= 360.0
       #
       threeD[1] = xs
       threeD[2] = ys
       threeD[3] = zmul*zs
       threeD[:color] = zmul*zs
       #
       twoD[1] = xs
       twoD[2] = ys
       twoD[:color] = cs
       #
       rotate_cam!(scene,-0.05, 0., 0.)
   end

   return scene
end

"""
    MakieScatter(scene::Scene,df)

Add a scatter plot of e.g. x,y,z

```
scene = MakieBase(θ,2:2:28;Tiso=15)
_, threeD, twoD = MakieScatter(scene,🔴_by_t[end])
```
"""
function MakieScatter(scene,df)

    zmul=1/5
    nmax=100000

    xs=fill(NaN,nmax)
    ys=fill(NaN,nmax)
    zs=fill(NaN,nmax)
    cs=fill(:black,nmax)
    nt=length(df[!, :lon])
    xs[1:nt] = deepcopy(df[!, :lon])
    #xs[xs.>180.0] .-= 360.0
    xs[xs.<20.0] .+= 360.0
    ys[1:nt] = deepcopy(df[!, :lat])
    zs[1:nt] = deepcopy(df[!, :z])
    z0=0*zs .- 200.0
    cs[1:nt]=deepcopy(df[!, :col])

    Makie.scatter!(scene, xs, ys, zmul*zs, markersize = 500.0,
    show_axis = false, color=zs, colormap = :turbo, strokewidth=0.0)
    threeD = scene[end]

    Makie.scatter!(scene, xs, ys, zmul*z0, markersize = 500.0,
    show_axis = false, color=cs, strokewidth=0.0)
    twoD = scene[end]

    return scene, threeD, twoD
end

"""
    MakieBase(θ,T; LONin=140.:0.5:250.,LATin=10.:0.5:50.,DEPin=0.:10:200.)

Contour plot of a gridded 2D array projected onto e.g. 200m depth plane.

```
scene = MakieBase(θ,2:2:28)
```
"""
function MakieBase(θ,T; Tiso=12, LONin=140.:0.5:250.,LATin=10.:0.5:50.,DEPin=0.:10:200.)

    sds=(x1=false,xn=true,y1=true,yn=false)
    zmul=1/5

    isa(T,AbstractRange) ? T=collect(T) : nothing

    lo=[i for i=LONin, j=LATin]
    la=[j for i=LONin, j=LATin]
    lon=deepcopy(lo); lat=deepcopy(la);

    xlims=extrema(Γ["XC"][1])
    lo[findall(lo.<xlims[1])]=lo[findall(lo.<xlims[1])].+360
    lo[findall(lo.>xlims[2])]=lo[findall(lo.>xlims[2])].-360
    (f,i,j,w,_,_,_)=InterpolationFactors(Γ,vec(lon),vec(lat))
    IntFac=(lon=lon,lat=lat,f=f,i=i,j=j,w=w)

    dMin,dMax=extrema(DEPin)
    d=isosurface(θ,Tiso,Γ["RC"])
    d[1][findall(isnan.(d[1]))].=0.
    dd=interp_to_lonlat(d,IntFac)
    dd[findall(dd.>=-dMin)].=NaN #-dMin
    dd[findall(dd.<=-dMax)].=NaN #-dMax

    kMax=maximum(findall(Γ["RC"].>-dMax))
    θbox=fill(NaN,(size(lo)...,kMax));
    [θbox[:,:,k].=interp_to_lonlat(θ[:,k],IntFac) for k=1:kMax];

    lim=FRect3D([minimum(lon) minimum(lat) -dMax*zmul],
            [maximum(lon)-minimum(lon) maximum(lat)-minimum(lat) (dMax-dMin)*zmul])
    #eyepos=Vec3f0(minimum(lon),mean(lat),-1.1*dMax*zmul)
    #lookat=Vec3f0(maximum(lon)+60,mean(lat)-40,-0.1*dMax*zmul)

    scene=Scene(camera = cam3d!,colormap = :blues)

    #bottom
    Makie.contour!(scene,vec(lon[:,1]), vec(lat[1,:]), θbox[:,:,kMax], levels = T,
    transformation = (:xy, -dMax*zmul), limits=lim, color=:black, linewidth = 2)#,show_axis = false)

    #sides
    sds.x1 ? Makie.contour!(scene,vec(lat[1,:]),vec(Γ["RC"][1:kMax])*zmul,θbox[1,:,:],
    levels=T, transformation = (:yz, lon[1,1]), color=:black, linewidth = 2) : nothing

    sds.xn ? Makie.contour!(scene,vec(lat[end,:]),vec(Γ["RC"][1:kMax])*zmul,θbox[end,:,:],
    levels=T, transformation = (:yz, lon[end,1]), color=:black, linewidth = 2) : nothing

    sds.y1 ? Makie.contour!(scene,vec(lon[:,1]),vec(Γ["RC"][1:kMax])*zmul,θbox[:,1,:],
    levels=T, transformation = (:xz, lat[1,1]), color=:black, linewidth = 2) : nothing

    sds.yn ? Makie.contour!(scene,vec(lon[:,end]),vec(Γ["RC"][1:kMax])*zmul,θbox[:,end,:],
    levels=T, transformation = (:xz, lat[1,end]), color=:black, linewidth = 2) : nothing

    #isotherm
    bw(i::Float64) = ColorSchemes.RGB{Float64}(i, i, i)
    wireframe!(scene,vec(lon[:,1]),vec(lat[1,:]), zmul*dd,
    linewidth=0.25, color=bw(0.7))

    #xlabel!("lon"); ylabel!("lat"); zlabel!("$zmul x depth")
    #xlabel!(""); ylabel!(""); zlabel!("")

    return scene
end


text:
export iauTf2d
"""
**  Convert hours, minutes, seconds to days.
**
**  This function is part of the International Astronomical Union's
**  SOFA (Standards of Fundamental Astronomy) software collection.
**
**  Status:  support function.
**
**  Given:
**     s         char    sign:  '-' = negative, otherwise positive
**     ihour     int     hours
**     imin      int     minutes
**     sec       double  seconds
**
**  Returned:
**     days      double  interval in days
**
**  Returned (function value):
**               int     status:  0 = OK
**                                1 = ihour outside range 0-23
**                                2 = imin outside range 0-59
**                                3 = sec outside range 0-59.999...
**
**  Notes:
**
**  1)  The result is computed even if any of the range checks fail.
**
**  2)  Negative ihour, imin and/or sec produce a warning status, but
**      the absolute value is used in the conversion.
**
**  3)  If there are multiple errors, the status value reflects only the
**      first, the smallest taking precedence.
**
**  This revision:  2013 June 18
**
**  SOFA release 2018-01-30
**
**  Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# int iauTf2d(char s, int ihour, int imin, double sec, double *days)
function iauTf2d(s::Char, ihour::Int, imin::Int, sec::Real)
    ref_days  = Ref{Float64}(0.0)

    status = ccall((:iauTf2d, libsofa_c), Cint,
        (Cchar, Cint, Cint, Cdouble, Ref{Cdouble}),
        convert(UInt8, s),
        convert(Int32, ihour),
        convert(Int32, imin),
        convert(Float64, sec),
        ref_days)

    return status, ref_days[]
end


text:
#### Tendency types

# Terminology:
#
# `∂_t Yᵢ + (∇•F₁(Y))ᵢ + (∇•F₂(Y,G)))ᵢ = (S(Y,G))ᵢ`
# `__Tᵤ__   ____T₁____   ______T₂______   ___S___`
#
#  - `Yᵢ` - the i-th prognostic variable
#  - `Y` - the prognostic state (column vector)
#  - `G = ∇Y` - the gradient of the prognostic state (rank 2 tensor)
#  - `F₁` - the first order tendency flux (rank 2 tensor)
#  - `F₂` - the second order tendency flux (rank 2 tensor)
#  - `Tᵤ` - the explicit time derivative (column vector)
#  - `T₁` - the first order flux divergence (column vector)
#  - `T₂` - the second order flux divergence (column vector)
#  - `S` - the non-conservative source (column vector)

using DispatchedTuples

export AbstractPrognosticVariable,
    AbstractMomentumVariable,
    AbstractEnergyVariable,
    AbstractMoistureVariable,
    AbstractPrecipitationVariable,
    AbstractTracersVariable

export FirstOrder, SecondOrder
export AbstractTendencyType, Flux, Source
export TendencyDef
export prognostic_var_source_map
export eq_tends, prognostic_vars

"""
    AbstractPrognosticVariable

Subtypes are used for specifying
each prognostic variable.
"""
abstract type AbstractPrognosticVariable end

abstract type AbstractMomentumVariable <: AbstractPrognosticVariable end
abstract type AbstractEnergyVariable <: AbstractPrognosticVariable end
abstract type AbstractMoistureVariable <: AbstractPrognosticVariable end
abstract type AbstractPrecipitationVariable <: AbstractPrognosticVariable end
abstract type AbstractTracersVariable{N} <: AbstractPrognosticVariable end

"""
    AbstractOrder

Subtypes are used for dispatching
on the flux order.
"""
abstract type AbstractOrder end

"""
    FirstOrder

A type for dispatching on first order fluxes
"""
struct FirstOrder <: AbstractOrder end

"""
    SecondOrder

A type for dispatching on second order fluxes
"""
struct SecondOrder <: AbstractOrder end

"""
    AbstractTendencyType

Subtypes are used for specifying a
tuple of tendencies to be accumulated.
"""
abstract type AbstractTendencyType end

"""
    Flux{O <: AbstractOrder}

A type for dispatching on flux tendency types
where `O` is an abstract order ([`FirstOrder`](@ref)
or [`SecondOrder`](@ref)).
"""
struct Flux{O <: AbstractOrder} <: AbstractTendencyType end

"""
    Source

A type for dispatching on source tendency types
"""
struct Source <: AbstractTendencyType end

"""
    TendencyDef

Subtypes are used for specifying
each tendency definition.
"""
abstract type TendencyDef{TT <: AbstractTendencyType} end

"""
    eq_tends(::AbstractPrognosticVariable, ::BalanceLaw, ::AbstractTendencyType)

A tuple of `TendencyDef`s given
 - `AbstractPrognosticVariable` prognostic variable
 - `AbstractTendencyType` tendency type
 - `BalanceLaw` balance law

i.e., a tuple of `TendencyDef`s corresponding
to `F₁`, `F₂`, **or** `S` for a single
prognostic variable in:

    `∂_t Yᵢ + (∇•F₁(Y))ᵢ + (∇•F₂(Y,G)))ᵢ = (S(Y,G))ᵢ`
"""
function eq_tends end

"""
    prognostic_vars(::BalanceLaw)

A tuple of `AbstractPrognosticVariable`s given
the `BalanceLaw`.

i.e., a tuple of `AbstractPrognosticVariable`s
corresponding to the column-vector `Yᵢ` in:

    `∂_t Yᵢ + (∇•F₁(Y))ᵢ + (∇•F₂(Y,G)))ᵢ = (S(Y,G))ᵢ`
"""
prognostic_vars(::BalanceLaw) = ()

"""
    projection(bl, ::TendencyDef, args, x)

Provide a hook to project individual tendencies.
Return identity by defualt
"""
projection(pv::PV, bl, ::TendencyDef{TT}, args, x) where {TT, PV} = x

"""
    var, name = get_prog_state(state::Union{Vars, Grad}, pv::AbstractPrognosticVariable)

Returns a tuple of two elements. `var` is a `Vars` or `Grad`
object, and `name` is a Symbol. They should be linked such that
`getproperty(var, name)` returns the corresponding prognostic
variable type `pv`.

# Example

```julia
get_prog_state(state, ::TotalMoisture) = (state.moisture, :ρq_tot)
var, name = get_prog_state(state, TotalMoisture())
@test getproperty(var, name) == state.moisture.ρq_tot
```
"""
function get_prog_state end

"""
    precompute(bl, args, ::AbstractTendencyType)

A nested NamedTuple of precomputed (cached) values
and or objects. This is useful for "expensive"
point-wise quantities that are used in multiple
tendency terms. For example, computing a quantity
that requires iteration.
"""
precompute(bl, args, ::AbstractTendencyType) = NamedTuple()

"""
    prognostic_var_source_map(driver_sources::Tuple)

A DispatchedTuple, given a Tuple
of the driver/experiment sources.

!!! note
    `prognostic_vars`, which returns a Tuple
    of prognostic variable types, must be
    defined for boundary condition types.
"""
function prognostic_var_source_map(driver_sources::Tuple)
    tup = map(driver_sources) do t
        map(prognostic_vars(t)) do pv
            (pv, t)
        end
    end
    tup = tuple_of_tuples(tup)
    return DispatchedTuple(tup)
end

# Flatten "tuple of tuple of tuples" to "tuple of tuples"
tuple_of_tuples(a::Tuple{AbstractPrognosticVariable, T}) where {T} = (a,)
tuple_of_tuples(a, b...) =
    tuple(tuple_of_tuples(a)..., tuple_of_tuples(b...)...)
tuple_of_tuples(a::Tuple) = tuple_of_tuples(a...)
tuple_of_tuples(a::Tuple{}) = a


text:
function KronGraph500NoPerm(SCALE,EdgePerVertex)
    N = 2^SCALE

    M = round(EdgePerVertex * N)

    A = 0.57
    B = 0.19
    C = 0.19
    D = 1-(A+B+C)

    ij = Array{Int}(ones(2,M))
    ab = A+B
    c_norm = C/(1-(A+B))
    a_norm = A/(A+B)

    for ib = 1:SCALE
        ii_bit = rand(1,M) .> ab
        jj_bit = rand(1,M) .> (c_norm * ii_bit + a_norm * !(ii_bit))
        ij = ij + 2^(ib-1) * [ii_bit;jj_bit]
    end
    return reshape(ij[1,:]',size(ij,2)), reshape(ij[2,:]',size(ij,2))
end


text:
@safetestset "Assets rendering" begin

  @safetestset "Embedded assets" begin
    using Genie
    using Genie.Renderer
    using Genie.Assets

    @test (Assets.js_settings() * Assets.embedded(joinpath("files", "embedded", "channels.js"))) == Assets.channels()

    @test Assets.channels()[1:18] == "window.Genie = {};"

    @test Assets.channels_script()[1:28] == "<script>\nwindow.Genie = {};\n"
  end;

end;


text:
"""
Subroutine dgpstep

This subroutine computes the gradient projection step

  s = P[x + alpha*w] - x,

where P is the projection on the n-dimensional interval [xl,xu].
"""
function dgpstep(n,x,xl,xu,alpha,w,s)
    # TODO
    # This computation of the gradient projection step avoids
    # rouding errors for the components that are feasible.

    @inbounds for i=1:n
        if x[i] + alpha*w[i] < xl[i]
            s[i] = xl[i] - x[i]
        elseif x[i] + alpha*w[i] > xu[i]
            s[i] = xu[i] - x[i]
        else
            s[i] = alpha*w[i]
        end
    end
    return
end


text:
#=
Base.@kwdef mutable struct SIPOptimizer <: MOI.AbstractOptimizer
    objective_sense::MOI.OptimizationSense = MOI.FEASIBLE_SENSE
    time_limit::Float64 = 1000.0
    is_silent::Bool = false
    x::Vector{Float64} = Float64[]
    obj_bound::Float64 = 0.0
    obj_value::Float64 = 0.0
    termination_status::MOI.TerminationStatusCode = MOI.OPTIMIZE_NOT_CALLED
    result_status_code::MOI.ResultStatusCode = MOI.UNKNOWN_RESULT_STATUS
    run_time::Float64 = 0.0
end

function MOI.empty!(m::SIPOptimizer)
end
function MOI.is_empty(m::Optimizer)
end

function MOI.set(m::SIPOptimizer, ::MOI.Silent, value)
    m.is_silent = value
    return nothing
end

function MOI.set(m::SIPOptimizer, ::MOI.TimeLimitSec, value::Nothing)
    m.time_limit = Inf
    return nothing
end
function MOI.set(m::SIPOptimizer, ::MOI.TimeLimitSec, value::Float64)
    m.time_limit = value
    return nothing
end

MOI.set(m::SIPOptimizer, p::MOI.RawParameter, value) = error("The SIPOptimizer does not have any raw parameters.")

function MOI.get(m::SIPOptimizer, ::MOI.ListOfVariableIndices)
end
function MOI.get(m::SIPOptimizer, ::MOI.NumberOfVariables)
end

function MOI.get(m::SIPOptimizer, ::MOI.ObjectiveValue)

end
function MOI.get(m::SIPOptimizer, ::MOI.ObjectiveBound)

end

MOI.get(m::SIPOptimizer, p::MOI.RawParameter) = error("The SIPOptimizer does not have any raw parameters.")

function MOI.get(m::SIPOptimizer, ::MOI.RelativeGap)
    abs(objective_bound() - objective_value())/objective_value()
end
MOI.get(m::SIPOptimizer, ::MOI.SolverName) = "EAGO: SIPOptimizer"
MOI.get(m::SIPOptimizer, ::MOI.TerminationStatus) = m._termination_status_code
MOI.get(m::SIPOptimizer, ::MOI.PrimalStatus) = m.result_status_code
MOI.get(m::SIPOptimizer, ::MOI.SolveTime) = m.run_time
MOI.get(m::SIPOptimizer, ::MOI.ResultCount) = (m.result_status_code === MOI.FEASIBLE_POINT) ? 1 : 0
MOI.get(m::SIPOptimizer, ::MOI.TimeLimitSec) = m.time_limit

MOI.supports(::SIPOptimizer, ::MOI.TimeLimitSec) = true
MOI.supports(::SIPOptimizer, ::MOI.ObjectiveSense) = true
function MOI.set(m::SIPOptimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    m.objective_sense = sense
    return nothing
end

function MOI.optimize!(m::SIPOptimizer)
end
=#


text:
module test_julia_strings

using Test

s = first(split("a,b", ","))

@test s isa AbstractString
@test s isa SubString{String}
@test !(s isa String)

@test repeat("123", 4) == "123123123123"
@test reshape(collect(repeat("123", 4)), (6,2)) ==
['1'  '1'
 '2'  '2'
 '3'  '3'
 '1'  '1'
 '2'  '2'
 '3'  '3']

@test rotl90(reshape(collect(repeat("123", 4)), (6,2))) ==
['1'  '2'  '3'  '1'  '2'  '3'
 '1'  '2'  '3'  '1'  '2'  '3']

@test occursin("and", "candy")
@test contains("candy", "and")
if VERSION >= v"1.6"
@test occursin("candy")("and")
@test contains("and")("candy")
end

buf = Base.SecretBuffer("hello")
@test read(buf, String) == "hello"
@test read(buf, String) == ""

seekstart(buf)
@test read(buf, String) == "hello"

Base.shred!(buf)
@test Base.isshredded(buf)

seekstart(buf)
@test read(buf, String) == ""

@test length("가") == 1
@test collect("가") == ['가']

@test ncodeunits("가") == 3
@test codeunits("가")    == b"가"    == Vector{UInt8}("가")    == [0xea, 0xb0, 0x80]
@test codeunits("apple") == b"apple" == Vector{UInt8}("apple") == [0x61, 0x70, 0x70, 0x6c, 0x65]
@test collect("apple") == ['a', 'p', 'p', 'l', 'e']

if VERSION >= v"1.3"
    x = pi
    @test repr(x; context=nothing) == sprint(show, x; context=nothing) == "π"

    s = IOBuffer()
    context = :compact => true
    ctx = IOContext(s, context)
    show(ctx, pi)
    @test String(take!(s)) == "π"
end

@test string(12345, base = 16) == "3039"
@test hex2bytes(string(12345, base = 16)) == [0x30, 0x39]
@test 12345 == 0x3039 == Int(0x3039) == (0x3000 + 0x0039) == (UInt16(0x30) << 8 + 0x39)
@test bytes2hex([0x30, 0x39]) == "3039"

end # module test_julia_strings


text:
using OperatorLearning, Flux

versioninfo()
include("bench_utils.jl")

@info "Benchmark FourierLayer"
include("benchFourierLayer.jl")


text:
using Distributions, LinearAlgebra, Compat, Expectations

# A default utility function

function u(c, σ)
    if c > 0
        return (c^(1 - σ) - 1) / (1 - σ)
    else
        return -10e6
    end
end

# default wage vector with probabilities

const n = 60                                           # n possible outcomes for wage
const default_w_vec = range(10, 20, length = n) # wages between 10 and 20
const a, b = 600, 400                                  # shape parameters
const dist = BetaBinomial(n-1, a, b)

mutable struct McCallModel{TF <: AbstractFloat,
                           TAV <: AbstractVector{TF}}
    α::TF         # Job separation rate
    β::TF         # Discount rate
    γ::TF         # Job offer rate
    c::TF         # Unemployment compensation
    σ::TF         # Utility parameter
    w_vec::TAV    # Possible wage values

    McCallModel(α::TF = 0.2,
                β::TF = 0.98,
                γ::TF = 0.7,
                c::TF = 6.0,
                σ::TF = 2.0,
                w_vec::TAV = default_w_vec,
                ) where {TF, TAV} =
        new{TF, TAV}(α, β, γ, c, σ, w_vec)
end

function update_bellman!(mcm, V, V_new, U, E)
    # Simplify notation
    α, β, σ, c, γ = mcm.α, mcm.β, mcm.σ, mcm.c, mcm.γ

    for (w_idx, w) in enumerate(mcm.w_vec)
        # w_idx indexes the vector of possible wages
        V_new[w_idx] = u(w, σ) + β * ((1 - α) * V[w_idx] + α * U)
    end

    U_new = u(c, σ) + β * (1 - γ) * U +
            β * γ * E*max.(U, V)
    return U_new
end

function solve_mccall_model(mcm; tol = 1e-5, max_iter = 2000)

    V = ones(length(mcm.w_vec))    # Initial guess of V
    V_new = similar(V)             # To store updates to V
    U = 1.0                        # Initial guess of U
    i = 0
    error = tol + 1
    E = expectation(dist, nodes = mcm.w_vec)

    while error > tol && i < max_iter
        U_new = update_bellman!(mcm, V, V_new, U, E)
        error_1 = maximum(abs, V_new - V)
        error_2 = abs(U_new - U)
        error = max(error_1, error_2)
        V[:] = V_new
        U = U_new
        i += 1
    end

    return V, U
end

max_stars_repo_name: yywingliang/lecture-source-jl

max_issues_repo_name: yywingliang/lecture-source-jl

max_forks_repo_name: yywingliang/lecture-source-jl


text:
module Read

import ..Git, ..Cache, ..Reqs
using ..Types

readstrip(path...) = strip(readall(joinpath(path...)))

url(pkg::String) = readstrip("METADATA", pkg, "url")
sha1(pkg::String, ver::VersionNumber) = readstrip("METADATA", pkg, "versions", string(ver), "sha1")

function available(names=readdir("METADATA"))
    pkgs = Dict{ByteString,Dict{VersionNumber,Available}}()
    for pkg in names
        isfile("METADATA", pkg, "url") || continue
        versdir = joinpath("METADATA", pkg, "versions")
        isdir(versdir) || continue
        for ver in readdir(versdir)
            ismatch(Base.VERSION_REGEX, ver) || continue
            isfile(versdir, ver, "sha1") || continue
            haskey(pkgs,pkg) || (pkgs[pkg] = eltype(pkgs)[2]())
            pkgs[pkg][convert(VersionNumber,ver)] = Available(
                readchomp(joinpath(versdir,ver,"sha1")),
                Reqs.parse(joinpath(versdir,ver,"requires"))
            )
        end
    end
    return pkgs
end
available(pkg::String) = get(available([pkg]),pkg,Dict{VersionNumber,Available}())

isinstalled(pkg::String) =
    pkg != "METADATA" && pkg != "REQUIRE" && pkg[1] != '.' && isdir(pkg)

function isfixed(pkg::String, avail::Dict=available(pkg))
    isinstalled(pkg) || error("$pkg is not an installed package.")
    isfile("METADATA", pkg, "url") || return true
    ispath(pkg, ".git") || return true
    Git.dirty(dir=pkg) && return true
    Git.attached(dir=pkg) && return true
    !Git.success(`cat-file -e HEAD:REQUIRE`, dir=pkg) && isfile(pkg,"REQUIRE") && return true
    head = Git.head(dir=pkg)
    for (ver,info) in avail
        head == info.sha1 && return false
    end
    cache = Cache.path(pkg)
    cache_has_head = isdir(cache) && Git.iscommit(head, dir=cache)
    for (ver,info) in avail
        if cache_has_head && Git.iscommit(info.sha1, dir=cache)
            Git.is_ancestor_of(head, info.sha1, dir=cache) && return false
        elseif Git.iscommit(info.sha1, dir=pkg)
            Git.is_ancestor_of(head, info.sha1, dir=pkg) && return false
        else
            Base.warn_once("unknown $pkg commit $(info.sha1[1:8]), metadata may be ahead of package cache")
        end
    end
    return true
end

function installed_version(pkg::String, avail::Dict=available(pkg))
    ispath(pkg,".git") || return typemin(VersionNumber)
    head = Git.head(dir=pkg)
    vers = [keys(filter((ver,info)->info.sha1==head, avail))...]
    !isempty(vers) && return maximum(vers)
    cache = Cache.path(pkg)
    cache_has_head = isdir(cache) && Git.iscommit(head, dir=cache)
    ancestors = VersionNumber[]
    descendants = VersionNumber[]
    for (ver,info) in avail
        sha1 = info.sha1
        base = if cache_has_head && Git.iscommit(sha1, dir=cache)
            Git.readchomp(`merge-base $head $sha1`, dir=cache)
        elseif Git.iscommit(sha1, dir=pkg)
            Git.readchomp(`merge-base $head $sha1`, dir=pkg)
        else
            Base.warn_once("unknown $pkg commit $(sha1[1:8]), metadata may be ahead of package cache")
            continue
        end
        base == sha1 && push!(ancestors,ver)
        base == head && push!(descendants,ver)
    end
    both = sort!(intersect(ancestors,descendants))
    isempty(both) || warn("$pkg: some versions are both ancestors and descendants of head: $both")
    if !isempty(descendants)
        v = minimum(descendants)
        return VersionNumber(v.major, v.minor, v.patch, ("",), ())
    elseif !isempty(ancestors)
        v = maximum(ancestors)
        return VersionNumber(v.major, v.minor, v.patch, (), ("",))
    else
        return typemin(VersionNumber)
    end
end

function requires_path(pkg::String, avail::Dict=available(pkg))
    pkgreq = joinpath(pkg,"REQUIRE")
    ispath(pkg,".git") || return pkgreq
    Git.dirty("REQUIRE", dir=pkg) && return pkgreq
    !Git.success(`cat-file -e HEAD:REQUIRE`, dir=pkg) && isfile(pkgreq) && return pkgreq
    head = Git.head(dir=pkg)
    for (ver,info) in avail
        if head == info.sha1
            return joinpath("METADATA", pkg, "versions", string(ver), "requires")
        end
    end
    return pkgreq
end

function requires_list(pkg::String, avail::Dict=available(pkg))
    reqs = filter!(Reqs.read(requires_path(pkg,avail))) do line
        isa(line,Reqs.Requirement)
    end
    map(req->req.package, reqs)
end
requires_dict(pkg::String, avail::Dict=available(pkg)) =
    Reqs.parse(requires_path(pkg,avail))

function installed(avail::Dict=available())
    pkgs = Dict{ByteString,(VersionNumber,Bool)}()
    for pkg in readdir()
        isinstalled(pkg) || continue
        ap = get(avail,pkg,Dict{VersionNumber,Available}())
        pkgs[pkg] = (installed_version(pkg,ap),isfixed(pkg,ap))
    end
    return pkgs
end

function fixed(avail::Dict=available(), inst::Dict=installed(avail),
    julia_version::VersionNumber=VERSION)
    pkgs = Dict{ByteString,Fixed}()
    for (pkg,(ver,fix)) in inst
        fix || continue
        ap = get(avail,pkg,Dict{VersionNumber,Available}())
        pkgs[pkg] = Fixed(ver,requires_dict(pkg,ap))
    end
    pkgs["julia"] = Fixed(julia_version)
    return pkgs
end

function free(inst::Dict=installed())
    pkgs = Dict{ByteString,VersionNumber}()
    for (pkg,(ver,fix)) in inst
        fix && continue
        pkgs[pkg] = ver
    end
    return pkgs
end

end # module

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia