text:
using DataFrames
using DataStructures
using Logging

Base.length(it::DataStructures.SDMIncludeLast) = it.last - it.first + 1

"""
Merge columns from a dataframe to another one which is considered main.
The firts data column in each data frame is for arguments.
So do no use in mapping object the first column numver (it is argument column).
"""
function merge_data_frame_to_main_by_columns(main::DataFrame, extra::DataFrame;
    mapping_columns::Nullable{Dict{Int, Int}}=Nullable{Dict{Int, Int}}(),
    extra_boundaries::Nullable{Tuple{Float64, Float64}}=Nullable{Tuple{Float64, Float64}}(),
    main_extrapolation::Symbol=:zero, extra_extrapolation::Symbol=:zero, extra_transform::Function=y->y,
    smoothing_minpoints::Tuple{Int, Int}=(5, 5), smoothing_sharpness::Float64=-1.0)
  Logging.configure(level=INFO)
  info("Merging the selected extra table data to main ones. The mapping rules are defined as:\n$mapping_columns")
  info("Main table has size $(size(main)); Extra table has size $(size(extra))")
  info("Checking the input argument...")
  mapping = Dict{Int, Int}()
  if (isnull(mapping_columns) || isempty(get(mapping_columns))) && size(extra, 2) == size(main, 2)
    for i = 1:size(main, 2)-1
      mapping[i] = i
    end
    info("Empty mapping is defined and the main column quantity equals to the extra one. So we create one to one mapping rule by defaut.")
  elseif isnull(mapping_columns) && size(extra, 2) ≠ size(main, 2)
    error("No mapping is not specified for unequal column quantities in 'main'($(size(main, 2)-1)) and 'columns'($(size(extra, 2)-1))")
  elseif !isnull(mapping_columns)
    mapping = get(mapping_columns)
    @assert all(colindex -> 1 <= colindex <= size(extra, 2)-1, keys(mapping))
      prerequisite("A columns index in mapping must be not greater than the column dataframe column size", ("indices", keys(mapping)));
    @assert all(mainindex -> 1 <= mainindex <= size(main, 2)-1, values(mapping))
      prerequisite("A main index in mapping must be not greater than the main dataframe column size", ("indices", values(mapping)));
  end
  info("Done...")
  # ----

  main_mapping = map(reverse, mapping)
  info("The maping rules are reversed: $main_mapping")
  extra_b = isnull(extra_boundaries) ? (extra[1,1], extra[end,1]) : get(extra_boundaries)
  info("Going to make merging for extra data in boundaries [$(extra_b[1]), $(extra_b[2])]")

  info("Making spline functions out of main and extra data tables...")
  FC = convert_functionDataFrame_to_Functions(extra; extrapolation=extra_extrapolation)
  Fᵉ = map(fcf->begin x->extra_transform(fcf(x)) end, FC)
  Fᵐ = convert_functionDataFrame_to_Functions(main; extrapolation=main_extrapolation)
  info("Testing the data functions:")
  info("Testing Fᵐ:")
  for i=1:length(Fᵐ)
    fᵐ = Fᵐ[i]
    xarg = (main[1,1] + main[end,1]) / 2
    info("Fᵐ[$i]($xarg) = $(fᵐ(xarg))")
  end
  info("Testing Fᵉ:")
  for i=1:length(Fᵉ)
    fᵉ = Fᵉ[i]
    xarg = (extra_b[1] + extra_b[2]) / 2
    info("Fᵉ[$i]($xarg) = $(fᵉ(xarg))")
  end
  info("Done...")
  info("Making a function array mixture according to the inverted mapping rules $(main_mapping)...")
  Lᵐ = length(Fᵐ)
  Fᵐᵉ = copy(Fᵐ)
  for i = 1:length(Fᵐ)
    if haskey(main_mapping, i)
      iᶜ = main_mapping[i]
      Fᵐᵉ[i] = Fᵉ[iᶜ]
      info("Mapping rule $i → $iᶜ; the extra spline function is substituted.")
    end
  end
  info("Testing Fᵐᵉ:")
  for i=1:length(Fᵐᵉ)
    fᵐᵉ = Fᵐᵉ[i]
    xarg = (extra_b[1] + extra_b[end]) / 2
    info("Fᵐᵉ[$i]($xarg) = $(fᵐᵉ(xarg))")
  end
  info("Done...")

  function_dict = Dict(:main=>Fᵐ, :extra=>Fᵉ)
  info("The most complicated part... Making a merged argument grid in order to reach maximum precision and smoothing...")
  intervals = calculate_intervals_with_smoothing(
    [
      :main=>collect(main[:,1]),
      :extra=>(isnull(extra_boundaries) ? collect(extra[:,1]) : collect(extra[extra_b[1] .<= extra[:,1] .<= extra_b[2], 1]))
    ];
    smoothing_minpoints=smoothing_minpoints
  )
  info("The following intervals are evaluated for a good new grid datatable:")
  for interval in intervals
    info("Interval: type - $(interval.first); table name - $(interval.second[1]); range - [$(interval.second[2][1]), $(interval.second[2][end])]")
  end

  info("Making new merged data table...")
  newmain = createDataFrameWithColumns(main.colindex.names; tp=Float64)
  for interval in intervals
    info("Working with interval: type - $(interval.first); table name - $(interval.second[1]); range - [$(interval.second[2][1]), $(interval.second[2][end])] ...")
    interval_type = interval.first
    related_intervals = interval.second[1]
    X = interval.second[2]

    @assert interval_type == :plain || interval_type == :smoothing
    if interval_type == :plain
      @assert length(related_intervals) >= 1
    elseif interval_type == :smoothing
      if length(related_intervals) ≠ 2
        warn("length(related_intervals) ≠ 2, related_intervals=$related_intervals")
      end
    end

    info("Selecting functions array...")
    funcs = begin
      if interval_type == :plain
        if length(related_intervals) == 1 && related_intervals[1] == :main
          info("Main data table functions are selected")
          Fᵐ
        elseif (length(related_intervals) == 1 && related_intervals[1] == :extra) ||
          (length(related_intervals) == 2 && :extra ∈ related_intervals)
          info("Merged data table functions are selected")
          Fᵐᵉ
        end
      elseif interval_type == :smoothing
        Fˢ = copy(Fᵐ)
        info("Smoothed data table functions are going to be evaluated...")
        for i = 1:length(Fˢ)
          if haskey(main_mapping, i) && length(related_intervals) == 2
            iᶜ = main_mapping[i]
            info("Making a smoothing function for the mapping rule $i → $iᶜ")
            xˡ = X[1]; xʳ = X[end]
            info("The smoothing interval is [$xˡ, $xʳ]")
            m_e = related_intervals[1] == :main && related_intervals[2] == :extra
            info("The sequence of smoothing intervals: $(related_intervals[1]), $(related_intervals[2]); $(m_e ? "main → extra" : "extra → main")")
            fˡ, fʳ = m_e ? (Fᵐ[i](xʳ), Fᵉ[iᶜ](xˡ)) : (Fᵉ[iᶜ](xʳ), Fᵐ[i](xˡ))

            if fˡ === NaN || fʳ === NaN
              warn("fˡ=$fˡ, fʳ=$fʳ, m_e=$m_e; xˡ=$xˡ, xʳ=$xʳ")
            end

            σf = sign(fʳ - fˡ)
            αˢʰᵃʳᵖ = if smoothing_sharpness < 0 10 else smoothing_sharpness end
            x₀ = (X[1] + X[end])/2; α = σf * abs(X[end]-X[1]) / αˢʰᵃʳᵖ
            info("Sigmoid parameters: σf=$σf, x₀=$x₀, α=$α; (σf=$σf)")
            smoothing = if σf > 0
                if m_e
                  info("σf > 0, Fᵐ[$i], Fᵉ[$iᶜ]")
                  sigmoid_of_name(Fᵐ[i], Fᵉ[iᶜ], x₀, α, "Fᵐ[$i], Fᵉ[$iᶜ]")
                else
                  info("σf > 0, Fᵉ[$iᶜ], Fᵐ[$i]")
                  sigmoid_of_name(Fᵉ[iᶜ], Fᵐ[i], x₀, α, "Fᵉ[$iᶜ], Fᵐ[$i]")
                end
              else
                if m_e
                  info("σf < 0, Fᵉ[$iᶜ], Fᵐ[$i]")
                  sigmoid_of_name(Fᵉ[iᶜ], Fᵐ[i], x₀, α, "Fᵉ[$iᶜ], Fᵐ[$i]")
                else
                  info("σf < 0, Fᵐ[$iᶜ], Fᵉ[$i]")
                  sigmoid_of_name(Fᵐ[i], Fᵉ[iᶜ], x₀, α, "Fᵐ[$i], Fᵉ[$iᶜ]")
                end
              end
            info("Smooting function for the rule $i → $iᶜ is done; σ($xˡ)=$(smoothing(xˡ)) ↔∿↔ σ($xʳ)=$(smoothing(xʳ))")
            Fˢ[i] = smoothing
          else
            if length(related_intervals) ≠ 2
              warn("No smoothing for the related intervals $related_intervals for function $i; using main function data...")
            end
          end
        end
        info("Smoothing functions done: $Fˢ")
        Fˢ
      else
        error("The following interval types are supported: $(:plain), $(:smoothing)")
      end
    end

    info("Filling the new data table in the interval [$(X[1]), $(X[end])]")
    info("Funcs array: $funcs")
    for x in X
      row = zeros(Float64, Lᵐ+1)
      row[1] = x
      for i = 1:Lᵐ
        row[i+1] = funcs[i](x)
      end
      push!(newmain, row)
    end
    info("Done...")

  end
  info("Done!")
  return newmain
end

# -----------------------------------------------------------------------------

type Edge
  value::Float64
  intervals::Dict{Symbol, Symbol}
  middles::Set{Symbol}
end

function calculate_intervals_with_smoothing(VX::Vector{Pair{Symbol, Vector{Float64}}}; smoothing_minpoints::Tuple{Int, Int}=(5, 5))
  Logging.configure(level=INFO)
  info("Calculation of intervals with smoothing...")
  info("The given interval names: $(map(it->it.first, VX))")
  info("Interval properties:")
  for ititem in VX
    info("Interval: name - $(ititem.first); range - [$(ititem.second[1]), $(ititem.second[end])]")
  end
  # ------
  # =====
  info("Checking the input data...")
  @assert !isempty(VX) && length(VX) > 1
  VX_denotions = map(vx -> vx.first, VX)
  @assert length(VX_denotions) == length(unique(VX_denotions)) prerequisite("Vectors symbols must be unique", ("symbols", VX_denotions))
  for (i, VXⁱ) in enumerate(VX)
    @assert !isempty(VXⁱ) && length(VXⁱ.second) > 1 prerequisite(
      "The $i-th interval $(VXⁱ.first) [$(VXⁱ.second[1]), $(VXⁱ.second[end]))] must have more than one element", ("length", length(VXⁱ.second)))
  end
  σVX = map(vx -> sign(vx.second[end] - vx.second[1]), VX)
  @assert abs(reduce((sign1, sign2) -> sign1 + sign2, σVX)) == length(VX) prerequisite(
    "All VX intervals must have the same non-zero sign.", ("σVX", σVX))
  for (i, VXⁱ) in enumerate(VX)
    @assert issorted(VXⁱ.second; rev = σVX[i] == -1) prerequisite(
      "The $i-th interval $(VXⁱ.first) [$(VXⁱ.second[1]), $(VXⁱ.second[end]))] must be sorted.", ("σVX[$i]", σVX[i]))
  end
  info("Done...")
  # =====
  σ = 1
  MX = Dict(VX)
  grid = unique(sort(collect(Base.flatten(map(vx -> vx.second, VX))); rev = σ == -1))
  info("Merged grid has length $(length(grid))")

  intervals = Dict(map(it->it.first=>(it.second[1], it.second[end]), VX))
  info("Intervals dictionary: $intervals")
  edges_raw = sort(collect(
      Base.flatten(
        map(
          it -> [(it.first, :begin, it.second[1]), (it.first, :end, it.second[end])],
          VX
        )
      )
    );
    lt = (e1, e2) -> e1[3] < e2[3], rev = σ == -1
  )
  edge_values = unique(sort(map(edge->edge[3], edges_raw); rev = σ == -1))
  info("Edge points: $edges_raw")
  info("Edge values: $edge_values")

  edges = Vector{Edge}()
  for edge_value in edge_values
    info("Making Edge object for the edge value $edge_value")
    iedges_by_intervals = collect(find(e->abs(e[3] - edge_value) < 1e-8, edges_raw))
    edges_by_intervals = getindex(edges_raw, iedges_by_intervals)
    info("Edges by intervals: indices - $iedges_by_intervals, edges - $edges_by_intervals")
    intervals_for_edge = Dict(map(e->e[1]=>e[2], edges_by_intervals))
    middles = Set{Symbol}(keys(filter((ik, iv)->iv[1] < edge_value < iv[2], intervals)))
    edge = Edge(edge_value, intervals_for_edge, middles)
    push!(edges, edge)
  end
  info("Edge objects: $edges")

  intervals_out = Vector{Pair{Symbol, Tuple{Vector{Symbol}, Vector{Float64}}}}()
  edge_previous = edges[1]; grid_previous = 1
  for ie = 1:length(edges)
    edge = edges[ie]
    info("EDGE: $edge")
    first_edge = all(edge_type->edge_type==:begin, values(edge.intervals)) && isempty(edge.middles)
    last_edge = all(edge_type->edge_type==:end, values(edge.intervals)) && isempty(edge.middles)
    if first_edge
      edge_previous = edge
      grid_previous = 1
    elseif last_edge
      Xˡ = grid[grid_previous:end]
      Xˡ_names = collect(
        keys(
          filter(
            (k, v) -> begin
              matches = v[1] <= Xˡ[1] < v[2] && v[1] < Xˡ[end] <= v[2]
              if matches
                info("The interval $k → $v matches the left grid interval [$(Xˡ[1]), $(Xˡ[end])].")
                true
              else
                ΔX = abs(Xˡ[end] - Xˡ[1]); Δiv = abs(v[2] - v[1])
                ϵ = (1 - min(ΔX, Δiv)/max(ΔX, Δiv))*100
                if ϵ < 5 && Xˡ[end] <= v[2]
                  warn("The interval $k → $v matches the left grid interval [$(Xˡ[1]), $(Xˡ[end])] partially; ϵ=$ϵ, Xˡ[end]=$(Xˡ[end]), v[2]=$(v[2]).")
                  true
                else
                  info("The interval $k → $v _does_not_ match the left grid interval [$(Xˡ[1]), $(Xˡ[end])]; ϵ=$ϵ, Xˡ[end]=$(Xˡ[end]), v[2]=$(v[2]).")
                  false
                end
              end
            end,
            intervals
          )
        )
      )

      info("LAST INTERVAL: [$(Xˡ[1]), $(Xˡ[end])], names: $(Xˡ_names)")

      push!(intervals_out, :plain=>(Xˡ_names, Xˡ))
    else
      info("Handling a middle edge: $(edge)")
      left_names = map(it->it.first, collect(filter((ik, iv)->iv==:end, edge.intervals)))
      right_names = map(it->it.first, collect(filter((ik, iv)->iv==:begin, edge.intervals)))
      info("On the left: $left_names")
      info("On the right: $right_names")
      middle_intervals = isempty(edge.middles) ? Dict{Symbol, Vector{Float64}}() : filter((ik, iv)->ik ∈ edge.middles, MX)
      info("Middles: $(keys(middle_intervals))")

      Δxˡ = isempty(left_names) ? 0 : maximum(
        map(
          it->begin
            Δxⁱᵗ = abs(it.second[end]-it.second[end-1])
            info("Collecting the left step interval for the interval [$(it.second[1]), $(it.second[end])] with $(length(it.second)) points, Δxⁱᵗ=$Δxⁱᵗ")
            Δxⁱᵗ
          end
        , collect(filter((ik, iv)->ik ∈ left_names, MX))))
      Δxʳ = isempty(right_names) ? 0 : maximum(
        map(
          it->begin
            Δxⁱᵗ = abs(it.second[1]-it.second[2])
            info("Collecting the right step interval for the interval [$(it.second[1]), $(it.second[end])] with $(length(it.second)) points, Δxⁱᵗ=$Δxⁱᵗ")
            Δxⁱᵗ
          end
        , collect(filter((ik, iv)->ik ∈ right_names, MX))))
      info("Δxˡ=$Δxˡ, Δxʳ=$Δxʳ")

      Δxᵐˡ = isempty(middle_intervals) ? 0 : maximum(map(it->begin
                            ixl = findlast(x->x < edge.value, it)
                            return abs(it[ixl+1] - it[ixl])
                          end,
                          values(middle_intervals)))
      Δxᵐʳ = isempty(middle_intervals) ? 0 : maximum(map(it->begin
                            ixl = findfirst(x->x > edge.value, it)
                            return abs(it[ixl] - it[ixl-1])
                          end,
                          values(middle_intervals)))
      info("Δxᵐˡ=$Δxᵐˡ, Δxᵐʳ=$Δxᵐʳ")
      Δxˢ = maximum([Δxˡ, Δxʳ, Δxᵐˡ, Δxᵐʳ])
      info("Δxˢ=$(Δxˢ)")

      ix = findfirst(x->abs(x-edge.value) < 1e-8, grid); @assert ix > 0
      xˡ = edge.value-Δxˢ
      info("Looking for xˡ=$xˡ (Δxˢ=$Δxˢ) in the grid [$(grid[1]), $(grid[end])]")
      ixˡ = findlast(x-> x < xˡ, grid)
      if ixˡ <= 0
        warn("ixˡ ≤ 0 → ixˡ=$ixˡ for edge $(edge.value) over the grid [$(grid[1]), $(grid[end])]; Assuming ixˡ=1...")
        ixˡ = 1
      end
      xʳ = edge.value+Δxˢ
      info("Looking for xʳ=$xʳ (Δxˢ=$Δxˢ) in the grid [$(grid[1]), $(grid[end])]")
      ixʳ = findfirst(x->x > xʳ, grid); @assert ixʳ > 0
      if ixʳ <= 0
        warn("ixʳ ≤ 0 → ixʳ=$ixʳ for edge $(edge.value) over the grid [$(grid[1]), $(grid[end])]; Assuming ixʳ=1...")
        ixʳ = 1
      end
      info("ix=$ix, ixˡ=$ixˡ, ixʳ=$ixʳ")

      pᵐⁱⁿ⁻ˡ = smoothing_minpoints[1]; pᵐⁱⁿ⁻ʳ = smoothing_minpoints[2]
      pˡ = abs(ix - ixˡ); pˡ = pˡ < pᵐⁱⁿ⁻ˡ ? pᵐⁱⁿ⁻ˡ : pˡ
      pʳ = abs(ixʳ - ix); pʳ = pʳ < pᵐⁱⁿ⁻ʳ ? pᵐⁱⁿ⁻ʳ : pʳ
      info("pˡ=$pˡ, pʳ=$pʳ")

      ixˡ = if ixˡ == 1 && pˡ > 0 ixˡ else ix - pˡ end
      ixʳ = if ixʳ == 1 && pʳ > 0 ixʳ elseif ixˡ == 1 && pˡ > 0 1 else ix + pʳ end
      info("ix=$ix, ixˡ=$ixˡ, ixʳ=$ixʳ; xˡ=$(grid[ixˡ]), xʳ=$(grid[ixʳ]), x=$(grid[ix])")

      Xˡ = if ixˡ == 1 && pˡ > 0 (Float64)[] else grid[grid_previous:ixˡ-1] end
      if isempty(Xˡ)
        warn("Left interval for the point x=$(grid[ix]) is empty - skipping it and a smoothing too...")
      else
        info("Looking for names for the interval [$(Xˡ[1]), $(Xˡ[end])] in the intervals dict: $intervals")
        Xˡ_names = collect(
          keys(
            filter(
              (k, v) -> begin
                matches = v[1] <= Xˡ[1] < v[2] && v[1] < Xˡ[end] <= v[2]
                if matches
                  info("The interval $k → $v matches the left grid interval [$(Xˡ[1]), $(Xˡ[end])].")
                  true
                else
                  ΔX = abs(Xˡ[end] - Xˡ[1]); Δiv = abs(v[2] - v[1])
                  ϵ = (1 - min(ΔX, Δiv)/max(ΔX, Δiv))*100
                  if ϵ < 5 && Xˡ[end] <= v[2]
                    warn("The interval $k → $v matches the left grid interval [$(Xˡ[1]), $(Xˡ[end])] partially; ϵ=$ϵ, Xˡ[end]=$(Xˡ[end]), v[2]=$(v[2]).")
                    true
                  else
                    info("The interval $k → $v _does_not_ match the left grid interval [$(Xˡ[1]), $(Xˡ[end])]; ϵ=$ϵ, Xˡ[end]=$(Xˡ[end]), v[2]=$(v[2]).")
                    false
                  end
                end
              end,
              intervals
            )
          )
        )
        info("LEFT INTERVAL: [$(Xˡ[1]), $(Xˡ[end])], names: $(Xˡ_names)")
        push!(intervals_out, :plain=>(unique(Xˡ_names), Xˡ))

        Xˢ = collect(linspace(grid[ixˡ], grid[ixʳ], 3*abs(ixʳ-ixˡ+1)))
        info("Smooting interval for the point $(grid[ix]) has $(length(Xˢ)) point(s).")
        Xˢ_names = Vector{Symbol}()
        append!(Xˢ_names, Xˡ_names)
        Xʳ_min = reduce((it1, it2)->begin
                  i1 = it1.second
                  i2 = it2.second
                  return abs(i1[end] - i1[1]) < abs(i2[end]-i2[1]) ? it1 : it2
                end, collect(filter((k, v)->v[1] < Xˢ[end] < v[2], intervals)))
        info("Right smoothing minimal interval: $(Xʳ_min)")
        Xʳ_name = Xʳ_min.first
        info("Smoothing names on the right: $(Xʳ_name)")
        append!(Xˢ_names, [Xʳ_name])
        info("SMOOTHING INTERVAL: [$(Xˢ[1]), $(Xˢ[end])], names: $(Xʳ_name)")
        push!(intervals_out, :smoothing=>(unique(Xˢ_names), Xˢ))
      end

      edge_previous = edge
      grid_previous = if ixˡ == 1 && pˡ > 0 ixʳ else ixʳ + 1 end
    end
  end
  info("Done ⌣")
  return intervals_out
end

# -----------------------------------------------------------------------------

function get_affected_interval_pair(interval::Tuple{Float64, Float64}, intervals::Dict{Symbol, Tuple{Float64, Float64}})
  Logging.configure(level=INFO)
  info("Evaluating closest interval pair for the interval values $interval")
  left = interval[1]; right = interval[2]
  left_diffs = map(it->(it.first, abs(it.second[2] - left)), intervals)
  right_diffs = map(it->(it.first, abs(right - it.second[1])), intervals)

  min_left = Nullable{Pair{Symbol, Tuple{Float64, Float64}}}()
  min_right = Nullable{Pair{Symbol, Tuple{Float64, Float64}}}()

  left_result = reduce((d1, d2) -> d1[2] < d2[2] ? d1 : d2, left_diffs)
  right_result = reduce((d1, d2) -> d1[2] < d2[2] ? d1 : d2, right_diffs)

  min_left = Nullable{Pair{Symbol, Tuple{Float64, Float64}}}(left_result.first=>intervals[left_result.first])
  min_right = Nullable{Pair{Symbol, Tuple{Float64, Float64}}}(right_result.first=>intervals[right_result.first])

  info("Done...")
  return min_left, min_right
end

function get_affected_interval_minimal(interval::Tuple{Float64, Float64}, intervals::Dict{Symbol, Tuple{Float64, Float64}})
  Logging.configure(level=INFO)
  info("Evaluating minimal interval for the interval values $interval")
  @assert interval[1] ≠ interval[2]
  entries = map(it->it[2], enumerate(intervals))
  min_interval = Nullable{Pair{Symbol, Tuple{Float64, Float64}}}()
  for entry in entries[2:end]
    name = entry.first
    b = entry.second
    if (b[1] <= interval[1] <= b[2] && b[1] <= interval[2] <= b[2]) || (b[2] <= interval[1] <= b[1] && b[2] <= interval[2] <= b[1])
      if isnull(min_interval)
        min_interval = Nullable{Pair{Symbol, Tuple{Float64, Float64}}}(entry)
      else
        mb = get(min_interval).second
        if abs(b[2] - b[1]) < abs(mb[2] - mb[1])
          min_interval = Nullable{Pair{Symbol, Tuple{Float64, Float64}}}(entry)
        end
      end
    end
  end
  info("Done...")
  return min_interval
end

function get_affected_intervals(interval::Tuple{Float64, Float64}, intervals::Dict{Symbol, Tuple{Float64, Float64}})
  @assert interval[1] ≠ interval[2]
  interval_names = Vector{Symbol}()
  entries = map(it->it[2], enumerate(intervals))
  for entry in entries
    name = entry.first
    b = entry.second
    if (b[1] <= interval[1] <= b[2] && b[1] <= interval[2] <= b[2]) || (b[2] <= interval[1] <= b[1] && b[2] <= interval[2] <= b[1])
      push!(interval_names, name)
    end
  end
  return interval_names
end

"""
Sort a given set of data frames and create a new data frame out of them.
"""
function merge_data_frames_sorted_sequencially(
  data::Vector{DataFrame}; points_to_smooth::Int = 5,
  piece_merging::Function=(df1::DataFrame, df2::DataFrame, x1::Float64, x2::Float64) ->
  begin createDataFrameWithColumns(df1.colindex.names; tp=typeof(df1[1,1])) end)
  Logging.configure(level=INFO)
  if isempty(data)
    info("The vector of data items is empty - nothing to do.")
    return DataFrame()
  elseif length(data) == 1
    info("The vector of data items has a single data frame simply - nothing to do.")
    return data[1]
  end
  data_map = Dict{Tuple{Float64, Float64}, DataFrame}()
  info("Merging the data pieces of size '$(length(data))'.")
  for data_piece in data
    # the first column is considered as an argument one
    Rᵃ = data_piece[1,1]
    Rᵇ = data_piece[end,1]
    @assert Rᵃ < Rᵇ prerequisite("Rᵃ < Rᵇ", ("Rᵃ", Rᵃ), ("Rᵇ", Rᵇ))
    info("A data piece is defined in the argument interval [$Rᵃ, $Rᵇ].")
    data_map[(Rᵃ, Rᵇ)] = data_piece
  end
  info("Data with the folowing intervals are going to be ordered and merged: $(sort(collect(keys(data_map))))")
  intervals = sort(collect(keys(data_map)))
  # ---
  for i = 1:length(intervals)-1
    interval_1 = intervals[i]; interval_2 = intervals[i+1]
    Rᵃ¹ = interval_1[1]; Rᵇ¹ = interval_1[2]
    Rᵃ² = interval_2[1]; Rᵇ² = interval_2[2]
    info("Checking data in intervals [$Rᵃ¹, $Rᵇ¹] and [$Rᵃ², $Rᵇ²]...")
    @assert Rᵃ¹ < Rᵇ² prerequisite("Rᵃ¹ < Rᵇ² for $i, $(i+1)", ("Rᵃ¹", Rᵃ¹), ("Rᵇ²", Rᵇ²))
    @assert Rᵇ¹ <= Rᵃ² prerequisite("Rᵇ¹ ≤ Rᵃ² for $i, $(i+1)", ("Rᵇ¹", Rᵇ¹), ("Rᵃ²", Rᵃ²))
    cols_1 = size(data_map[interval_1], 2); cols_2 = size(data_map[interval_2], 2)
    @assert cols_1 == cols_1 prerequisite("cols¹ == cols² for $i, $(i+1); Data are not mergeable.", ("cols¹", cols_1), ("cols²", cols_2))
    col_names_1 = data_map[interval_1].colindex.names; col_names_2 = data_map[interval_2].colindex.names
    @assert col_names_1 == col_names_2 prerequisite("col_names_1 == col_names_2 for $i, $(i+1)", ("col_names_1", col_names_1), ("col_names_2", col_names_2))
  end
  info("Seems that the data are mergeable. Let's try it...")
  # ---
  # ---
  data_type = typeof(data_map[intervals[1]][1,1]) # TODO: check all column data types before merge
  data_frame = createDataFrameWithColumns(data_map[intervals[1]].colindex.names; tp=data_type)
  info("Going to do a merge of $(length(data)) pieces...")
  for i = 1:length(intervals)-1
    interval_1 = intervals[i]; interval_2 = intervals[i+1]
    Rᵃ¹ = interval_1[1]; Rᵇ¹ = interval_1[2]
    Rᵃ² = interval_2[1]; Rᵇ² = interval_2[2]
    info("Merging data intervals [$Rᵃ¹, $Rᵇ¹] and [$Rᵃ², $Rᵇ²]...")
    data_1 = data_map[interval_1][1:end-points_to_smooth-1,:]
    info("First data block: $(size(data_1))")
    data_2 = data_map[interval_2][points_to_smooth+1:end,:]
    info("Second data block: $(size(data_2))")
    x1 = data_map[interval_1][end-points_to_smooth,1]
    x2 = data_map[interval_2][points_to_smooth,1]
    info("A smoothing area is [$x1, $x2].")
    data_12 = piece_merging(data_map[interval_1], data_map[interval_2], x1, x2)
    append!(data_frame, data_1)
    append!(data_frame, data_12)
    append!(data_frame, data_2)
    info("Done merging [$Rᵃ¹, $Rᵇ¹] and [$Rᵃ², $Rᵇ²].")
  end
  info("Done merging of $(length(data)) pieces.")
  return data_frame
end

"""
We assume that the first column in data frame object contains data points.
"""
function merge_function_table_piece_to_main_one(main::DataFrame, piece::DataFrame;
    argframe::Tuple{Float64, Float64}=(-1.0, -1.0),
    piece_ycolumns::Vector{Int}=Vector{Int}(1:size(piece, 2)))
  Nₘ = size(main, 1); L = size(main, 2) - 1
  Nₚ = size(piece, 1); Lₚ = size(piece, 2) - 1
  Xᵖ = convert(Vector{Float64}, piece[:,1])
  Xᵐ = convert(Vector{Float64}, main[:,1])
  xₐ₁ = argframe[1]; xₐ₂ = argframe[2]
  x₁ₚ = first(Xᵖ); x₂ₚ = last(Xᵖ)
  x₁ₘ = first(Xᵐ); x₂ₘ = last(Xᵐ)
  @assert Lₚ == L prerequisite("Lₚ == L", ("Lₚ", Lₚ), ("L", L))
  @assert x₁ₚ <= xₐ₁ && xₐ₂ <= x₂ₚ prerequisite("x₁ₚ <= xₐ₁ && xₐ₂ <= x₂ₚ", ("x₁ₚ", x₁ₚ), ("xₐ₁", xₐ₁), ("x₂ₚ", x₂ₚ), ("xₐ₂", xₐ₂))
  @assert x₁ₘ <= x₁ₚ && x₂ₚ <= x₂ₘ prerequisite("x₁ₘ <= x₁ₚ && x₂ₚ <= x₂ₘ", ("x₁ₘ", x₁ₘ), ("x₁ₚ", x₁ₚ), ("x₂ₚ", x₂ₚ), ("x₂ₘ", x₂ₘ))

  @assert length(piece_ycolumns) <= L prerequisite("length(piece_ycolumns) <= L", ("length(piece_ycolumns)", length(piece_ycolumns)), ("L", L))
  for col in piece_ycolumns
    @assert 1 <= col <= L prerequisite("1 ≤ col ≤ L in piece_ycolumns", ("col", col), ("L", L))
  end
  pycols = isempty(piece_ycolumns) ? Vector{Int}(1:1:L) : Vector{Int}(piece_ycolumns)

  ixₚˡᵉᶠᵗ = findlast(x -> x < xₐ₁, Xᵖ) + 1; ixₚʳⁱᵍʰᵗ = findfirst(x -> x > xₐ₂, Xᵖ) - 1
  xₚˡᵉᶠᵗ = Xᵖ[ixₚˡᵉᶠᵗ]; xₚʳⁱᵍʰᵗ = Xᵖ[ixₚʳⁱᵍʰᵗ]
  ixₘˡᵉᶠᵗ = findlast(x -> x < xₚˡᵉᶠᵗ, Xᵐ); ixₘʳⁱᵍʰᵗ = findfirst(x -> x > xₚʳⁱᵍʰᵗ, Xᵐ)
  xₘˡᵉᶠᵗ = Xᵐ[ixₘˡᵉᶠᵗ]; xₘʳⁱᵍʰᵗ = Xᵐ[ixₘʳⁱᵍʰᵗ]
  @assert xₘˡᵉᶠᵗ < xₚˡᵉᶠᵗ prerequisite("xₘˡᵉᶠᵗ < xₚˡᵉᶠᵗ", ("xₘˡᵉᶠᵗ", xₘˡᵉᶠᵗ), ("xₚˡᵉᶠᵗ", xₚˡᵉᶠᵗ))
  @assert xₚʳⁱᵍʰᵗ < xₘʳⁱᵍʰᵗ prerequisite("xₚʳⁱᵍʰᵗ < xₘʳⁱᵍʰᵗ", ("xₚʳⁱᵍʰᵗ", xₚʳⁱᵍʰᵗ), ("xₘʳⁱᵍʰᵗ", xₘʳⁱᵍʰᵗ))

  Fᵐ = convert_functionDataFrame_to_Functions(main)
  Fᵖ = convert_functionDataFrame_to_Functions(piece; skip=Vector{Int}(filter(ix -> ix ∉ pycols, 1:1:L)))

  newdata = createDataFrameWithColumns(main.colindex.names)
  for k = 1:ixₘˡᵉᶠᵗ
    row = data_row_to_vector(main, k)
    push!(newdata, row)
  end

  # --- boundary area
  Δxₘˡᵉᶠᵗ = (xₚˡᵉᶠᵗ - xₘˡᵉᶠᵗ) / 10.0
  for x = xₘˡᵉᶠᵗ+Δxₘˡᵉᶠᵗ:Δxₘˡᵉᶠᵗ:xₚˡᵉᶠᵗ-Δxₘˡᵉᶠᵗ
    row = Vector{Float64}(L+1); row[1] = x
    for l = 1:L
      row[l+1] = l ∈ pycols ? (Fᵐ[l](x) + Fᵖ[l](x))/2.0 : Fᵐ[l](x)
    end
    push!(newdata, row)
  end
  # ---

  for kₚ = ixₚˡᵉᶠᵗ:ixₚʳⁱᵍʰᵗ
    x = Xᵖ[kₚ]
    row = Vector{Float64}(L+1); row[1] = x
    for lₚ = 1:Lₚ
      row[lₚ+1] = lₚ ∈ pycols ? Fᵖ[lₚ](x) : Fᵐ[lₚ](x)
    end
    push!(newdata, row)
  end

  # --- boundary area
  Δxₘʳⁱᵍʰᵗ = (xₘʳⁱᵍʰᵗ - xₚʳⁱᵍʰᵗ) / 10.0
  for x = xₚʳⁱᵍʰᵗ+Δxₘʳⁱᵍʰᵗ:Δxₘʳⁱᵍʰᵗ:xₘʳⁱᵍʰᵗ-Δxₘʳⁱᵍʰᵗ
    row = Vector{Float64}(L+1); row[1] = x
    for l = 1:L
      row[l+1] = l ∈ pycols ? (Fᵐ[l](x) + Fᵖ[l](x))/2.0 : Fᵐ[l](x)
    end
    push!(newdata, row)
  end
  # ---

  for k = ixₘʳⁱᵍʰᵗ:Nₘ
    row = data_row_to_vector(main, k)
    push!(newdata, row)
  end

  return newdata
end


text:
@testset "Construction" begin
    @test StackSet() == StackSet()
    @test StackSet([1, 3, 2]) == StackSet([1,2,3])

    @test length(StackSet([101, 115, 151, 103, 115])) == 4
    @test length(StackSet([-50, -10, -1, -60, -5, -50])) == 5
    @test_throws ArgumentError StackSet([-10, 60])
    @test_throws ArgumentError StackSet([1001, 1070])

    # Equality between when converting to Set
    @test Set(StackSet([501, 515, 540])) == Set(StackSet([540, 515, 501]))
end

@testset "Basic" begin
    @test empty(StackSet([900, 941])) == StackSet()
    @test !isempty(StackSet([0]))
    @test isempty(StackSet([]))
    @test isempty(StackSet())
    @test length(StackSet()) == 0
    @test length(empty(StackSet())) == 0
end

@testset "Iteration" begin
    s = StackSet([577])
    v, st = iterate(s)
    @test v == 577
    @test iterate(s, st) === nothing

    @test Set([i for i in StackSet([-30, 11, 0, 1])]) == Set([1, 0, 11, -30])
end

const vectors = [
[0],
[-5, 6],
[-401, -350, -401],
collect(-4:50),
[503913, 503889, 503936],
]

@testset "Misc" begin
    @test_throws ArgumentError maximum(StackSet())
    @test_throws ArgumentError minimum(StackSet())
    for v in vectors
        @test maximum(StackSet(v)) == maximum(Set(v))
        @test minimum(StackSet(v)) == minimum(Set(v))
        @test first(iterate(v)) in StackSet(v)
        @test !(555 in StackSet(v))
    end
end

@testset "Filter" begin
    for F in [isodd, iseven, x -> true, x -> false, isequal(6)]
        for v in vectors
            @test Set(filter(F, StackSet(v))) == filter(F, Set(v))
        end
    end
end

@testset "Push" begin
    @test push(empty(StackSet()), 99) == StackSet([99])
    @test first(push(StackSet(), -501)) == -501
    for v in vectors
        s = StackSet(v)
        s´ = Set(v)
        for i in [0, -5, 11, 503889, 10040342, -9]
            if max(maximum(s), i) - min(minimum(s), i) > 63
                @test_throws ArgumentError s2 = push(s, i)
            else
                @test Set(push(s, i)) == push!(copy(s´), i)
                @test push(s, i, StackCollections.unsafe) == push(s, i)
            end
        end
    end
end

@testset "Pop & Delete" begin
    @test_throws KeyError pop(StackSet(), 0)
    @test_throws KeyError pop(StackSet([-9, 9]), 8)
    @test delete(StackSet(), 0) == StackSet()
    @test delete(StackSet([-9, 9]), 8) == StackSet([-9, 9])
    @test pop(StackSet([19, 8, 22]), 8) == StackSet([19, 22])

    for v in vectors
        s = StackSet(v)
        s´ = Set(v)
        for i in [0, -5, 11, 503889, 10040342, -9]
            if i in s
                s2´ = copy(s´)
                pop!(s2´, i)
                @test Set(pop(s, i)) == s2´
                @test Set(delete(s, i)) == delete!(copy(s´), i)
            else
                @test_throws KeyError pop(s, i)
                @test delete(s, i) == s
            end
        end
    end
end

@testset "Misc methods" begin
    vs = [[(4, 1), ()],
    [(0, 11), (11, 12)],
    [(7, 9, 51, 7), (51, 9, 9, 9, 9)],
    [(5, 3, 1, 9), (3, 5, 1, 9, 9, 38, 41)],
    [(-101, -71), (-71, -105, -89)],
    [(0, 9, -5), (1, -10, 19)],
    ]

    for F in [union, intersect, symdiff, setdiff]
        for (v1, v2) in vs
            @test F(Set{Int}(v1), Set{Int}(v2)) == Set(F(StackSet(v1), StackSet(v2)))
            @test F(Set{Int}(v2), Set{Int}(v1)) == Set(F(StackSet(v2), StackSet(v1)))

            @test F(StackSet(v1), Set(v2)) == StackSet(F(Set(v1), Set(v2)))
            @test F(StackSet(v2), Set(v1)) == StackSet(F(Set(v2), Set(v1)))
        end
    end
    for (v1, v2) in vs
        for F in [⊊, issubset]
            @test F(Set{Int}(v1), Set{Int}(v2)) == F(StackSet(v1), StackSet(v2))
            @test F(Set{Int}(v2), Set{Int}(v1)) == F(StackSet(v2), StackSet(v1))
        end
        @test allunique(StackSet(v1))
        @test allunique(StackSet(v2))
    end

    for (v1, v2) in vs
        @test issubset(Set{Int}(v1), Set{Int}(v2)) == issubset(StackSet(v1), StackSet(v2))
        @test issubset(Set{Int}(v2), Set{Int}(v1)) == issubset(StackSet(v2), StackSet(v1))
    end
end

@testset "Convert to/from DigitSet" begin
    for bad in [(-1,), (44, 10, -3), (55, 65), (101, 102, 103)]
        s = StackSet(bad)
        @test_throws ArgumentError DigitSet(s)
    end
    for good in [(0, 1, 2, 3), 14:21, 1:4:60, (61, 62, 1, 9), (), (0,)]
        s = StackSet(good)
        d = DigitSet(good)
        @test s == d
        @test collect(s) == collect(d)

        s2 = StackSet(d)
        d2 = DigitSet(s)
        @test s == d == s2 == d2
        @test collect(s) == collect(s2) == collect(d2)
    end

    @test DigitSet([1,2,39]) != StackSet([1, 2, 40])
    @test DigitSet([0,1,2]) != StackSet([5,6,7])
end

ext: jl


text:
config = MOI.Test.Config()
optimize!_max_bridged(mock) = MOIU.mock_optimize!(mock, [ 6.0, 1.0, 9.0, -3.0*√2],
    (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[3, 1, 1/3]],
    (MOI.VectorOfVariables, MOI.RotatedSecondOrderCone) => [[3, 1/3, √2]]
)
optimize!_min_bridged(mock) = MOIU.mock_optimize!(mock, [-6.0, 1.0, 9.0, 3.0*√2],
    (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[3, -1, 1/3]],
    (MOI.VectorOfVariables, MOI.RotatedSecondOrderCone) => [[3, 1/3, -√2]]
)
optimize!_max_cached(mock) = MOIU.mock_optimize!(mock, [1.0, 9.0, -3.0*√2, 6.0],
    (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[3, 1, 1/3]],
    (MOI.VectorOfVariables, MOI.RotatedSecondOrderCone) => [[3, 1/3, √2]]
)
optimize!_min_cached(mock) = MOIU.mock_optimize!(mock, [1.0, 9.0, 3.0*√2, -6.0],
    (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[3, -1, 1/3]],
    (MOI.VectorOfVariables, MOI.RotatedSecondOrderCone) => [[3, 1/3, -√2]]
)
for mock in [bridged_mock(optimize!_max_bridged, optimize!_min_bridged), cached_mock(optimize!_max_cached, optimize!_min_cached)]
    Tests.BPT12e399_rem_test(mock, config)
end
optimize!_max_bridged(mock) = MOIU.mock_optimize!(mock, [ 10.0, 4.0, 0.0, 5.0, 0.0, -5.0, 5.0],
    (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[0.0, 1, 0, 1, 1]]
)
optimize!_min_bridged(mock) = MOIU.mock_optimize!(mock, [-10.0, 4.0, 0.0, 5.0, 0.0,  5.0, 5.0],
    (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[0.0, 1, 0, -1, 1]]
)
optimize!_max_cached(mock) = MOIU.mock_optimize!(mock, [4.0, 0.0, 5.0, 0.0, -5.0, 5.0, 10.0],
    (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[0.0, 1, 0, 1, 1]]
)
optimize!_min_cached(mock) = MOIU.mock_optimize!(mock, [4.0, 0.0, 5.0, 0.0,  5.0, 5.0, -10.0],
    (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[0.0, 1, 0, -1, 1]]
)
for mock in [bridged_mock(optimize!_max_bridged, optimize!_min_bridged), cached_mock(optimize!_max_cached, optimize!_min_cached)]
    Tests.BPT12e399_maxdegree_test(mock, config)
end

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
function yalparser(target::String)

    fp = open(target, "r")
    lines = readlines(fp)

    block_x = Float64[]
    block_y = Float64[]
    block_area = Float64[]
    block_name = String[]
    arr = Vector{String}[]
    ids = Set{String}()
    W = 0.0
    H = 0.0

    it = 1
    while it <= length(lines)
        line = rstrip(chomp(lines[it]),collect(";"))
        if contains(line, "MODULE ") # start of module
            spl = split(line)
            push!(block_name, spl[2])
            it += 1
            line = rstrip(chomp(lines[it]),collect(";"))
            spl = split(line)
            typ = spl[2]
            it += 1
            line = rstrip(chomp(lines[it]),collect(";"))
            spl = split(line)
            nums = float(spl[3:end])
            x = abs(nums[1] - nums[5])
            y = abs(nums[2] - nums[6])
            push!(block_x, mean([nums[1],nums[5]]))
            push!(block_y, mean([nums[2],nums[6]]))
            push!(block_area, x*y)
            if typ == "GENERAL"
                # do nothing
            elseif typ == "PARENT"
                W = y
                H = x
                arr = Array(Vector{String}, length(block_name)-1)
                while !contains(line, "NETWORK")
                    it += 1
                    line = lines[it]
                end
                it += 1
                line = rstrip(chomp(lines[it]),collect(";"))
                while !contains(line, "ENDNETWORK")
                    spl = split(line)
                    idx = find(block_name .== spl[2])[1]
                    arr[idx] = spl[3:end]
                    for elem in spl[3:end]
                        push!(ids, elem)
                    end
                    it += 1
                    line = rstrip(chomp(lines[it]),collect(";"))
                end
            end
        end
        it += 1
    end

    # compute nets now
    nets = Vector{Int}[]
    for elem in ids
        this_net = Int[]
        for it in 1:length(arr)
            if elem in arr[it]
                push!(this_net, it)
            end
        end
        push!(nets, this_net)
    end

    return (W,H), block_area[1:end-1], nets

end

function load_yal(filesource::String, β::Real)
    (W,H), a, nets = yalparser(joinpath(Pkg.dir("FloorLayout"), "examples", "data", filesource*".yal"))
    N = length(a)

    b = β*ones(N)
    wub = min(sqrt(a.*b), W)
    hub = min(sqrt(a.*b), H)
    wlb = a ./ hub
    hlb = a ./ wub

    c = zeros(N,N)
    for net in nets
        k = length(net)
        for i in net, j in net
            if i < j
                c[i,j] += 1/(k-1)
                c[j,i] += 1/(k-1)
            end
        end
    end

    Problem(N, W, H, wlb, hlb, wub, hub, a, b, c)
end


text:
using LaserTypes
using LaserTypes: pochhammer, immutable_cache
using StaticArrays
using UnPack
using Test

x₀ = SVector{3}(0,0,0)
t₀ = 0

@testset "(:y, :x) orientation" begin
    s = setup_laser(LaguerreGaussLaser, :SI, m = 1, p = 1, profile=ConstantProfile, orientation=(:y,:x))
    Ex, Ey, Ez = E(x₀, t₀, s)
    @test iszero(Ex)
    @test iszero(Ey)
    @test iszero(Ez)

    @unpack E₀, Nₚₘ, z_R = immutable_cache(s)
    c = s.constants.c
    Bx, By, Bz = B(x₀, t₀, s)
    @test Bx ≈ - factorial(s.p)/pochhammer(abs(s.m)+1,s.p)*E₀*Nₚₘ*(√2*s.w₀)/(c*z_R)
    @test iszero(By)
    @test iszero(Bz)
end


text:
"""
    Original Code Taken from GraphLayout.jl
    Use the spring/repulsion model of Fruchterman and Reingold (1991):
        Attractive force:  f_a(d) =  d^2 / k
        Repulsive force:  f_r(d) = -k^2 / d
    where d is distance between two vertices and the optimal distance
    between vertices k is defined as C * sqrt( area / num_vertices )
    where C is a parameter we can adjust
    Arguments:
    adj_matrix Adjacency matrix of some type. Non-zero of the eltype
               of the matrix is used to determine if a link exists,
               but currently no sense of magnitude
    C          Constant to fiddle with density of resulting layout
    MAXITER    Number of iterations we apply the forces
    INITTEMP   Initial "temperature", controls movement per iteration
"""
function layout_spring{T}(adj_matrix::Array{T,2}, dim=1; C=2.0, MAXITER=100, INITTEMP=2.0)

    size(adj_matrix, 1) != size(adj_matrix, 2) && error("Adj. matrix must be square.")
    const N = size(adj_matrix, 1)

    # Initial layout is random on the square [-1,+1]^2
    locs_x = 2*rand(N) .- 1.0
    locs_y = 2*rand(N) .- 1.0
    locs_z = 2*rand(N) .- 1.0

    # The optimal distance bewteen vertices
    const K = C * sqrt(4.0 / N)

    # Store forces and apply at end of iteration all at once
    force_x = zeros(N)
    force_y = zeros(N)
    force_z = zeros(N)

    # Iterate MAXITER times
    @inbounds for iter = 1:MAXITER
        # Calculate forces
        for i = 1:N
            force_vec_x = 0.0
            force_vec_y = 0.0
            force_vec_z = 0.0
            for j = 1:N
                i == j && continue
                d_x = locs_x[j] - locs_x[i]
                d_y = locs_y[j] - locs_y[i]
                d_z = locs_z[j] - locs_z[i]
                d   = sqrt(d_x^2 + d_y^2 + d_z^2)
                if adj_matrix[i,j] != zero(eltype(adj_matrix)) || adj_matrix[j,i] != zero(eltype(adj_matrix))
                    F_d = d / K - K^2 / d^2
                else
                    F_d = -K^2 / d^2
                end

                force_vec_x += F_d*d_x
                force_vec_y += F_d*d_y
                force_vec_z += F_d*d_z
            end
            force_x[i] = force_vec_x
            force_y[i] = force_vec_y
            force_z[i] = force_vec_z
        end

        TEMP = INITTEMP / iter
        # Now apply them, but limit to temperature
        for i = 1:N
            force_mag  = sqrt(force_x[i]^2 + force_y[i]^2 + force_z[i]^2)
            scale      = min(force_mag, TEMP)/force_mag
            locs_x[i] += force_x[i] * scale
            locs_y[i] += force_y[i] * scale
            locs_z[i] += force_y[i] * scale
        end
    end

    # Scale to unit square
    min_x, max_x = minimum(locs_x), maximum(locs_x)
    min_y, max_y = minimum(locs_y), maximum(locs_y)
    min_z, max_z = minimum(locs_z), maximum(locs_z)
    function scaler(z, a, b)
        2.0*((z - a)/(b - a)) - 1.0
    end
    map!(z -> scaler(z, min_x, max_x), locs_x)
    map!(z -> scaler(z, min_y, max_y), locs_y)
    map!(z -> scaler(z, min_z, max_z), locs_z)

    if dim == 0
        locs_z = zeros(size(locs_x))
    end

    return locs_x,locs_y,locs_z
end


text:
# Build a network from a GTFS file name
using ZipFile
using CSV
using DataFrames
using Logging
using Tables
using Dates
using .OSRM

function strip_colnames!(df)
    rename!(strip, df)
end

# refactor needed - split all of these load methods for individual files out into functions
function build_network(gtfs_filenames::Vector{String}, osrm::Union{OSRMInstance, Missing}=missing, max_transfer_distance_meters::Float64=1000.0)::TransitNetwork
    # initialize a new, empty transit network
    net::TransitNetwork = TransitNetwork()

    for gtfs_filename in gtfs_filenames
        stop_number_for_stop_id::Dict{String,UInt16} = Dict()

        @info "Reading $gtfs_filename..."

        # first read stops
        r = ZipFile.Reader(gtfs_filename)
        @debug "opened file"
        filename_map = Dict([f.name=>f for f in r.files]...)
        kys = keys(filename_map)
        @debug "indexed file, keys $kys"

        @info "..stops.txt"
        stop_df = DataFrame(CSV.File(read(filename_map["stops.txt"])))
        strip_colnames!(stop_df)
        nstops = nrow(stop_df)

        # pre-allocate space for the new stops TODO can I sizehint a dict?
        # TODO should this just be nstops? are we over-allocating?
        sizehint!(net.stops, length(net.stops) + nstops)

        for srow in Tables.rows(stop_df)
            if (ismissing(srow.stop_lat) | ismissing(srow.stop_lon))
                @warn "Stop $(srow.stop_id) is missing coordinates, skipping"
            else
                stop = Stop(srow.stop_lat, srow.stop_lon)
                push!(net.stops, stop)
                net.stopidx_for_id["$gtfs_filename:$(srow.stop_id)"] = length(net.stops)
            end
        end

        @info "...Read $nstops stops"

        @info "..routes.txt"
        route_df = DataFrame(CSV.File(read(filename_map["routes.txt"]); types=Dict(:route_short_name => String, :route_long_name => String)))
        strip_colnames!(route_df)
        nroutes = nrow(route_df)

        sizehint!(net.routes, length(net.routes) + nroutes)

        for rrow in Tables.rows(route_df)
            rte = Route(rrow.route_short_name, rrow.route_long_name, rrow.route_type)
            push!(net.routes, rte)
            net.routeidx_for_id["$gtfs_filename:$(rrow.route_id)"] = length(net.routes)
        end

        @info "...Read $nroutes routes"

        if haskey(filename_map, "calendar.txt")
            @info "..calendar.txt"

            cal_df = DataFrame(CSV.File(read(filename_map["calendar.txt"])))
            strip_colnames!(cal_df)
            ncals = nrow(cal_df)

            sizehint!(net.services, length(net.services) + ncals)

            for crow in Tables.rows(cal_df)
                svc = Service(
                    crow.monday::Integer == 1,
                    crow.tuesday::Integer == 1,
                    crow.wednesday::Integer == 1,
                    crow.thursday::Integer == 1,
                    crow.friday::Integer == 1,
                    crow.saturday::Integer == 1,
                    crow.sunday::Integer == 1,
                    parse_gtfsdate(crow.start_date),
                    parse_gtfsdate(crow.end_date),
                    # calendar_dates processed later
                    Vector{Date}(),
                    Vector{Date}()
                )
                push!(net.services, svc)
                net.serviceidx_for_id["$gtfs_filename:$(crow.service_id)"] = length(net.services)
            end

            @info "...Read $ncals calendars"
        else
            @info "..calendar.txt (not present)"
        end

        if haskey(filename_map, "calendar_dates.txt")
            @info "..calendar_dates.txt"

            cal_df = DataFrame(CSV.File(read(filename_map["calendar_dates.txt"])))
            strip_colnames!(cal_df)
            ncals = nrow(cal_df)

            for crow in Tables.rows(cal_df)
                local svc::Service
                namespaced_service_id = "$gtfs_filename:$(crow.service_id)"

                if haskey(net.serviceidx_for_id, namespaced_service_id)
                    svc = net.services[net.serviceidx_for_id[namespaced_service_id]]
                else
                    svc = Service(
                        false,
                        false,
                        false,
                        false,
                        false,
                        false,
                        false,
                        Date("19700101", "yyyymmdd"),
                        Date("22500101", "yyyymmdd"),
                        Vector{Date}(),
                        Vector{Date}()
                    )
                    push!(net.services, svc)
                    net.serviceidx_for_id[namespaced_service_id] = length(net.services)
                end

                date = parse_gtfsdate(crow.date)

                if (crow.exception_type::Integer == 1)
                    push!(svc.added_dates, date)
                elseif (crow.exception_type::Integer == 2)
                    push!(svc.removed_dates, date)
                end
            end

        else
            @info "..calendar_dates.txt (not present)"
        end

        @info "..trips.txt"
        trip_df = DataFrame(CSV.File(filename_map["trips.txt"]))
        strip_colnames!(trip_df)
        ntrips = nrow(trip_df)

        sizehint!(net.trips, length(net.trips) + ntrips)

        for trow in Tables.rows(trip_df)
            service = net.serviceidx_for_id["$gtfs_filename:$(trow.service_id)"]
            route = net.routeidx_for_id["$gtfs_filename:$(trow.route_id)"]
            trp = Trip(Vector{StopTime}(), route, service, -1)
            push!(net.trips, trp)
            net.tripidx_for_id["$gtfs_filename:$(trow.trip_id)"] = length(net.trips)
        end

        @info "...loaded $ntrips trips"

        @info "..stop_times.txt"

        st_df = DataFrame(CSV.File(filename_map["stop_times.txt"], typemap=Dict(Dates.Time=>String)))
        strip_colnames!(st_df)
        nst = nrow(st_df)

        for strow in Tables.rows(st_df)
            trp = net.trips[net.tripidx_for_id["$gtfs_filename:$(strow.trip_id)"]]
            stopidx = net.stopidx_for_id["$gtfs_filename:$(strow.stop_id)"]
            st = StopTime(
                stopidx,
                strow.stop_sequence,
                parse_gtfstime(strow.arrival_time),
                parse_gtfstime(strow.departure_time)
            )
            push!(trp.stop_times, st)
        end

        @info "Loaded $nst stop times"

        if haskey(filename_map, "frequencies.txt")
            @warn "frequencies.txt present but not yet supported - frequency trips will run once at the time they are scheduled in stop times (often midnight)"
        end
    end

    @info "Loaded all GTFS files"

    @info "Sorting stop times..."
    sort_stoptimes!(net)
    @info "Done sorting stop times"

    @info "Interpolating stop times..."
    interpolate_stoptimes!(net)
    @info "Done interpolating stop times"

    @info "Finding trip patterns..."
    find_trip_patterns!(net)
    @info "Done finding trip patterns"

    if ismissing(osrm)
        @info "Finding transfers within $(max_transfer_distance_meters)m crow-flies distance..."
        find_transfers_distance!(net, max_transfer_distance_meters)
    else
        @info "Finding transfers within $(max_transfer_distance_meters)m using network distances"
        find_transfers_osrm!(net, osrm, max_transfer_distance_meters)
    end

    @info "Indexing patterns..."
    index_network!(net)

    @info "Network build completed."

    return net
end


text:
using MLStyle.Modules.Cond

usign(x) = @cond begin
    x < 0 => -1
    x == 0 => 0
    x > 0 => 1
end

@testset "cond" begin
    @test usign(-1) == -1
    @test usign(-100) == -1

    @test usign(1) == 1
    @test usign(100) == 1

    @test usign(0) == 0

    @test :ohhh == @cond begin
        false => :emmm
        _ => :ohhh
    end
end


text:
"""
    acquire_lock(lock_key[, lock_value=string(uuid4()); client::Client=get_global_client(), timeout=nothing, seconds_between_checks=0.1])

Creates a redis lock key, blocks if the lock already exists, returns the `lock_value`.

It is preferred that `redis_lock` is used instead as it will handle the acquiring and releasing of
locks within a do-block context.

# Arguments
- `lock_key`: Name of the redis lock key.
- `lock_value=string(uuid4())`: Token value of the redis lock, ensures that only owners of a lock can release it, defaults to UUID.
- `client::Client=get_global_client()`: Redis client instance, defualts to global instance.
- `timeout=nothing`: Timeout of the lock in seconds, if `nothing` then lock will not timeout.
- `seconds_between_checks=0.1`: Sleep time (seconds) between each lock exists check.

# Examples
```julia-repl
julia> acquire_lock("example_lock", "lock_token")
"lock_token"

julia> release_lock("example_lock", "wrong_token")
false  # Returns false if the lock value does not match

julia> release_lock("example_lock", "lock_token")
true  # Returns true if the lock value matches and lock was successfully released
```
"""
function acquire_lock(lock_key, lock_value=string(uuid4()); client::Client=get_global_client(), timeout=nothing, seconds_between_checks=0.1)
    while true
        px = isnothing(timeout) ? [] : ["PX", timeout * 1000]
        was_set = !isnothing(execute(["SET", lock_key, lock_value, "NX", px...], client))

        if was_set
            return lock_value
        end

        sleep(seconds_between_checks)
    end
end

"""
    release_lock(lock_key, lock_value[; client::Client=get_global_client()])::Bool

Returns `true` if the lock value matches and lock was successfully released, `false` otherwise.

It is preferred that `redis_lock` is used instead as it will handle the acquiring and releasing of
locks within a do-block context.

# Examples
```julia-repl
julia> acquire_lock("example_lock", "lock_token")
"lock_token"

julia> release_lock("example_lock", "wrong_token")
false  # Returns false if the lock value does not match

julia> release_lock("example_lock", "lock_token")
true  # Returns true if the lock value matches and lock was successfully released
```
"""
function release_lock(lock_key, lock_value; client::Client=get_global_client())::Bool
    if get(lock_key; client=client) == lock_value
        del(lock_key; client=client)
        return true
    else
        return false
    end
end

"""
    redis_lock(fn::Function, lock_key[, lock_value=string(uuid4()); client::Client=get_global_client(), timeout=nothing, seconds_between_checks=0.1])

Enters a redis lock context, blocks if the lock already exists.

# Arguments
- `lock_key`: Name of the redis lock key.
- `lock_value=string(uuid4())`: Token value of the redis lock, ensures that only owners of a lock can release it, defaults to UUID.
- `client::Client=get_global_client()`: Redis client instance, defualts to global instance.
- `timeout=nothing`: Timeout of the lock in seconds, if `nothing` then lock will not timeout.
- `seconds_between_checks=0.1`: Sleep time (seconds) between each lock exists check.

# Examples
```julia-repl
julia> @async redis_lock("example_lock") do
           sleep(3)  # Lock will exist for 3 seconds
       end

julia> while !isredislocked("example_lock")
           sleep(0.1)  # Ensure async lock is active before proceeding
       end

julia> redis_lock("example_lock") do
           println("This message will be delayed by 3 seconds!")  # Blocked by first lock
       end
```
"""
function redis_lock(fn::Function, lock_key, lock_value=string(uuid4()); client::Client=get_global_client(), timeout=nothing, seconds_between_checks=0.1)
    acquire_lock(lock_key, lock_value; client=client, timeout=timeout, seconds_between_checks=seconds_between_checks)

    try
        fn()
    finally
        release_lock(lock_key, lock_value; client=client)
    end
end

"""
    isredislocked(lock_key[; client::Client=get_global_client()])

Returns `true` if `lock_key` exists, otherwise `false`.
"""
isredislocked(lock_key; client::Client=get_global_client()) = exists(lock_key; client=client) == 1


text:
module MyTestPackage

include("scripts.jl")
# Write your package code here.

export sayhi

end


text:
@testset "Geometry" begin

    Random.seed!(5) # reproducability

    @testset "generate_point" begin
        SAMPLES = 10000

        @testset "$geom" for geom in (Box, BoxPBC)
            lengths = [1,2,3]
            b = geom(lengths)
            testpoints = hcat([XXZNumerics.generate_point(b) for _ in 1:SAMPLES]...)

            for i in 1:length(lengths)
                @test reduce(min, testpoints[i,:]) > 0
                @test reduce(max, testpoints[i,:]) < lengths[i]
                @test abs(mean(testpoints[i,:]) - lengths[i]/2) < sqrt(SAMPLES)*lengths[i]
            end
        end

        @testset "$geom" for geom in (NoisyChain, NoisyChainPBC)
            L = 5
            # no noise
            chain = geom(L, 2, 0)
            @test sample_blockaded(chain, L) ≈ 2 .* collect(1:L)

            # some noise - but no overlaps
            chain = geom(L, 2, 1)
            testpoints = hcat([sample_blockaded(chain, L) for _ in 1:SAMPLES]...)

            for i in 1:L
                @test reduce(min, testpoints[i,:]) > 2i-1
                @test reduce(max, testpoints[i,:]) < 2i+1
                @test abs(mean(testpoints[i,:]) - 2i) < sqrt(SAMPLES)
            end
        end
    end

    @testset "distance" begin

        @testset "box" begin
            # note that box size does not matter for distances...
            box = Box([10,10,10,10])
            @test distance(box, [4,0,0,0]  , [1,0,0,0]) ≈ 3 # simple
            @test distance(box, [3,4,0,0]  , [0,0,0,0]) ≈ 5
            @test distance(box, [1,2,3,4.0], [2,1,2,5]) ≈ 2 # more complicated (and float)
            @test distance(box, [2,1,2,5], [1,2,3,4.0]) ≈ 2 # check symmetry
        end

        @testset "boxpbc" begin
            boxpbc = BoxPBC([1,2,3])
            @test distance(boxpbc, [0.4, 0.4, 0.4], [0.5, 0.5, 0.5]) ≈ sqrt(3)*0.1  # normal distance
            @test distance(boxpbc, [0.1, 0  ,   0], [0.9, 0  ,   0]) ≈ 0.2          # wrap around
            @test distance(boxpbc, [0.4, 1.9,   0], [0.7, 0.3,   0]) ≈ 0.5          # mixed case
            @test distance(boxpbc, [1  ,   1,   2], [0,     0,   0]) ≈ sqrt(2)      # more complicated
            @test distance(boxpbc, [0,     0,   0], [1  ,   1,   2]) ≈ sqrt(2)      # more complicated (check symmetry)
        end

        @testset "NoisyChain" begin
            chain = NoisyChain(5,1,0)
            p1 = XXZNumerics.generate_point(chain, 1)
            p2 = XXZNumerics.generate_point(chain, 2)
            p3 = XXZNumerics.generate_point(chain, 3)
            p4 = XXZNumerics.generate_point(chain, 4)
            @test distance(chain, p1, p2) ≈ 1
            @test distance(chain, p2, p1) ≈ 1
            @test distance(chain, p1, p3) ≈ 2
            @test distance(chain, p1, p4) ≈ 3
        end

        @testset "NoisyChainPBC" begin
            chain = NoisyChainPBC(5,1,0)
            p1 = XXZNumerics.generate_point(chain, 1)
            p2 = XXZNumerics.generate_point(chain, 2)
            p3 = XXZNumerics.generate_point(chain, 3)
            p4 = XXZNumerics.generate_point(chain, 4)
            @test distance(chain, p1, p2) ≈ 1
            @test distance(chain, p2, p1) ≈ 1
            @test distance(chain, p1, p3) ≈ 2
            @test distance(chain, p1, p4) ≈ 2 # PBC!
            @test distance(chain, p4, p1) ≈ 2 # PBC! (check symmetry)
        end
    end

    @testset "blockade" begin

        @testset "Box" begin
            b = Box([2,])
            @test_throws ErrorException sample_blockaded(b, 3)
            points = sample_blockaded(b, 2)
            p1, p2 = points
            @test distance_matrix(b, points) == [0 distance(b, p1, p2); distance(b, p1, p2) 0]
            @test distance(b, p1, p2) > 1
        end

        @testset "BoxPBC" begin
            b = BoxPBC([2.2,])
            @test_throws ErrorException sample_blockaded(b, 3)
            points = sample_blockaded(b, 2)
            p1, p2 = points
            @test distance_matrix(b, points) == [0 distance(b, p1, p2); distance(b, p1, p2) 0]
            @test distance(b, p1, p2) > 1
        end
    end

    @testset "distance_matrix" begin
        N = 10
        @testset "$geom" for geom in (Box, BoxPBC)
            @testset "$dim" for dim in 1:3
                g = geom(4*ones(dim)*N^(1/dim)) # just make the volume big enough to fit all particles in 4x the blockaded volume should be enough
                points = sample_blockaded(g, N)
                dists = distance_matrix(g, points)
                @test size(dists) == (N, N)
                @test dists ≈ distance_matrix(g, hcat(points...))
                @test all(>(1), (dists[i,j] for i in 2:N for j in i+1:N))
            end
        end
        @testset "$geom" for geom in (NoisyChain, NoisyChainPBC)
            @testset "$dim" for dim in 1:3
                g = geom(N, 1.5, 0.6) # just make the volume big enough to fit all particles in 4x the blockaded volume should be enough
                points = sample_blockaded(g, N)
                dists = distance_matrix(g, points)
                @test size(dists) == (N, N)
                @test dists ≈ distance_matrix(g, hcat(points...))
                @test all(>(1), (dists[i,j] for i in 2:N for j in i+1:N))
            end
        end
    end
end


text:
using ParticleGridMapper
using StaticArrays
#using PyPlot
using Test

const N = 3 #spatial dimension
const T = Float64
const Npart = 1000 #number of particles

using Random
Random.seed!(1114) #happy birthday! <3

println("nthreads = ", Threads.nthreads())

const BOXSIZE = 0.807 #use something != 1.0 to test round-off error

function get_spherical_cloud(Ngas, rmax)
    X = Vector{SVector{N,T}}(undef, Ngas)
    for i in 1:Ngas
        r3d = rmax * 2 #s.t. the while-loop will be evaluated at least once
        while r3d > rmax
            m_frac = rand() * 0.99;      #m_frac = 0 ~ 0.99
            r3d = rmax * m_frac;
        end
        phi       = ( 2. * rand() - 1. ) * pi
        cos_theta = ( 2. * rand() - 1. )
        x = r3d * sqrt( 1. - cos_theta^2 ) * cos(phi)
        y = r3d * sqrt( 1. - cos_theta^2 ) * sin(phi)
        z = r3d * cos_theta
        X[i] = @SVector [x,y,z]
    end
    return X
end

let
@testset "spherical cloud with AMR" begin

#randomly distributing particles
#X = [@SVector rand(N) for _ in 1:Npart]
X = get_spherical_cloud(Npart, 0.5*BOXSIZE)
topnode_length = @SVector(ones(N)) * BOXSIZE  #actual length of tree

center = @SVector zeros(N)

radius = ParticleGridMapper.norm.(X./BOXSIZE)
Mtot = 1.0
mass = ones(Npart) .* (Mtot / Npart) #particle mass
rho = 2*Mtot/(4*pi) .* radius.^(-2)

volume = mass ./ rho #only needed for SPH

#large enough s.t. there are no grids with 0 ngb particles and we can test all(image./BOXSIZE .≈  1.) (see below)
hsml = 0.7*BOXSIZE .* radius.^(2/3)

part = [DataP2G{N,T}(SVector(X[i]), i, hsml[i], mass[i], 0.0, T[]) for i in eachindex(X)]

#build the tree
treeNgb = buildtree(part, center, topnode_length);
grid_volumes = get_AMRgrid_volumes(treeNgb)
tree_max_depth = get_max_tree_depth(treeNgb)
#println( "true maximum depth of tree = ", log2( BOXSIZE / minimum(grid_volumes)^(1/3) ) )

gridAMR = get_AMRgrid(treeNgb)
num_nodes  = length(gridAMR[gridAMR.==1])
num_leaves = length(gridAMR[gridAMR.==0])
@test num_nodes + num_leaves == length(gridAMR) #gridAMR is either 0 or 1
@test Npart <= num_leaves <=  (2^N)*Npart

boxsizes = @SVector(ones(N)) * BOXSIZE  #for periodic B.C.

#map a const. field to AMR gird with NGP, SPH and MFM with different max_depth
for idepth in 1:tree_max_depth

    max_depth = idepth

    treeAMR = deepcopy(treeNgb)
    set_max_depth_AMR!(treeAMR, max_depth)
    balance_all_level!(treeAMR)

    gridAMR = get_AMRgrid(treeAMR)
    num_nodes  = length(gridAMR[gridAMR.==1])
    num_leaves = length(gridAMR[gridAMR.==0])
    @test num_nodes + num_leaves == length(gridAMR) #gridAMR is either 0 or 1

    grid_volumes = get_AMRgrid_volumes(treeAMR)
    @test sum(grid_volumes) ≈ prod(topnode_length) #particle of unity

    map_particle_to_AMRgrid_NGP!(treeAMR, ones(Npart), serial=true)
    onesAMR_NGP_s = get_AMRfield(treeAMR)
    map_particle_to_AMRgrid_NGP!(treeAMR, ones(Npart))
    onesAMR_NGP = get_AMRfield(treeAMR)
    @test all(onesAMR_NGP .≈ onesAMR_NGP_s)

    map_particle_to_AMRgrid_SPH!(treeAMR, ones(Npart), volume, X, hsml, boxsizes, treeNgb=treeNgb, serial=true)
    onesAMR_SPH_s = get_AMRfield(treeAMR)
    map_particle_to_AMRgrid_SPH!(treeAMR, ones(Npart), volume, X, hsml, boxsizes, treeNgb=treeNgb)
    onesAMR_SPH = get_AMRfield(treeAMR)
    @test all(onesAMR_SPH .≈ onesAMR_SPH_s)

    treeAMR = deepcopy(treeNgb)
    set_max_depth_AMR!(treeAMR, max_depth)
    balance_all_level!(treeAMR)
    @time for i in 1:length(X)
        map_particle_to_AMRgrid_loopP_recursive!(treeAMR, one(T), volume[i], hsml[i], X[i], boxsizes)
    end
    onesAMR_SPH_loopP = get_AMRfield(treeAMR)
    @test all(onesAMR_SPH .≈ onesAMR_SPH_loopP)

    map_particle_to_AMRgrid_MFM!(treeAMR, ones(Npart), X, hsml, boxsizes, treeNgb=treeNgb, serial=true)
    onesAMR_MFM_s = get_AMRfield(treeAMR)
    map_particle_to_AMRgrid_MFM!(treeAMR, ones(Npart), X, hsml, boxsizes, treeNgb=treeNgb)
    onesAMR_MFM = get_AMRfield(treeAMR)
    @test all(onesAMR_MFM .≈ onesAMR_MFM_s)
    @test all(onesAMR_MFM .≈ 1.) #const. field should give all ones for MFM as it is indep. of volume estimates (unless there are grids with no ngb particles)

    #this has to be right after map_particle_to_AMRgrid_MFM_thread! s.t. all(image./BOXSIZE .≈  1.)
    for p in 1:10
        nx = ny = 2^p
        image = zeros(T, nx, ny)
        image = project_AMRgrid_to_image(nx, ny, 1, 2, treeAMR, center, boxsizes, serial=true)
        image2 = zeros(T, nx, ny);
        image2 = project_AMRgrid_to_image(nx, ny, 1, 2, treeAMR, center, boxsizes)
        @test all(image .≈ image2) #parallel == serial
        @test all(image./BOXSIZE .≈  1.) #projection working properly
    end

end #idepth
end #@testset
end #let


text:
# Using user recipes from Plots.

using RecipesBase

include("utility.jl")

# Build a recipe which acts on a custom type.
@recipe function f(data::Data, var::AbstractString;
   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1)

   ndim = data.head.ndim

   if startswith(data.head.headline, "normalized")
      hasunits = false
   else
      hasunits = true
      unitw = getunit(data, var)
   end

   if ndim == 1
      VarIndex_ = findindex(data, var)
      if hasunits
         unitx = getunit(data, data.head.variables[1])
         x = data.x .* unitx
         w = data.w
         y = w[:,VarIndex_] .* unitw
      else
         x, w = data.x, data.w
         y = w[:,VarIndex_]
      end

      @series begin
         seriestype --> :path
         x, y
      end
   elseif ndim == 2
      x, y, w = getdata(data, var, plotrange, plotinterval)

      unitx = getunit(data, data.head.variables[1])
      unity = getunit(data, data.head.variables[2])

      if unitx isa UnitfulBatsrus.Unitlike
         x *= unitx
         y *= unity
         w *= unitw
      end

      @series begin
         seriestype --> :contourf  # use := if you want to force it
         x, y, w
      end
   end
end


text:
function CC.builtin_tfunction(interp::XInterpreter, @nospecialize(f), argtypes::Array{Any,1},
                           sv::Union{InferenceState,Nothing})
    ret = @invoke builtin_tfunction(interp::AbstractInterpreter, f, argtypes::Array{Any,1},
                                    sv::Union{InferenceState,Nothing})

    return ret
end

function CC.return_type_tfunc(interp::XInterpreter, argtypes::Vector{Any}, sv::InferenceState)
    ret = @invoke return_type_tfunc(interp::AbstractInterpreter, argtypes::Vector{Any}, sv::InferenceState)

    return ret
end


text:
export protein_bonds

"A list of bonded pairs representing the covalent bonds in a standard protein."
const protein_bonds = [
        # Standard residues
        BondedPair("ALA", "N", "CA"),
        BondedPair("ALA", "N", "1H"),
        BondedPair("ALA", "N", "2H"),
        BondedPair("ALA", "N", "3H"),
        BondedPair("ALA", "CA", "C"),
        BondedPair("ALA", "CA", "CB"),
        BondedPair("ALA", "CA", "HA"),
        BondedPair("ALA", "C", "O"),
        BondedPair("ALA", "CB", "1HB"),
        BondedPair("ALA", "CB", "2HB"),
        BondedPair("ALA", "CB", "3HB"),
        BondedPair("ARG", "N", "CA"),
        BondedPair("ARG", "N", "1H"),
        BondedPair("ARG", "N", "2H"),
        BondedPair("ARG", "N", "3H"),
        BondedPair("ARG", "CA", "C"),
        BondedPair("ARG", "CA", "CB"),
        BondedPair("ARG", "CA", "HA"),
        BondedPair("ARG", "C", "O"),
        BondedPair("ARG", "CB", "CG"),
        BondedPair("ARG", "CB", "1HB"),
        BondedPair("ARG", "CB", "2HB"),
        BondedPair("ARG", "CG", "CD"),
        BondedPair("ARG", "CG", "1HG"),
        BondedPair("ARG", "CG", "2HG"),
        BondedPair("ARG", "CD", "NE"),
        BondedPair("ARG", "CD", "1HD"),
        BondedPair("ARG", "CD", "2HD"),
        BondedPair("ARG", "NE", "CZ"),
        BondedPair("ARG", "NE", "HE"),
        BondedPair("ARG", "CZ", "NH1"),
        BondedPair("ARG", "CZ", "NH2"),
        BondedPair("ARG", "NH1", "1HH1"),
        BondedPair("ARG", "NH1", "2HH1"),
        BondedPair("ARG", "NH2", "1HH2"),
        BondedPair("ARG", "NH2", "2HH2"),
        BondedPair("ASN", "N", "CA"),
        BondedPair("ASN", "N", "1H"),
        BondedPair("ASN", "N", "2H"),
        BondedPair("ASN", "N", "3H"),
        BondedPair("ASN", "CA", "C"),
        BondedPair("ASN", "CA", "CB"),
        BondedPair("ASN", "CA", "HA"),
        BondedPair("ASN", "C", "O"),
        BondedPair("ASN", "CB", "CG"),
        BondedPair("ASN", "CB", "1HB"),
        BondedPair("ASN", "CB", "2HB"),
        BondedPair("ASN", "CG", "OD1"),
        BondedPair("ASN", "CG", "ND2"),
        BondedPair("ASN", "ND2", "1HD1"),
        BondedPair("ASN", "ND2", "2HD1"),
        BondedPair("ASP", "N", "CA"),
        BondedPair("ASP", "N", "1H"),
        BondedPair("ASP", "N", "2H"),
        BondedPair("ASP", "N", "3H"),
        BondedPair("ASP", "CA", "C"),
        BondedPair("ASP", "CA", "CB"),
        BondedPair("ASP", "CA", "HA"),
        BondedPair("ASP", "C", "O"),
        BondedPair("ASP", "CB", "CG"),
        BondedPair("ASP", "CB", "1HB"),
        BondedPair("ASP", "CB", "2HB"),
        BondedPair("ASP", "CG", "OD1"),
        BondedPair("ASP", "CG", "OD2"),
        BondedPair("CYS", "N", "CA"),
        BondedPair("CYS", "N", "1H"),
        BondedPair("CYS", "N", "2H"),
        BondedPair("CYS", "N", "3H"),
        BondedPair("CYS", "CA", "C"),
        BondedPair("CYS", "CA", "CB"),
        BondedPair("CYS", "CA", "HA"),
        BondedPair("CYS", "C", "O"),
        BondedPair("CYS", "CB", "SG"),
        BondedPair("CYS", "CB", "1HB"),
        BondedPair("CYS", "CB", "2HB"),
        BondedPair("CYS", "SG", "HG"),
        BondedPair("GLN", "N", "CA"),
        BondedPair("GLN", "N", "1H"),
        BondedPair("GLN", "N", "2H"),
        BondedPair("GLN", "N", "3H"),
        BondedPair("GLN", "CA", "C"),
        BondedPair("GLN", "CA", "CB"),
        BondedPair("GLN", "CA", "HA"),
        BondedPair("GLN", "C", "O"),
        BondedPair("GLN", "CB", "CG"),
        BondedPair("GLN", "CB", "1HB"),
        BondedPair("GLN", "CB", "2HB"),
        BondedPair("GLN", "CG", "CD"),
        BondedPair("GLN", "CG", "1HG"),
        BondedPair("GLN", "CG", "2HG"),
        BondedPair("GLN", "CD", "OE1"),
        BondedPair("GLN", "CD", "NE2"),
        BondedPair("GLN", "NE2", "1HE2"),
        BondedPair("GLN", "NE2", "2HE2"),
        BondedPair("GLU", "N", "CA"),
        BondedPair("GLU", "N", "1H"),
        BondedPair("GLU", "N", "2H"),
        BondedPair("GLU", "N", "3H"),
        BondedPair("GLU", "CA", "C"),
        BondedPair("GLU", "CA", "CB"),
        BondedPair("GLU", "CA", "HA"),
        BondedPair("GLU", "C", "O"),
        BondedPair("GLU", "CB", "CG"),
        BondedPair("GLU", "CB", "1HB"),
        BondedPair("GLU", "CB", "2HB"),
        BondedPair("GLU", "CG", "CD"),
        BondedPair("GLU", "CG", "1HG"),
        BondedPair("GLU", "CG", "2HG"),
        BondedPair("GLU", "CD", "OE1"),
        BondedPair("GLU", "CD", "OE2"),
        BondedPair("GLY", "N", "CA"),
        BondedPair("GLY", "N", "1H"),
        BondedPair("GLY", "N", "2H"),
        BondedPair("GLY", "N", "3H"),
        BondedPair("GLY", "CA", "C"),
        BondedPair("GLY", "CA", "1HA"),
        BondedPair("GLY", "CA", "2HA"),
        BondedPair("GLY", "C", "O"),
        BondedPair("HIS", "N", "CA"),
        BondedPair("HIS", "N", "1H"),
        BondedPair("HIS", "N", "2H"),
        BondedPair("HIS", "N", "3H"),
        BondedPair("HIS", "CA", "C"),
        BondedPair("HIS", "CA", "CB"),
        BondedPair("HIS", "CA", "HA"),
        BondedPair("HIS", "C", "O"),
        BondedPair("HIS", "CB", "CG"),
        BondedPair("HIS", "CB", "1HB"),
        BondedPair("HIS", "CB", "2HB"),
        BondedPair("HIS", "CG", "ND1"),
        BondedPair("HIS", "CG", "CD2"),
        BondedPair("HIS", "ND1", "CE1"),
        BondedPair("HIS", "ND1", "1HD1"),
        BondedPair("HIS", "CD2", "NE2"),
        BondedPair("HIS", "CD2", "1HD2"),
        BondedPair("HIS", "CE1", "NE2"),
        BondedPair("HIS", "CE1", "1HE1"),
        BondedPair("HIS", "NE2", "1HE2"),
        BondedPair("ILE", "N", "CA"),
        BondedPair("ILE", "N", "1H"),
        BondedPair("ILE", "N", "2H"),
        BondedPair("ILE", "N", "3H"),
        BondedPair("ILE", "CA", "C"),
        BondedPair("ILE", "CA", "CB"),
        BondedPair("ILE", "CA", "HA"),
        BondedPair("ILE", "C", "O"),
        BondedPair("ILE", "CB", "CG1"),
        BondedPair("ILE", "CB", "CG2"),
        BondedPair("ILE", "CB", "HB"),
        BondedPair("ILE", "CG1", "CD1"),
        BondedPair("ILE", "CG1", "1HG1"),
        BondedPair("ILE", "CG1", "2HG1"),
        BondedPair("ILE", "CG2", "1HG2"),
        BondedPair("ILE", "CG2", "2HG2"),
        BondedPair("ILE", "CG2", "3HG2"),
        BondedPair("ILE", "CD1", "1HD1"),
        BondedPair("ILE", "CD1", "2HD1"),
        BondedPair("ILE", "CD1", "3HD1"),
        BondedPair("LEU", "N", "CA"),
        BondedPair("LEU", "N", "1H"),
        BondedPair("LEU", "N", "2H"),
        BondedPair("LEU", "N", "3H"),
        BondedPair("LEU", "CA", "C"),
        BondedPair("LEU", "CA", "CB"),
        BondedPair("LEU", "CA", "HA"),
        BondedPair("LEU", "C", "O"),
        BondedPair("LEU", "CB", "CG"),
        BondedPair("LEU", "CB", "1HB"),
        BondedPair("LEU", "CB", "2HB"),
        BondedPair("LEU", "CG", "CD1"),
        BondedPair("LEU", "CG", "CD2"),
        BondedPair("LEU", "CG", "HG"),
        BondedPair("LEU", "CD1", "1HD1"),
        BondedPair("LEU", "CD1", "2HD1"),
        BondedPair("LEU", "CD1", "3HD1"),
        BondedPair("LEU", "CD2", "1HD2"),
        BondedPair("LEU", "CD2", "2HD2"),
        BondedPair("LEU", "CD2", "3HD2"),
        BondedPair("LYS", "N", "CA"),
        BondedPair("LYS", "N", "1H"),
        BondedPair("LYS", "N", "2H"),
        BondedPair("LYS", "N", "3H"),
        BondedPair("LYS", "CA", "C"),
        BondedPair("LYS", "CA", "CB"),
        BondedPair("LYS", "CA", "HA"),
        BondedPair("LYS", "C", "O"),
        BondedPair("LYS", "CB", "CG"),
        BondedPair("LYS", "CB", "1HB"),
        BondedPair("LYS", "CB", "2HB"),
        BondedPair("LYS", "CG", "CD"),
        BondedPair("LYS", "CG", "1HG"),
        BondedPair("LYS", "CG", "2HG"),
        BondedPair("LYS", "CD", "CE"),
        BondedPair("LYS", "CD", "1HD"),
        BondedPair("LYS", "CD", "2HD"),
        BondedPair("LYS", "CE", "NZ"),
        BondedPair("LYS", "CE", "1HE"),
        BondedPair("LYS", "CE", "2HE"),
        BondedPair("LYS", "NZ", "1HZ"),
        BondedPair("LYS", "NZ", "2HZ"),
        BondedPair("LYS", "NZ", "3HZ"),
        BondedPair("MET", "N", "CA"),
        BondedPair("MET", "N", "1H"),
        BondedPair("MET", "N", "2H"),
        BondedPair("MET", "N", "3H"),
        BondedPair("MET", "CA", "C"),
        BondedPair("MET", "CA", "CB"),
        BondedPair("MET", "CA", "HA"),
        BondedPair("MET", "C", "O"),
        BondedPair("MET", "CB", "CG"),
        BondedPair("MET", "CB", "1HB"),
        BondedPair("MET", "CB", "2HB"),
        BondedPair("MET", "CG", "SD"),
        BondedPair("MET", "CG", "1HG"),
        BondedPair("MET", "CG", "2HG"),
        BondedPair("MET", "SD", "CE"),
        BondedPair("MET", "CE", "1HE"),
        BondedPair("MET", "CE", "2HE"),
        BondedPair("MET", "CE", "3HE"),
        BondedPair("PHE", "N", "CA"),
        BondedPair("PHE", "N", "1H"),
        BondedPair("PHE", "N", "2H"),
        BondedPair("PHE", "N", "3H"),
        BondedPair("PHE", "CA", "C"),
        BondedPair("PHE", "CA", "CB"),
        BondedPair("PHE", "CA", "HA"),
        BondedPair("PHE", "C", "O"),
        BondedPair("PHE", "CB", "CG"),
        BondedPair("PHE", "CB", "1HB"),
        BondedPair("PHE", "CB", "2HB"),
        BondedPair("PHE", "CG", "CD1"),
        BondedPair("PHE", "CG", "CD2"),
        BondedPair("PHE", "CD1", "CE1"),
        BondedPair("PHE", "CD1", "HD1"),
        BondedPair("PHE", "CD2", "CE2"),
        BondedPair("PHE", "CD2", "HD2"),
        BondedPair("PHE", "CE1", "CZ"),
        BondedPair("PHE", "CE1", "HE1"),
        BondedPair("PHE", "CE2", "CZ"),
        BondedPair("PHE", "CE2", "HE2"),
        BondedPair("PHE", "CZ", "HZ"),
        BondedPair("PRO", "N", "CA"),
        BondedPair("PRO", "N", "CD"),
        BondedPair("PRO", "CA", "C"),
        BondedPair("PRO", "CA", "CB"),
        BondedPair("PRO", "CA", "HA"),
        BondedPair("PRO", "C", "O"),
        BondedPair("PRO", "CB", "CG"),
        BondedPair("PRO", "CB", "1HB"),
        BondedPair("PRO", "CB", "2HB"),
        BondedPair("PRO", "CG", "CD"),
        BondedPair("PRO", "CG", "1HG"),
        BondedPair("PRO", "CG", "2HG"),
        BondedPair("PRO", "CD", "1HD"),
        BondedPair("PRO", "CD", "2HD"),
        BondedPair("SER", "N", "CA"),
        BondedPair("SER", "N", "1H"),
        BondedPair("SER", "N", "2H"),
        BondedPair("SER", "N", "3H"),
        BondedPair("SER", "CA", "C"),
        BondedPair("SER", "CA", "CB"),
        BondedPair("SER", "CA", "HA"),
        BondedPair("SER", "C", "O"),
        BondedPair("SER", "CB", "OG"),
        BondedPair("SER", "CB", "1HB"),
        BondedPair("SER", "CB", "2HB"),
        BondedPair("SER", "OG", "HG"),
        BondedPair("THR", "N", "CA"),
        BondedPair("THR", "N", "1H"),
        BondedPair("THR", "N", "2H"),
        BondedPair("THR", "N", "3H"),
        BondedPair("THR", "CA", "C"),
        BondedPair("THR", "CA", "CB"),
        BondedPair("THR", "CA", "HA"),
        BondedPair("THR", "C", "O"),
        BondedPair("THR", "CB", "OG1"),
        BondedPair("THR", "CB", "CG2"),
        BondedPair("THR", "CB", "HB"),
        BondedPair("THR", "OG1", "HG1"),
        BondedPair("THR", "CG2", "1HG2"),
        BondedPair("THR", "CG2", "2HG2"),
        BondedPair("THR", "CG2", "3HG2"),
        BondedPair("TRP", "N", "CA"),
        BondedPair("TRP", "N", "1H"),
        BondedPair("TRP", "N", "2H"),
        BondedPair("TRP", "N", "3H"),
        BondedPair("TRP", "CA", "C"),
        BondedPair("TRP", "CA", "CB"),
        BondedPair("TRP", "CA", "HA"),
        BondedPair("TRP", "C", "O"),
        BondedPair("TRP", "CB", "CG"),
        BondedPair("TRP", "CB", "1HB"),
        BondedPair("TRP", "CB", "2HB"),
        BondedPair("TRP", "CG", "CD1"),
        BondedPair("TRP", "CG", "CD2"),
        BondedPair("TRP", "CD1", "NE1"),
        BondedPair("TRP", "CD1", "HD1"),
        BondedPair("TRP", "CD2", "CE2"),
        BondedPair("TRP", "CD2", "CE3"),
        BondedPair("TRP", "NE1", "CE2"),
        BondedPair("TRP", "NE1", "HE1"),
        BondedPair("TRP", "CE2", "CZ2"),
        BondedPair("TRP", "CE3", "CZ3"),
        BondedPair("TRP", "CE3", "HE3"),
        BondedPair("TRP", "CZ2", "CH2"),
        BondedPair("TRP", "CZ2", "HZ2"),
        BondedPair("TRP", "CZ3", "CH2"),
        BondedPair("TRP", "CZ3", "HZ3"),
        BondedPair("TRP", "CH2", "HH2"),
        BondedPair("TYR", "N", "CA"),
        BondedPair("TYR", "N", "1H"),
        BondedPair("TYR", "N", "2H"),
        BondedPair("TYR", "N", "3H"),
        BondedPair("TYR", "CA", "C"),
        BondedPair("TYR", "CA", "CB"),
        BondedPair("TYR", "CA", "HA"),
        BondedPair("TYR", "C", "O"),
        BondedPair("TYR", "CB", "CG"),
        BondedPair("TYR", "CB", "1HB"),
        BondedPair("TYR", "CB", "2HB"),
        BondedPair("TYR", "CG", "CD1"),
        BondedPair("TYR", "CG", "CD2"),
        BondedPair("TYR", "CD1", "CE1"),
        BondedPair("TYR", "CD1", "HD1"),
        BondedPair("TYR", "CD2", "CE2"),
        BondedPair("TYR", "CD2", "HD2"),
        BondedPair("TYR", "CE1", "CZ"),
        BondedPair("TYR", "CE1", "HE1"),
        BondedPair("TYR", "CE2", "CZ"),
        BondedPair("TYR", "CE2", "HE2"),
        BondedPair("TYR", "CZ", "OH"),
        BondedPair("TYR", "OH", "HH"),
        BondedPair("VAL", "N", "CA"),
        BondedPair("VAL", "N", "1H"),
        BondedPair("VAL", "N", "2H"),
        BondedPair("VAL", "N", "3H"),
        BondedPair("VAL", "CA", "C"),
        BondedPair("VAL", "CA", "CB"),
        BondedPair("VAL", "CA", "HA"),
        BondedPair("VAL", "C", "O"),
        BondedPair("VAL", "CB", "CG1"),
        BondedPair("VAL", "CB", "CG2"),
        BondedPair("VAL", "CB", "HB"),
        BondedPair("VAL", "CG1", "1HG1"),
        BondedPair("VAL", "CG1", "2HG1"),
        BondedPair("VAL", "CG1", "3HG1"),
        BondedPair("VAL", "CG2", "1HG2"),
        BondedPair("VAL", "CG2", "2HG2"),
        BondedPair("VAL", "CG2", "3HG2"),

        # N-H bonds
        BondedPair("ALA", "N", "H"),
        BondedPair("ARG", "N", "H"),
        BondedPair("ASN", "N", "H"),
        BondedPair("ASP", "N", "H"),
        BondedPair("CYS", "N", "H"),
        BondedPair("GLN", "N", "H"),
        BondedPair("GLU", "N", "H"),
        BondedPair("GLY", "N", "H"),
        BondedPair("HIS", "N", "H"),
        BondedPair("ILE", "N", "H"),
        BondedPair("LEU", "N", "H"),
        BondedPair("LYS", "N", "H"),
        BondedPair("MET", "N", "H"),
        BondedPair("PHE", "N", "H"),
        BondedPair("PRO", "N", "H"),
        BondedPair("SER", "N", "H"),
        BondedPair("THR", "N", "H"),
        BondedPair("TRP", "N", "H"),
        BondedPair("TYR", "N", "H"),
        BondedPair("VAL", "N", "H"),

        # Terminal oxygen
        BondedPair("ALA", "C", "OXT"),
        BondedPair("ARG", "C", "OXT"),
        BondedPair("ASN", "C", "OXT"),
        BondedPair("ASP", "C", "OXT"),
        BondedPair("CYS", "C", "OXT"),
        BondedPair("GLN", "C", "OXT"),
        BondedPair("GLU", "C", "OXT"),
        BondedPair("GLY", "C", "OXT"),
        BondedPair("HIS", "C", "OXT"),
        BondedPair("ILE", "C", "OXT"),
        BondedPair("LEU", "C", "OXT"),
        BondedPair("LYS", "C", "OXT"),
        BondedPair("MET", "C", "OXT"),
        BondedPair("PHE", "C", "OXT"),
        BondedPair("PRO", "C", "OXT"),
        BondedPair("SER", "C", "OXT"),
        BondedPair("THR", "C", "OXT"),
        BondedPair("TRP", "C", "OXT"),
        BondedPair("TYR", "C", "OXT"),
        BondedPair("VAL", "C", "OXT"),

        # A few special residues
        BondedPair("MSE", "N", "CA"),
        BondedPair("MSE", "N", "1H"),
        BondedPair("MSE", "N", "2H"),
        BondedPair("MSE", "N", "3H"),
        BondedPair("MSE", "CA", "C"),
        BondedPair("MSE", "CA", "CB"),
        BondedPair("MSE", "CA", "HA"),
        BondedPair("MSE", "C", "O"),
        BondedPair("MSE", "CB", "CG"),
        BondedPair("MSE", "CB", "1HB"),
        BondedPair("MSE", "CB", "2HB"),
        BondedPair("MSE", "CG", "1HG"),
        BondedPair("MSE", "CG", "2HG"),
        BondedPair("MSE", "CE", "1HE"),
        BondedPair("MSE", "CE", "2HE"),
        BondedPair("MSE", "CE", "3HE"),
        BondedPair("MSE", "CG", "SE"),
        BondedPair("MSE", "SE", "CE"),
        BondedPair("MSO", "N", "CA"),
        BondedPair("MSO", "N", "1H"),
        BondedPair("MSO", "N", "2H"),
        BondedPair("MSO", "N", "3H"),
        BondedPair("MSO", "CA", "C"),
        BondedPair("MSO", "CA", "CB"),
        BondedPair("MSO", "CA", "HA"),
        BondedPair("MSO", "C", "O"),
        BondedPair("MSO", "CB", "CG"),
        BondedPair("MSO", "CB", "1HB"),
        BondedPair("MSO", "CB", "2HB"),
        BondedPair("MSO", "CG", "1HG"),
        BondedPair("MSO", "CG", "2HG"),
        BondedPair("MSO", "CE", "1HE"),
        BondedPair("MSO", "CE", "2HE"),
        BondedPair("MSO", "CE", "3HE"),
        BondedPair("MSO", "CG", "SE"),
        BondedPair("MSO", "SE", "CE"),
        BondedPair("MSO", "SE", "OE"),
        BondedPair("SEP", "N", "CA"),
        BondedPair("SEP", "N", "1H"),
        BondedPair("SEP", "N", "2H"),
        BondedPair("SEP", "N", "3H"),
        BondedPair("SEP", "CA", "C"),
        BondedPair("SEP", "CA", "CB"),
        BondedPair("SEP", "CA", "HA"),
        BondedPair("SEP", "C", "O"),
        BondedPair("SEP", "CB", "OG"),
        BondedPair("SEP", "CB", "1HB"),
        BondedPair("SEP", "CB", "2HB"),
        BondedPair("SEP", "OG", "P"),
        BondedPair("SEP", "P", "O1P"),
        BondedPair("SEP", "P", "O2P"),
        BondedPair("SEP", "P", "O3P"),
        BondedPair("TPO", "N", "CA"),
        BondedPair("TPO", "N", "1H"),
        BondedPair("TPO", "N", "2H"),
        BondedPair("TPO", "N", "3H"),
        BondedPair("TPO", "CA", "C"),
        BondedPair("TPO", "CA", "CB"),
        BondedPair("TPO", "CA", "HA"),
        BondedPair("TPO", "C", "O"),
        BondedPair("TPO", "CB", "OG1"),
        BondedPair("TPO", "CB", "CG2"),
        BondedPair("TPO", "CB", "HB"),
        BondedPair("TPO", "CG2", "1HG2"),
        BondedPair("TPO", "CG2", "2HG2"),
        BondedPair("TPO", "CG2", "3HG2"),
        BondedPair("TPO", "OG1", "P"),
        BondedPair("TPO", "P", "O1P"),
        BondedPair("TPO", "P", "O2P"),
        BondedPair("TPO", "P", "O3P"),
        BondedPair("PTR", "N", "CA"),
        BondedPair("PTR", "N", "1H"),
        BondedPair("PTR", "N", "2H"),
        BondedPair("PTR", "N", "3H"),
        BondedPair("PTR", "CA", "C"),
        BondedPair("PTR", "CA", "CB"),
        BondedPair("PTR", "CA", "HA"),
        BondedPair("PTR", "C", "O"),
        BondedPair("PTR", "CB", "CG"),
        BondedPair("PTR", "CB", "1HB"),
        BondedPair("PTR", "CB", "2HB"),
        BondedPair("PTR", "CG", "CD1"),
        BondedPair("PTR", "CG", "CD2"),
        BondedPair("PTR", "CD1", "CE1"),
        BondedPair("PTR", "CD1", "HD1"),
        BondedPair("PTR", "CD2", "CE2"),
        BondedPair("PTR", "CD2", "HD2"),
        BondedPair("PTR", "CE1", "CZ"),
        BondedPair("PTR", "CE1", "HE1"),
        BondedPair("PTR", "CE2", "CZ"),
        BondedPair("PTR", "CE2", "HE2"),
        BondedPair("PTR", "CZ", "OH"),
        BondedPair("PTR", "OH", "P"),
        BondedPair("PTR", "P", "O1P"),
        BondedPair("PTR", "P", "O2P"),
        BondedPair("PTR", "P", "O3P"),
        BondedPair("KPI", "N", "CA"),
        BondedPair("KPI", "N", "1H"),
        BondedPair("KPI", "N", "2H"),
        BondedPair("KPI", "N", "3H"),
        BondedPair("KPI", "CA", "C"),
        BondedPair("KPI", "CA", "CB"),
        BondedPair("KPI", "CA", "HA"),
        BondedPair("KPI", "C", "O"),
        BondedPair("KPI", "CB", "CG"),
        BondedPair("KPI", "CB", "1HB"),
        BondedPair("KPI", "CB", "2HB"),
        BondedPair("KPI", "CG", "CD"),
        BondedPair("KPI", "CG", "1HG"),
        BondedPair("KPI", "CG", "2HG"),
        BondedPair("KPI", "CD", "CE"),
        BondedPair("KPI", "CD", "1HD"),
        BondedPair("KPI", "CD", "2HD"),
        BondedPair("KPI", "CE", "NZ"),
        BondedPair("KPI", "CE", "1HE"),
        BondedPair("KPI", "CE", "2HE"),
        BondedPair("KPI", "NZ", "CX1"),
        BondedPair("KPI", "CX1", "CX2"),
        BondedPair("KPI", "CX1", "C1"),
        BondedPair("KPI", "CX2", "O1"),
        BondedPair("KPI", "CX2", "O2"),
        BondedPair("KCX", "N", "CA"),
        BondedPair("KCX", "N", "1H"),
        BondedPair("KCX", "N", "2H"),
        BondedPair("KCX", "N", "3H"),
        BondedPair("KCX", "CA", "C"),
        BondedPair("KCX", "CA", "CB"),
        BondedPair("KCX", "CA", "HA"),
        BondedPair("KCX", "C", "O"),
        BondedPair("KCX", "CB", "CG"),
        BondedPair("KCX", "CB", "1HB"),
        BondedPair("KCX", "CB", "2HB"),
        BondedPair("KCX", "CG", "CD"),
        BondedPair("KCX", "CG", "1HG"),
        BondedPair("KCX", "CG", "2HG"),
        BondedPair("KCX", "CD", "CE"),
        BondedPair("KCX", "CD", "1HD"),
        BondedPair("KCX", "CD", "2HD"),
        BondedPair("KCX", "CE", "NZ"),
        BondedPair("KCX", "CE", "1HE"),
        BondedPair("KCX", "CE", "2HE"),
        BondedPair("KCX", "NZ", "CX"),
        BondedPair("KCX", "CX", "OQ1"),
        BondedPair("KCX", "CX", "OQ2"),
        BondedPair("CME", "N", "CA"),
        BondedPair("CME", "N", "1H"),
        BondedPair("CME", "N", "2H"),
        BondedPair("CME", "N", "3H"),
        BondedPair("CME", "CA", "C"),
        BondedPair("CME", "CA", "CB"),
        BondedPair("CME", "CA", "HA"),
        BondedPair("CME", "C", "O"),
        BondedPair("CME", "CB", "1HB"),
        BondedPair("CME", "CB", "2HB"),
        BondedPair("CME", "SG", "SD"),
        BondedPair("CME", "SD", "CE"),
        BondedPair("CME", "CE", "1HE"),
        BondedPair("CME", "CE", "2HE"),
        BondedPair("CME", "CE", "CZ"),
        BondedPair("CME", "CZ", "1HZ"),
        BondedPair("CME", "CZ", "2HZ"),
        BondedPair("CME", "CZ", "OH"),
        BondedPair("CME", "OH", "HH"),
        BondedPair("OMH", "N", "CA"),
        BondedPair("OMH", "N", "1H"),
        BondedPair("OMH", "N", "2H"),
        BondedPair("OMH", "N", "3H"),
        BondedPair("OMH", "CA", "C"),
        BondedPair("OMH", "CA", "CB"),
        BondedPair("OMH", "CA", "HA"),
        BondedPair("OMH", "C", "O"),
        BondedPair("OMH", "CB", "OG"),
        BondedPair("OMH", "CB", "1HB"),
        BondedPair("OMH", "CB", "2HB"),
        BondedPair("OMH", "OG", "P"),
        BondedPair("OMH", "P", "O1"),
        BondedPair("OMH", "P", "O2"),
        BondedPair("OMH", "P", "O3"),
        BondedPair("OMH", "O3", "C31"),
    ]


text:
using Thebes, Luxor

include(string(pathof(Thebes), "../../../data/evenmoreobjects.jl"))

evenmoreobjects = [:Apple, :Sword01, :asu, :dragon, :goblet,
:head, :heart, :helm, :house, :king, :klingon, :mushroom, :pear, :r2,
:seashell, :space_shuttle, :space_station, :teapot, :volks,
:x29_plane]

function anothergfunction(vertices, faces, labels, action=:fill)
    if !isempty(faces)
        @layer begin
            for (n, p) in enumerate(faces)
                poly(p, action)
            end
        end
    end
end

function main()
    Drawing(1200, 1200, "/tmp/evenmoreobjects.png")
    origin()
    background("grey20")
    setopacity(0.75)

    axes3D()
    eyepoint(Point3D(500, 500, 500))
    centerpoint(Point3D(0, 0, 0))
    perspective(1000)
    setline(1)

    K = 120

    m = 1
    for x in -500:K:500
        for y in -500:K:500
            o = evenmoreobjects[mod1(m, length(evenmoreobjects))]
            object = make(eval(o))

            scalefactor = distance(extrema(object.vertices)...)

            loc = Point3D(x, y, 0)
            moveby!(object, loc)

            scaleby!(object, 1/scalefactor * K, 1/scalefactor * K, 1/scalefactor * K)
            randomhue()
            setopacity(0.5)
            pin(object, gfunction= anothergfunction)
            m += 1
        end
    end
    finish()
end

main()

preview()


text:
include("api.jl")
include("generation.jl")
include("evaluation.jl")
include("spconstructs.jl")
include("util.jl")
include("creation.jl")
include("SAASolver.jl")


text:
is_diagonal_noise{uType,tType,isinplace,ND}(prob::AbstractRODEProblem{uType,tType,isinplace,ND}) = ND <: Void

isinplace{uType,tType,iip}(prob::AbstractODEProblem{uType,tType,iip}) = iip
isinplace{uType,iip}(prob::AbstractSteadyStateProblem{uType,iip}) = iip


text:
using Bonobo
using Test

using JuMP
using Cbc

const BB = Bonobo

@testset "Bonobo.jl" begin
    include("end2end/mip.jl")
end


text:
using Base.LinAlg.BLAS

function scale{T}(alpha::T, a::Matrix{T}, result::Matrix{T})
    @inbounds for j=1:size(a,2)
        for i=1:size(a,1)
            result[i,j] = alpha*a[i,j]
        end
    end
    return nothing
end

N = 4000
srand(0)
A = rand(N,N)
result = zeros(N,N)

function f0{T}(alpha::T, A::Matrix{T}, result::Matrix{T}, repeat=1)
    for i=1:repeat
        scale(alpha, A, result)
    end
    return nothing
end

@time f1(0.1, A, result)
@time f1(0.1, A, result)

A_vec = reshape(A, length(A))
result_vec = reshape(result, length(result))
N = length(A_vec)
@time BLAS.blascopy!(N, A_vec, 1, result_vec, 1)
@time BLAS.blascopy!(N, A_vec, 1, result_vec, 1)
@time BLAS.scale!(0.1, A)
@time BLAS.scale!(0.1, A)


text:
include("WaveFunctionCollapse.jl")

using ArgParse
using FileIO
using .WaveFunctionCollapse

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "filename"
            required = true
        "--patternsize"
            arg_type = Int
            default = 3
        "--width"
            arg_type = Int
            default = 32
        "--height"
            arg_type = Int
            default = 32
        "--periodic_input"
            action = :store_true
        "--periodic_output"
            action = :store_true
        "--mirror_input_horizontally"
            action = :store_true
        "--mirror_input_vertically"
            action = :store_true
        "--rotate_input_clockwise"
            action = :store_true
        "--rotate_input_anticlockwise"
            action = :store_true
        "--seed"
            arg_type = Int
            default = 0
        "--save_to_gif"
            action = :store_true
        "--fast"
            action = :store_true
    end

    return parse_args(s, as_symbols=true)
end

function main()
    parsed_args = parse_commandline()
    filename = parsed_args[:filename]
    delete!(parsed_args, :filename)
    img = generate(filename; parsed_args...)
    FileIO.save("output.png", img)
    println("Done! Saved the output image as output.png.")
end

main()


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the hyperbolic diffusion equations

equations = HyperbolicDiffusionEquations1D(nu=1.25)

"""
    initial_condition_poisson_nonperiodic(x, t, equations::HyperbolicDiffusionEquations1D)

A non-priodic harmonic function used in combination with
[`source_terms_poisson_nonperiodic`](@ref) and [`boundary_condition_poisson_nonperiodic`](@ref).

!!! note
    The only harmonic functions in 1D have the form phi(x) = A + Bx
"""
function initial_condition_harmonic_nonperiodic(x, t, equations::HyperbolicDiffusionEquations1D)
  # elliptic equation: -νΔϕ = f
  if t == 0.0
    phi = 5.0
    q1  = 0.0
  else
    A = 3
    B = exp(1)
    phi = A + B * x[1]
    q1  = B
  end
  return SVector(phi, q1)
end
initial_condition = initial_condition_harmonic_nonperiodic

boundary_conditions = BoundaryConditionDirichlet(initial_condition)

solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)

coordinates_min = -1.0
coordinates_max =  2.0
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=2,
                n_cells_max=30_000,
                periodicity=false)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    boundary_conditions=boundary_conditions,
                                    source_terms=source_terms_harmonic)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 30.0)
ode = semidiscretize(semi, tspan);

summary_callback = SummaryCallback()

resid_tol = 5.0e-12
steady_state_callback = SteadyStateCallback(abstol=resid_tol, reltol=0.0)

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

stepsize_callback = StepsizeCallback(cfl=1.75)

callbacks = CallbackSet(summary_callback, steady_state_callback,
                        analysis_callback, alive_callback,
                        save_solution,
                        stepsize_callback)

###############################################################################
# run the simulation

sol = Trixi.solve(ode, Trixi.HypDiffN3Erk3Sstar52(),
                  dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
                  save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
import StatsFuns: logsumexp

struct AnnealingEstimate
    subsample::Int
    num_temps::Int
    num_samples::Int
end

name(x::AnnealingEstimate) = "AIS"

function annealed_importance_sampling(initial, ch::MarkovChain, subsample::Int, num_temps::Int)
    weights = zeros(Float64, num_temps)
    temps = range(0, 1; length=num_temps + 1)
    acceptance = zeros(Float64, num_temps)
    acceptance[1] = 1.0

    e_prev = energy(initial, ch.system, temps[1])
    e_cur = energy(initial, ch.system, temps[2])
    weights[1] = e_prev - e_cur

    ch.θ = temps[2]
    sampler = MetropolisSampler(0, subsample, e_cur, initial, ch)
    for (i, acc) in zip(2:num_temps, sampler)
        e_prev = sampler.current_energy
        e_cur = energy(sampler.state, ch.system, temps[i + 1])

        weights[i] = weights[i - 1] + e_prev - e_cur
        acceptance[i] = acc / (subsample + 1)

        # change the temperature for the next iteration
        sampler.chain.θ = temps[i + 1]
        sampler.current_energy = e_cur
    end
    temps, weights, acceptance
end

struct AnnealingEstimationResult <: SimulationResult
    estimate::AnnealingEstimate
    inv_temps::Vector{Float64}
    weights::Array{Float64,2}
    acceptance::Array{Float64,2}
    initial_conditionals::Vector{Float64}
end
log_marginal(result::AnnealingEstimationResult) =  logmeanexp(result.weights[end, :])

function simulate(algorithm::AnnealingEstimate, initial, system; kwargs...)
    ch = chain(system; kwargs...)

    all_weights = zeros(Float64, algorithm.num_temps, algorithm.num_samples)
    acc = zeros(Float64, algorithm.num_temps, algorithm.num_samples)
    inv_temps = nothing
    initial_conditionals = zeros(algorithm.num_samples)
    for i in 1:algorithm.num_samples
        initial_energy = Inf
        signal = initial
        while isinf(initial_energy)
            signal = sample(initial, system)
            initial_energy = energy(signal, system, 1.0)
        end

        initial_conditionals[i] = initial_energy
        (temps, weights, acceptance) = annealed_importance_sampling(signal, ch, algorithm.subsample, algorithm.num_temps)
        inv_temps = temps
        all_weights[:, i] = weights
        acc[:, i] = acceptance
    end

    AnnealingEstimationResult(algorithm, collect(inv_temps), all_weights, acc, initial_conditionals)
end

function Statistics.var(result::AnnealingEstimationResult)
    max_weight = maximum(result.weights[end, :])
    log_mean_weight = max_weight + log(mean(exp.(result.weights[end,:] .- max_weight)))
    log_var = log(var(exp.(result.weights[end, :] .- max_weight))) + 2 * max_weight
    log_var -= log(size(result.weights, 2))

    exp(-2log_mean_weight + log_var)
end

function summary(results::AnnealingEstimationResult...)
    DataFrame(Weights=[r.weights[end, :] for r in results])
end


text:
"""
Update the parameters of the model using the gardients (∇) and the learning rate (η).
"""
function update_model_weights(parameters, ∇, η)
    L = trunc(Int, length(parameters) / 2)

    # update the parameters (weights and biases) for all the layers
    for l = 1:L
        parameters[string("W_", l)] = min.(max.(parameters[string("W_", l)] - η .* ∇[string("δW_", l)], -3), 3)
        # try not editign biases, and get swish working
        parameters[string("b_", l)] = min.(max.(parameters[string("b_", l)] - η .* ∇[string("δb_", l)], -3), 3)

        # Debugging printlns
        #println(∇[string("δW_", l)][1:5])
        #println(parameters[string("W_", l)][1:5])
    end # for

    return parameters
end # function

max_stars_repo_name: Squalm/NeuralNetworks

max_issues_repo_name: Squalm/NeuralNetworks

max_forks_repo_name: Squalm/NeuralNetworks


text:
#####
##### Methods for WorldView{CodeCache}
#####

function Core.Compiler.haskey(wvc::WorldView{CodeCache}, mi::MethodInstance)
    return Core.Compiler.get(wvc, mi, nothing) !== nothing
end

function Core.Compiler.get(wvc::WorldView{CodeCache}, mi::MethodInstance, default)
    cache = wvc.cache
    for ci in get!(cache.dict, mi, CodeInstance[])
        if ci.min_world <= wvc.worlds.min_world && wvc.worlds.max_world <= ci.max_world
            # TODO: if (code && (code == jl_nothing || jl_ir_flag_inferred((jl_array_t*)code)))
            return ci
        end
    end

    return default
end

function Core.Compiler.getindex(wvc::WorldView{CodeCache}, mi::MethodInstance)
    r = Core.Compiler.get(wvc, mi, nothing)
    r === nothing && throw(KeyError(mi))
    return r::CodeInstance
end

function Core.Compiler.setindex!(wvc::WorldView{CodeCache}, ci::CodeInstance,
                                 mi::MethodInstance)
    return Core.Compiler.setindex!(wvc.cache, ci, mi)
end


text:
--- test/libdl.jl.orig	2016-09-20 05:54:22.000000000 +0300
+++ test/libdl.jl
@@ -177,7 +177,7 @@ let dl = C_NULL
     end
 end

-if Sys.KERNEL in (:Linux, :FreeBSD)
+if Sys.KERNEL in (:Linux, :DragonFly, :FreeBSD)
     ccall(:jl_read_sonames, Void, ())
 end

max_stars_repo_name: danielfojt/DeltaPorts

max_issues_repo_name: danielfojt/DeltaPorts

max_forks_repo_name: danielfojt/DeltaPorts