text:
# Autogenerated wrapper script for libcgal_julia_jll for x86_64-apple-darwin
export libcgal_julia_exact, libcgal_julia_inexact

using CGAL_jll
using libcxxwrap_julia_jll
JLLWrappers.@generate_wrapper_header("libcgal_julia")
JLLWrappers.@declare_library_product(libcgal_julia_exact, "@rpath/libcgal_julia_exact.0.dylib")
JLLWrappers.@declare_library_product(libcgal_julia_inexact, "@rpath/libcgal_julia_inexact.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header(CGAL_jll, libcxxwrap_julia_jll)
    JLLWrappers.@init_library_product(
        libcgal_julia_exact,
        "lib/libcgal_julia_exact.0.17.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcgal_julia_inexact,
        "lib/libcgal_julia_inexact.0.17.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@tm1def mutable struct Chore
  Name::Union{String,Nothing}
  StartTime::Union{Dates.DateTime,Nothing}
  DSTSensitive::Union{Bool,Nothing}
  Active::Union{Bool,Nothing}
  ExecutionMode::Union{String,Nothing}
  # this might take some unpicking
  Frequency::Union{String,Nothing}
  Attributes::Union{Dict,Nothing}
end

Chore(name::AbstractString) = Chore(Dict("Name" => name))

namefield(chore::Chore) = chore.Name

@api_default function get_all_chores(api::TM1API; options...)
  result, page_data = tm1_get_paged_json(api, "Chores"; options...)
  map(Chore, get(result, "value", [])), page_data
end

@api_default function get_chore(api::TM1API, chore_name::AbstractString; options...)
  result = tm1_get_json(api, "Chores('" * chore_name * "')"; options...)
  Chore(result)
end

@api_default function delete_chore(api::TM1API, chore_name::AbstractString; options...)
  tm1_delete(api, "Chores('" * chore_name * "')"; options...)
end

ext: jl


text:
using Test
using InteractiveUtils
using MagneticReadHead: moduleof, functiontypeof

# Define an extra method of eps in this module, so we can test methods of
Base.eps(::typeof(moduleof)) = "dummy"

@testset "moduleof" begin
    for meth in methods(detect_ambiguities)
        @test moduleof(meth) == Test
    end

    # We define a verion of eps in this module
    # but we expect that it is still counted as being in `Base`
    for meth in methods(eps)
        @test moduleof(meth) == Base
    end

    for meth in methods(Vector)  # this is a UnionAll
        @test moduleof(meth) == Core
    end
end

@testset "functiontypeof" begin
    for meth in methods(sum)
        @test functiontypeof(meth) <: typeof(sum)
    end
    for meth in methods(+)  # This includes some parametric types
        @test functiontypeof(meth) <: typeof(+)
    end
    for meth in methods(detect_ambiguities)
        @test functiontypeof(meth) <: typeof(detect_ambiguities)
    end
end


text:
using QuantLib

function theVegaBumps(factorwiseBumping::Bool, marketModel::AbstractMarketModel, doCaps::Bool)
  multiplierCutoff = 50.0
  projectionTolerance = 1e-4
  numberRates = marketModel.numberOfRates
  caps = VolatilityBumpInstrumentJacobianCap[]

  if doCaps
    capStrike = marketModel.initialRates[1]
    for i = 1:numberRates - 1
      nextCap = VolatilityBumpInstrumentJacobianCap(i, i, capStrike)
      push!(caps, nextCap)
    end
  end

  swaptions = Vector{VolatilityBumpInstrumentJacobianSwaption}(undef, numberRates)

  for i in eachindex(swaptions)
    swaptions[i] = VolatilityBumpInstrumentJacobianSwaption(i, numberRates)
  end

  possibleBumps = VegaBumpCollection(marketModel, factorwiseBumping)

  bumpFinder = OrthogonalizedBumpFinder(possibleBumps, swaptions, caps, multiplierCutoff, projectionTolerance)

  theBumps = Vector{Vector{Matrix{Float64}}}()
  get_vega_bumps!(bumpFinder, theBumps)
end

function InverseFloater(rateLevel::Float64)
  numberOfRates = 20
  accrual = 0.5
  firstTime = 0.5

  strike = 0.15
  fixedMultiplier = 2.0
  floatingSpread = 0.0
  payer = true

  rateTimes = Float64[firstTime + (i-1) * accrual for i = 1:numberOfRates + 1]

  accruals = fill(accrual, numberOfRates)
  fixedStrikes = fill(strike, numberOfRates)
  floatingSpreads = fill(floatingSpread, numberOfRates)
  fixedMultipliers = fill(fixedMultiplier, numberOfRates)

  paymentTimes = Float64[firstTime + i * accrual for i = 1:numberOfRates]

  inverseFloater = MultiStepInverseFloater(rateTimes, accruals, accruals, fixedStrikes, fixedMultipliers, floatingSpreads, paymentTimes, payer)

  # exercise schedule, we can exercise on any rate time except the last one
  exerciseTimes = rateTimes[1:end-1]

  # naive exercise strategy, exercise above a trigger level
  trigger = 0.05
  swapTriggers = fill(trigger, length(exerciseTimes))
  naifStrategy = SwapRateTrigger(rateTimes, swapTriggers, exerciseTimes)

  # Longstaff-Schwartz exercise strategy
  collectedData = Vector{Vector{NodeData}}()
  basisCoefficients = Vector{Vector{Float64}}()

  # control that does nothing, need it because some control is required
  control = NothingExerciseValue(rateTimes)

  basisSystem = SwapForwardBasisSystem(rateTimes, exerciseTimes)

  # rebate that does nothing, need it because some rebate is expected
  # when you break a swap, nothing happens
  nullRebate = NothingExerciseValue(rateTimes)

  dummyProduct = CallSpecifiedMultiProduct(inverseFloater, naifStrategy, ExerciseAdapter(nullRebate))

  evolution = dummyProduct.evolution

  # params for models
  seed = 12332
  trainingPaths = 65536
  paths = 65536
  vegaPaths = 16384

  println("inverse floater")
  println("fixed strikes: ", strike)
  println("number rates: ", numberOfRates)
  println("training paths: ", trainingPaths)
  println("paths: ", paths)
  println("vega paths: ", vegaPaths)

  # set up a callibration, this would typically be done by using a callibrator
  println("rate level: ", rateLevel)

  initialNumeraireValue = 0.95
  volLevel = 0.11
  beta = 0.2
  gamma = 1.0
  numberOfFactors = min(5, numberOfRates)
  displacementLevel = 0.02

  # set up vectors
  initialRates = fill(rateLevel, numberOfRates)
  volatilities = fill(volLevel, numberOfRates)
  displacements = fill(displacementLevel, numberOfRates)

  correlations = ExponentialForwardCorrelation(rateTimes, volLevel, beta, gamma)

  calibration = FlatVol(volatilities, QuantLib.clone(correlations), evolution, numberOfFactors, initialRates, displacements)
  marketModel = QuantLib.clone(calibration)

  # use a factory
  generatorFactory = SobolBrownianGeneratorFactory(SobolDiagonalOrdering(), seed)

  numeraires = money_market_measure(evolution)
  evolver = LogNormalFwdRatePc(marketModel, generatorFactory, numeraires)
  clonedEvolver = QuantLib.clone(evolver)

  t1 = time()

  collect_node_data!(evolver, inverseFloater, basisSystem, nullRebate, control, trainingPaths, collectedData)

  t2 = time()

  # calcualte the exercise stategy's coefficients
  generic_longstaff_schwartz_regression!(collectedData, basisCoefficients)

  # turn the coefficients into an exercise strategy
  exerciseStrategy = LongstaffSchwartzExerciseStrategy(basisSystem, basisCoefficients, evolution, numeraires, nullRebate, control)

  # callable receiver swap
  callableProduct = CallSpecifiedMultiProduct(inverseFloater, exerciseStrategy, ExerciseAdapter(nullRebate))

  allProducts = MarketModelComposite()
  add_product!(allProducts, inverseFloater)
  add_product!(allProducts, callableProduct)
  finalize!(allProducts)

  accounter = AccountingEngine(clonedEvolver, QuantLib.clone(allProducts), initialNumeraireValue)

  stats = QuantLib.Math.GenericSequenceStats()
  multiple_path_values!(accounter, stats, paths)

  t3 = time()

  means = QuantLib.Math.stats_mean(stats)

  println("Means: ", means)
  println("Time to build strategy: ", t2 - t1, " seconds")
  println("Time to price: ", t3 - t2, " seconds")

  # vegas

  # do it twice once with factorwise bumping, one without
  pathsToVegas = vegaPaths

  for i = 0:3
    allowFactorwiseBumping = i % 2 > 0
    doCaps = round(Int, i / 2) > 0

    evolverEuler = LogNormalFwdRateEuler(marketModel, generatorFactory, numeraires)

    pathwiseInverseFloater = MarketModelPathwiseInverseFloater(rateTimes, accruals, accruals, fixedStrikes, fixedMultipliers, floatingSpreads, paymentTimes, payer)

    pathwiseInverseFloaterClone = QuantLib.clone(pathwiseInverseFloater)

    callableProductPathwise = CallSpecifiedPathwiseMultiProduct(pathwiseInverseFloaterClone, exerciseStrategy)

    theBumps = theVegaBumps(allowFactorwiseBumping, marketModel, doCaps)

    accountingEngineVegas = PathwiseVegasOuterAccountingEngine(QuantLib.clone(evolverEuler), QuantLib.clone(callableProductPathwise), marketModel, theBumps, initialNumeraireValue)

    values = Vector{Float64}()
    errors = Vector{Float64}()

    multiple_path_values!(accountingEngineVegas, values, errors, pathsToVegas)

    println("vega output")
    println("factorwise bumping: ", allowFactorwiseBumping)
    println("doCaps: ", doCaps)

    r = 1
    println("price estimate: ", values[r])
    r += 1

    for i = 1:numberOfRates
      println("Delta ", i, ", ", values[r], ", ", errors[r])
      r += 1
    end

    totalVega = 0.0

    for t = r:length(values)
      println("vega, ", t - 1 - numberOfRates, ", ", values[t], ", ", errors[t])
      totalVega += values[t]
    end

    println("total vega: ", totalVega)
  end

  ## TODO Upper Bound
end

function main()
  for i = 5:9
    InverseFloater(i / 100.0)
  end
end


text:
# BAM Reader
# ==========

"""
    BAM.Reader(input::IO; index=nothing)

Create a data reader of the BAM file format.

# Arguments
* `input`: data source
* `index=nothing`: filepath to a random access index (currently *bai* is supported)
"""
mutable struct Reader{T} <: BioGenerics.IO.AbstractReader
    stream::BGZFStreams.BGZFStream{T}
    header::SAM.Header
    start_offset::BGZFStreams.VirtualOffset
    refseqnames::Vector{String}
    refseqlens::Vector{Int}
    index::Union{Nothing, BAI}
end

function Base.eltype(::Type{Reader{T}}) where T
    return Record
end

function BioGenerics.IO.stream(reader::Reader)
    return reader.stream
end

function Reader(input::IO; index=nothing)
    if isa(index, AbstractString)
        index = BAI(index)
    else
        if index != nothing
            error("unrecognizable index argument")
        end
    end
    reader = init_bam_reader(input)
    reader.index = index
    return reader
end

function Base.show(io::IO, reader::Reader)
    println(io, summary(reader), ":")
      print(io, "  number of contigs: ", length(reader.refseqnames))
end

"""
    header(reader::Reader; fillSQ::Bool=false)::SAM.Header

Get the header of `reader`.

If `fillSQ` is `true`, this function fills missing "SQ" metainfo in the header.
"""
function BioGenerics.header(reader::Reader; fillSQ::Bool=false)::SAM.Header
    header = reader.header
    if fillSQ
        if !isempty(findall(reader.header, "SQ"))
            throw(ArgumentError("SAM header already has SQ records"))
        end
        header = copy(header)
        for (name, len) in zip(reader.refseqnames, reader.refseqlens)
            push!(header, SAM.MetaInfo("SQ", ["SN" => name, "LN" => len]))
        end
    end
    return header
end

# function BioGenerics.header(reader::Reader)
#     return header(reader)
# end

function Base.seek(reader::Reader, voffset::BGZFStreams.VirtualOffset)
    seek(reader.stream, voffset)
end

function Base.seekstart(reader::Reader)
    seek(reader.stream, reader.start_offset)
end

function Base.iterate(reader::Reader, rec=Record())
    if eof(reader)
        return nothing
    end
    read!(reader, rec)
    return copy(rec), rec
end

# Initialize a BAM reader by reading the header section.
function init_bam_reader(input::BGZFStreams.BGZFStream)
    # magic bytes
    B = read(input, UInt8)
    A = read(input, UInt8)
    M = read(input, UInt8)
    x = read(input, UInt8)
    if B != UInt8('B') || A != UInt8('A') || M != UInt8('M') || x != 0x01
        error("input was not a valid BAM file")
    end

    # SAM header
    textlen = read(input, Int32)
    samreader = SAM.Reader(IOBuffer(read(input, textlen)))

    # reference sequences
    refseqnames = String[]
    refseqlens = Int[]
    n_refs = read(input, Int32)
    for _ in 1:n_refs
        namelen = read(input, Int32)
        data = read(input, namelen)
        seqname = unsafe_string(pointer(data))
        seqlen = read(input, Int32)
        push!(refseqnames, seqname)
        push!(refseqlens, seqlen)
    end

    voffset = isa(input.io, Base.AbstractPipe) ?
        BGZFStreams.VirtualOffset(0, 0) :
        BGZFStreams.virtualoffset(input)
    return Reader(
        input,
        samreader.header,
        voffset,
        refseqnames,
        refseqlens,
        nothing)
end

function init_bam_reader(input::IO)
    return init_bam_reader(BGZFStreams.BGZFStream(input))
end

function _read!(reader::Reader, record)
    unsafe_read(
        reader.stream,
        pointer_from_objref(record),
        FIXED_FIELDS_BYTES)
    dsize = data_size(record)
    if length(record.data) < dsize
        resize!(record.data, dsize)
    end
    unsafe_read(reader.stream, pointer(record.data), dsize)
    record.reader = reader
    return record
end


text:
# ---
# title: 1542. Find Longest Awesome Substring
# id: problem1542
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: String, Bit Manipulation
# link: <https://leetcode.com/problems/find-longest-awesome-substring/description/>
# hidden: true
# ---
#
# Given a string `s`. An _awesome_ substring is a non-empty substring of `s`
# such that we can make any number of swaps in order to make it palindrome.
#
# Return the length of the maximum length **awesome substring** of `s`.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "3242415"
#     Output: 5
#     Explanation: "24241" is the longest awesome substring, we can form the palindrome "24142" with some swaps.
#
#
# **Example 2:**
#
#
#
#     Input: s = "12345678"
#     Output: 1
#
#
# **Example 3:**
#
#
#
#     Input: s = "213123"
#     Output: 6
#     Explanation: "213123" is the longest awesome substring, we can form the palindrome "231132" with some swaps.
#
#
# **Example 4:**
#
#
#
#     Input: s = "00"
#     Output: 2
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 10^5`
#   * `s` consists only of digits.
#
#
## @lc code=start
using LeetCode

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


text:
# ---
# title: 743. Network Delay Time
# id: problem743
# author: Indigo
# date: 2021-07-01
# difficulty: Medium
# categories: Heap, Depth-first Search, Breadth-first Search, Graph
# link: <https://leetcode.com/problems/network-delay-time/description/>
# hidden: true
# ---
#
# There are `N` network nodes, labelled `1` to `N`.
#
# Given `times`, a list of travel times as **directed** edges `times[i] = (u, v,
# w)`, where `u` is the source node, `v` is the target node, and `w` is the time
# it takes for a signal to travel from source to target.
#
# Now, we send a signal from a certain node `K`. How long will it take for all
# nodes to receive the signal? If it is impossible, return `-1`.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2019/05/23/931_example_1.png)
#
#
#
#     Input: times = [[2,1,1],[2,3,1],[3,4,1]], N = 4, K = 2
#     Output: 2
#
#
#
#
# **Note:**
#
#   1. `N` will be in the range `[1, 100]`.
#   2. `K` will be in the range `[1, N]`.
#   3. The length of `times` will be in the range `[1, 6000]`.
#   4. All edges `times[i] = (u, v, w)` will have `1 <= u, v <= N` and `0 <= w <= 100`.
#
#
## @lc code=start
using LeetCode

function network_delay_time(times::Vector{Vector{Int}}, n::Int, k::Int)
    graph = [Tuple{Int,Int}[] for _ in 1:n]
    for (u, v, t) in times
        push!(graph[u], (v, t))
    end
    pq = PriorityQueue{Int,Int}([(k, 0)])
    res = Int[]
    while !isempty(pq)
        node, d = dequeue_pair!(pq)
        for (neb, d1) in graph[node]
            pq[neb] = min(get(pq, neb, typemax(Int)), d + d1)
        end
        push!(res, d)
    end
    return length(res) == n ? maximum(res) : -1
end
## @lc code=end


text:
using HyperdimensionalComputing
using Test

include("vectors.jl")
include("operations.jl")
include("encoding.jl")
include("inference.jl")


text:
@testset "Counters" begin

    @testset "count and probabilities" begin
        seq = res"ARNDCQEGHILKMFPSTWYV-"

        for alphabet in (UngappedAlphabet(),
                         GappedAlphabet(),
                         ReducedAlphabet("(AILMV)(NQST)(RHK)(DE)(FWY)CGP"))
            for N in 1:3
                seqs = ((seq for i in 1:N)...,)::NTuple{N,Vector{Residue}}

                table = ContingencyTable(Float64, Val{N}, alphabet) # zeros in MIToS 1.0
                if N == 1
                    @test table[1] == 0.0
                elseif N == 2
                    @test table[1,1] == 0.0
                else
                    @test table[1,1,1] == 0.0
                end
                @test getmarginals(table)[1,1] == 0.0
                @test gettotal(table) == 0.0

                count!(table, NoClustering(), NoPseudocount(), seqs...)

                @test table == count(seqs..., alphabet=alphabet)

                if isa(alphabet, ReducedAlphabet)
                    @test table[1] == 5.0
                    @test getmarginals(table)[1] == 5.0
                    @test gettotal(table) ≈ 20.0 # Reduced alphabet without gap
                else
                    @test table[1] == 1.0
                    @test getmarginals(table)[1] == 1.0
                    @test gettotal(table) ≈ length(alphabet)
                    if N == 2
                        len = length(alphabet)
                        @test gettablearray(table) == Matrix{Float64}(I, len, len)
                        @test getmarginalsarray(table)[:,1] == [1.0 for i in 1:len]
                    end
                end

                normalize!(table)
                @test table == probabilities(seqs..., alphabet=alphabet)
            end
        end

        @testset "Using clustering" begin

            clusters = Clusters([21], ones(Int,21), Weights(ones(21)/21))

            for alphabet in (UngappedAlphabet(),
                            GappedAlphabet(),
                            ReducedAlphabet("(AILMV)(NQST)(RHK)(DE)(FWY)CGP"))

                for N in 1:3
                    seqs = ((seq for i in 1:N)...,)::NTuple{N,Vector{Residue}}

                    table = ContingencyTable(Float64, Val{N}, alphabet)

                    count!(table, clusters, NoPseudocount(), seqs...)

                    @test table == count(seqs..., alphabet=alphabet, weights=clusters)

                    len = length(alphabet)
                    if isa(alphabet, ReducedAlphabet)
                        @test table[1] == 5.0/21
                        @test getmarginals(table)[1] == 5.0/21
                        @test gettotal(table) ≈ (1.0/21) * 20.0 # ReducedAlphabet without gap
                    else
                        @test table[1] == 1.0/21
                        @test getmarginals(table)[1] == 1.0/21
                        @test gettotal(table) ≈ (1.0/21) * len
                        if N == 2
                            @test gettablearray(table) == (1.0/21) .* Matrix{Float64}(I, len, len)
                            @test getmarginalsarray(table)[:,1] == [1.0/21 for i in 1:len]
                        end
                    end

                    normalize!(table)
                    @test table == probabilities(seqs..., alphabet=alphabet, weights=clusters)
                end
            end
        end

        @testset "Using pseudocount" begin

            for alphabet in (UngappedAlphabet(),
                            GappedAlphabet(),
                            ReducedAlphabet("(AILMV)(NQST)(RHK)(DE)(FWY)CGP"))

                for N in 1:3
                    seqs = ((seq for i in 1:N)...,)::NTuple{N,Vector{Residue}}

                    table = ContingencyTable(Float64, Val{N}, alphabet)

                    count!(table, NoClustering(), AdditiveSmoothing(1.0), seqs...)

                    @test table == count(seqs..., alphabet = alphabet,
                                                  pseudocounts = AdditiveSmoothing(1.0))

                    len = Float64(length(alphabet))
                    if isa(alphabet, ReducedAlphabet)
                        @test table[1] == 5.0 + 1.0
                        @test getmarginals(table)[1] == 5. + (N==1 ? 1. : N==2 ? len : len^2)
                        @test gettotal(table) ≈ length(gettable(table)) + 20.0 # without gap
                    else
                        @test table[1] == 2.0
                        @test getmarginals(table)[1] == 1. + (N==1 ? 1. : N==2 ? len : len^2)
                        @test gettotal(table) ≈ len + length(gettable(table))
                        if N == 2
                            @test gettablearray(table) == Matrix{Float64}(I, Int(len), Int(len)) .+ 1.0
                            @test getmarginalsarray(table)[:,1] == [ len + 1.0 for i in 1:len ]
                        end
                    end

                    normalize!(table)
                    @test table == probabilities(seqs..., alphabet = alphabet,
                                                 pseudocounts = AdditiveSmoothing(1.0))
                end
            end
        end

        @testset "pseudofrequencies" begin

            table = ContingencyTable(Float64, Val{2}, UngappedAlphabet())

            probabilities!(table, NoClustering(), NoPseudocount(),
                                  BLOSUM_Pseudofrequencies(1.,0.), seq, seq)

            @test table == probabilities(seq, seq,
                                         pseudofrequencies = BLOSUM_Pseudofrequencies(1.,0.))
            @test table == probabilities(seq, seq)

            @test table[1] == 1.0/20.
            @test getmarginals(table)[1] == 1.0/20.
            @test gettotal(table) ≈ 1.0
            @test gettablearray(table) == Matrix{Float64}(I, 20, 20) ./ 20.0

            @test table != probabilities(seq, seq,
                                         pseudofrequencies = BLOSUM_Pseudofrequencies(0.,1.))
        end

        @testset "BigFloat" begin
            table = ContingencyTable(BigFloat, Val{2}, UngappedAlphabet())
            clusters = Clusters([21], ones(Int,21), Weights(ones(21)/21))
            probabilities!(table, clusters, AdditiveSmoothing(one(BigFloat)),
                           BLOSUM_Pseudofrequencies(1.0,1.0), seq, seq)

            @test sum(table) ≈ one(BigFloat)
            @test eltype(table) == BigFloat
            @test isa(table[1], BigFloat)
        end

        @testset "delete_dimensions" begin
            Pxyz = probabilities(seq, seq, seq)
            @test delete_dimensions(Pxyz, 3) == probabilities(seq, seq)
            @test delete_dimensions(Pxyz, 3, 2) == probabilities(seq)
            Pxy = delete_dimensions(Pxyz, 3)
            @test delete_dimensions!(Pxy, Pxyz, 1) == probabilities(seq, seq)
            @test sum(Pxy) ≈ 1.0
            Nxyz = count(seq, seq, seq)
            Nxy = delete_dimensions(Nxyz, 3)
            @test Nxy == count(seq, seq)
            @test delete_dimensions(Nxyz, 3, 2) == count(seq)
            @test delete_dimensions!(Nxy, Nxyz, 1) == count(seq, seq)
        end
    end
end


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

# This code was ported from a MATLAB implementation available at
# http://exolete.com/code/lbm/
# with permission of the author.
# Original MATLAB implementation is copyright (c) 2006 Iain Haslam.

using DocOpt

function boltzmann_kernel(nx, ny, numactivenodes, F, ON, UX, UY)
    avu = 1.
    prevavu = 1.
    ts = 0
    deltaU = 1e-7
    t1 = 4/9
    t2 = 1/9
    t3 = 1/36
    c_squ = 1/3
    omega = 1.0
    while (ts<4000 && 1e-10<abs((prevavu-avu)/avu)) || ts<100
        # Propagate
        F[4]=F[4][[2:nx; 1], [ny; 1:ny-1]]
        F[3]=F[3][:, [ny; 1:ny-1]]
        F[2]=F[2][[nx; 1:nx-1], [ny; 1:ny-1]]
        F[5]=F[5][[2:nx; 1], :]
        F[1]=F[1][[nx; 1:nx-1], :]
        F[6]=F[6][[2:nx; 1], [2:ny; 1]]
        F[7]=F[7][:, [2:ny; 1]]
        F[8]=F[8][[nx; 1:nx-1], [2:ny; 1]]
        DENSITY=sum(F)
        UX=(sum(F[[1;2;8]])-sum(F[[4;5;6]]))./DENSITY
        UY=(sum(F[[2;3;4]])-sum(F[[6;7;8]]))./DENSITY
        UX[1,1:ny]=UX[1,1:ny].+deltaU # Increase inlet pressure
        UX[ON]=0.
        UY[ON]=0.
        DENSITY[ON]=0.
        U_SQU=UX.^2+UY.^2
        U_C2=UX+UY
        U_C4=-UX+UY
        U_C6=-U_C2
        U_C8=-U_C4
        FEQ = Array(Any, 9)
        # Calculate equilibrium distribution: stationary
        FEQ[9]=t1.*DENSITY.*(1.-U_SQU./(2*c_squ))
        # nearest-neighbours
        FEQ[1]=t2.*DENSITY.*(1.+UX/c_squ+0.5*(UX./c_squ).^2-U_SQU./(2*c_squ))
        FEQ[3]=t2.*DENSITY.*(1.+UY/c_squ+0.5*(UY./c_squ).^2-U_SQU./(2*c_squ))
        FEQ[5]=t2.*DENSITY.*(1.-UX/c_squ+0.5*(UX./c_squ).^2-U_SQU./(2*c_squ))
        FEQ[7]=t2.*DENSITY.*(1.-UY/c_squ+0.5*(UY./c_squ).^2-U_SQU./(2*c_squ))
        # next-nearest neighbours
        FEQ[2]=t3.*DENSITY.*(1.+U_C2/c_squ+0.5*(U_C2./c_squ).^2-U_SQU./(2*c_squ))
        FEQ[4]=t3.*DENSITY.*(1.+U_C4/c_squ+0.5*(U_C4./c_squ).^2-U_SQU./(2*c_squ))
        FEQ[6]=t3.*DENSITY.*(1.+U_C6/c_squ+0.5*(U_C6./c_squ).^2-U_SQU./(2*c_squ))
        FEQ[8]=t3.*DENSITY.*(1.+U_C8/c_squ+0.5*(U_C8./c_squ).^2-U_SQU./(2*c_squ))
        FEQ = omega .* FEQ .+ (1. - omega) .* F
        FEQ[1][ON] .= F[5][ON]
        FEQ[2][ON] .= F[6][ON]
        FEQ[3][ON] .= F[7][ON]
        FEQ[4][ON] .= F[8][ON]
        FEQ[5][ON] .= F[1][ON]
        FEQ[6][ON] .= F[2][ON]
        FEQ[7][ON] .= F[3][ON]
        FEQ[8][ON] .= F[4][ON]
        F = FEQ
        prevavu = avu
        avu = sum(UX) / numactivenodes
        ts = ts + 1
    end
    return F
end

function boltzmann(nx, ny)
    density = 1.0
    F  = [ fill!(zeros(nx, ny), density/9) for i = 1:9 ]
    UX = zeros(nx, ny)
    UY = zeros(nx, ny)
    ON = rand(nx, ny) .> 0.7 # extremely porous random domain
    println("length of ON = ", sum(ON))
    OFF = ~ON
    numactivenodes = sum(OFF)
    boltzmann_kernel(nx, ny, numactivenodes, F, ON, UX, UY)
end

function main()
    doc = """boltzmann.jl

2D lattice Boltzmann fluid flow model.

Usage:
  boltzmann.jl -h | --help
  boltzmann.jl [--nx=<nx>] [--ny=<ny>]

Options:
  -h --help  Show this screen.
  --nx=<nx>  Specify size of simulation grid in the x dimension [default: 200].
  --ny=<ny>  Specify size of simulation grid in the y dimension [default: 200].
"""
    arguments = docopt(doc)
    nx = parse(Int, arguments["--nx"])
    ny = parse(Int, arguments["--ny"])

    srand(0)

    println("grid size: ", nx, " x ", ny)

    tic()
    F = boltzmann(1, 1)
    println("SELFPRIMED ", toq())

    tic()
    F = boltzmann(nx, ny)
    println("SELFTIMED ", toq())

    println("checksum: ", sum(F[1]), " ", sum(F[2]), " ", sum(F[3]), " ", sum(F[4]), " ", sum(F[5]), " ", sum(F[6]), " ", sum(F[7]), " ", sum(F[8]), " ", sum(F[9]))
end

main()


text:
const _V_OPEN = "["
const _H_OPEN = "{"
_is_open(s) = s == _V_OPEN || s == _H_OPEN

# We have a non-capturing and a capturing version of each regex.
# The non-capturing are used for checking validity, the capturing are
# used for parsing. Ideally we should only use the capturing version
# and cache the result but this is not performance-sensitive code.
# So our trade-off here is just not doing the capturing (slower) version
# two times for each token.
const _PIECE_REGEX_NC = r"^[1-9]\d*p[1-9]\d*x[1-9]\d*$"
const _PIECE_REGEX_C = r"^([1-9]\d*)p([1-9]\d*)x([1-9]\d*)$"
const _PLATE_REGEX_NC = Regex(
	"^P[1-9]\\d*x[1-9]\\d*[$(_V_OPEN)$(_H_OPEN)]?\$"
)
const _PLATE_REGEX_C = Regex(
	"^P([1-9]\\d*)x([1-9]\\d*)([$(_V_OPEN)$(_H_OPEN)]?)\$"
)

const _V_CLOSE = "]"
const _H_CLOSE = "}"
_is_close(s) = s == _V_CLOSE || s == _H_CLOSE

_is_valid_piece(s) = occursin(_PIECE_REGEX_NC, s)
_is_valid_plate(s) = occursin(_PLATE_REGEX_NC, s)
# The valid tokens do not include _is_open because these characters are
# part of a valid plate token.
_is_valid_token(s) = _is_close(s) || _is_valid_piece(s) || _is_valid_plate(s)

# Only guarantee correctness and not throwing an exception if the token
# returns true for `_is_valid_token`
_is_piece(token) = isdigit(first(token))
_is_plate(token) = first(token) === 'P'

function _parse_piece(token, ::Type{D}, ::Type{S}) where {D, S}
	m = match(_PIECE_REGEX_C, token)
	i, l, w =
		parse(D, m.captures[1]), parse(S, m.captures[2]), parse(S, m.captures[3])

	return CutPattern(l, w, i)
end

function _are_matching_brackets(open, close)
	open == "[" && close == "]" && return true
	open == "{" && close == "}" && return true
	return false
end

function _parse_plate!(
	token, pattern_stack :: Vector{Vector{CutPattern{D, S}}}, close_stack
) where {D, S}
	m = match(_PLATE_REGEX_C, token)
	L, W = parse(S, m.captures[1]), parse(S, m.captures[2])
	open = m.captures[3]
	if isempty(open) # i.e., the plate is waste
		subpatterns = CutPattern{D, S}[]
		is_vertical = true
	else # i.e., the plate has subpatterns
		subpatterns = pop!(pattern_stack)
		@assert _is_open(open)
		close = pop!(close_stack)
		if !_are_matching_brackets(open, close)
			error(
				"A subpattern started by a $(open) bracket ends with a" *
				" $(close) bracket.")
		end
		is_vertical = open == _V_OPEN
	end

	return CutPattern(L, W, is_vertical, subpatterns)
end

function parse_pattern(
	str :: String, ::Type{D}, ::Type{S}
) :: Vector{CutPattern{D, S}} where {D, S}
	tokens = split(str)

	# Note: pattern_stack starts with an empty level, while close_stack starts
	# completely empty.
	pattern_stack = Vector{CutPattern{D, S}}[CutPattern{D, S}[]]
	close_stack = String[]

	for token in reverse(tokens)
		if !_is_valid_token(token)
			error("Token \"$token\" was not recognized as a valid token.")
		end

		if _is_close(token)
			push!(pattern_stack, CutPattern{D, S}[])
			push!(close_stack, token)
			continue
		end

		if _is_piece(token)
			push!(last(pattern_stack), _parse_piece(token, D, S))
			continue
		end

		@assert _is_plate(token)

		# A plate may have subpatterns or just be waste. If it is waste,
		# then pattern_stack and close_stack are left untouched. If it has
		# subpatterns then it check the correctness by popping close_stack
		# and it gets its subpatterns by popping pattern_stack.
		plate = _parse_plate!(token, pattern_stack, close_stack)
		push!(last(pattern_stack), plate)
	end

	if !isempty(close_stack)
		error("The pattern finished without closing all subpatterns.")
	end

	return only(pattern_stack)
end


text:
using BenchmarkTools, Shuffle

const SUITE = BenchmarkGroup()

const DECK = collect(1:52)

## Faro ##

SUITE["Faro"] = BenchmarkGroup(["deterministic"])
SUITE["Faro"]["in"] = BenchmarkGroup()
SUITE["Faro"]["out"] = BenchmarkGroup()

SUITE["Faro"]["in"]["shuffle"] = @benchmarkable shuffle($DECK, Faro{:in}())

SUITE["Faro"]["in"]["shuffle!"] = @benchmarkable shuffle!($DECK, Faro{:in}())

SUITE["Faro"]["in"]["nshuffle"] = @benchmarkable nshuffle($DECK, 26, Faro{:in}())

SUITE["Faro"]["in"]["nshuffle!"] = @benchmarkable nshuffle!($DECK, 26, Faro{:in}())

SUITE["Faro"]["out"]["shuffle"] = @benchmarkable shuffle($DECK, Faro{:out}())

SUITE["Faro"]["out"]["shuffle!"] = @benchmarkable shuffle!($DECK, Faro{:out}())

SUITE["Faro"]["out"]["nshuffle"] = @benchmarkable nshuffle($DECK, 26, Faro{:out}())

SUITE["Faro"]["out"]["nshuffle!"] = @benchmarkable nshuffle!($DECK, 26, Faro{:out}())

## Cut ##

SUITE["Cut"] = BenchmarkGroup(["deterministic"])

SUITE["Cut"]["shuffle"] = @benchmarkable shuffle($DECK, Cut(n)) setup=(n=rand(0:52))

SUITE["Cut"]["nshuffle"] = @benchmarkable nshuffle($DECK, 10, Cut(n)) setup=(n=rand(0:52))

## RandomShuffle ##

SUITE["Random"] = BenchmarkGroup(["random"])

SUITE["Random"]["shuffle"] = @benchmarkable shuffle($DECK, RandomShuffle())

SUITE["Random"]["shuffle!"] = @benchmarkable shuffle!($DECK, RandomShuffle())

## GilbertShannonReeds ##

SUITE["GilbertShannonReeds"] = BenchmarkGroup(["random"])

SUITE["GilbertShannonReeds"]["shuffle"] =
    @benchmarkable shuffle($DECK, GilbertShannonReeds())

SUITE["GilbertShannonReeds"]["nshuffle"] =
    @benchmarkable nshuffle($DECK, 7, GilbertShannonReeds())

SUITE["GilbertShannonReeds"]["nshuffle!"] =
    @benchmarkable nshuffle!($DECK, 7, GilbertShannonReeds())


text:
module TestFactCheck

using FactCheck

# We have to define macros around certain functions that return expressions
# so that the :escape Exprs inside them are expanded.
#
# See: https://gist.github.com/zachallaun/5257930
#
macro throws_pred(ex) FactCheck.throws_pred(ex) end
macro fact_pred(x, y) FactCheck.fact_pred(x, y) end

type Foo a end
type Bar a end
==(x::Foo, y::Foo) = x.a == y.a

facts("FactCheck core functions") do

    context("throws_pred is true on error") do
        @fact @throws_pred(1 + 1)   => (false, "no error")
        @fact @throws_pred(error()) => (true, "error")
    end

    context("fact_pred tests equality on values") do
        @fact @fact_pred(1, 1) => (true, 1)
        @fact @fact_pred("foo", "foo") => (true, "foo")
        @fact @fact_pred(Foo(1), Foo(1)) => (x) -> x[1]
        @fact @fact_pred(Bar(1), Bar(1)) => (x) -> !x[1]
    end

    context("fact_pred applies Function assertions") do
        @fact @fact_pred(2, iseven) => (true, 2)
        @fact @fact_pred(2, isodd)  => (false, 2)

        isone(x) = x == 1
        @fact @fact_pred(1, isone)  => (true, 1)
    end

end

type Baz end
type Bazz a end

facts("FactCheck assertion helper functions") do

    context("`not` works for values and functions") do
        notone = not(1)
        @fact notone(2) => true
        @fact notone(1) => false

        noteven = not(iseven)
        @fact noteven(3) => true
        @fact noteven(2) => false
        @fact not(iseven)(2) => false
    end

    context("`truthy` is anything other than nothing or false (which is 0)") do
        @fact truthy(-1) => true
        @fact truthy("") => true
        @fact truthy([]) => true
        @fact truthy(Dict())  => true
        @fact truthy(nothing) => false
        @fact truthy(false)   => false
        @fact truthy(0)       => false
    end

    context("`anything` and `irrelevant` are always true") do
        @fact anything => exactly(irrelevant)
        @fact anything(false)   => true
        @fact irrelevant(false) => true
    end

    context("`exactly` can be used to check object equality") do
        @fact exactly(exactly)(exactly) => true

        x() = ()
        @fact exactly(x)(x) => true

        # types with no fields return a singleton object when instantiated
        @fact exactly(Baz())(Baz()) => true

        @fact exactly(Bazz(1))(Bazz(1)) => false
    end

    context("`roughly` compares numbers... roughly") do
        @fact 2.4999999999999 => roughly(2.5) #roughly(2.5)(2.4999) => true
        @fact 9.5 => roughly(10; atol=1.0) #roughly(10, 1)(9.5) => true
        @fact 10.5 => roughly(10; atol=1.0) #roughly(10, 1)(10.5) => true
    end

    context("`roughly` compares matrixes... roughly") do
        X = [1.1 1.2; 2.1 2.2]
        Y = X + [0 0.000001; -0.00000349 0.00001]
        Z = [1 1; 2 2]
        @fact X => roughly(Y)
        @fact X => roughly(Z; atol=0.2)
    end

end

exitstatus()

end # module


text:
export getSensTMatVec

function getSensTMatVec(v::Vector,m::Vector,pFor::FWIparam)

    # extract pointers
    M  				= pFor.Mesh
    omega 			= pFor.omega
    gamma 			= pFor.gamma
    Q     			= pFor.Sources
    P     			= pFor.Receivers
    Ainv			= pFor.ForwardSolver

	batchSize 		= pFor.forwardSolveBatchSize;
	select  		= pFor.sourceSelection;

	nsrc = size(Q,2); nrec = size(P,2);

	if length(select) > 0
		nsrc = length(select);
	end

	if batchSize > nsrc
		batchSize = nsrc;
	end

	numBatches 	= ceil(Int64,nsrc/batchSize);
	n = prod(M.n.+1);

	H = spzeros(ComplexF64,n,n);
	if isa(Ainv,ShiftedLaplacianMultigridSolver)
		H = GetHelmholtzOperator(M,m,omega, gamma, true,useSommerfeldBC);
		Ainv.helmParam = HelmholtzParam(M,gamma,m,omega,true,useSommerfeldBC);
		H = H + GetHelmholtzShiftOP(m, omega,Ainv.shift[1]);
		H = sparse(H');
		# H is actually shifted laplacian now...
	elseif isa(Ainv,JuliaSolver)
		H = GetHelmholtzOperator(M,m,omega, gamma, true,useSommerfeldBC);
	end

	JTv = zeros(Float64,n)
	Vdatashape = reshape(v,nrec,nsrc);

	for k_batch = 1:numBatches
		batchIdxs = (k_batch-1)*batchSize + 1 : min(k_batch*batchSize,nsrc);
		V = P*Vdatashape[:,batchIdxs];
		V,Ainv = solveLinearSystem(H,V,Ainv,1);
		if pFor.useFilesForFields
			filename = getFieldsFileName(omega);
			file     = matopen(filename);
			U        = read(file,string("Ubatch_",k_batch));
			close(file);
		else
			U = pFor.Fields[:,batchIdxs];
		end
		# V = conj((1.0.+1im*vec(gamma)).*U).*V;
		V = conj((1.0.-1im*vec(gamma./omega)).*U).*V;

		JTv .+= omega^2*vec(real(sum(V,dims=2)));
	end

	if isa(Ainv,ShiftedLaplacianMultigridSolver)
		# println("Clearing!!!");
		clear!(Ainv.MG);
	end

	# V = P*Vdatashape;
	# V,Ainv = solveLinearSystem(H,V,Ainv,1)   # Lam = ForwardSolver\(P*V);

	# # JTv    = conj((1-1im*vec(gamma./omega)).*U).*V;
	# # V      = 0;
	# # JTv	 = omega^2*vec(sum(real(JTv),2));

    return JTv
end


text:
import Rimu.BitStringAddresses: update_component
function update_component(fs::CompositeFS{2}, new_add_a, new_add_b, ::Val{1}, ::Val{2})
    return typeof(fs)(new_add_a, new_add_b)
end
function update_component(fs::CompositeFS{2}, new_add_a, new_add_b, ::Val{2}, ::Val{1})
    return typeof(fs)(new_add_b, new_add_a)
end
function update_component(
    fs::CompositeFS, new_add_a, new_add_b, ::Val{A}, ::Val{B}
) where {A,B}
    fs = update_component(fs, new_add_a, Val(A))
    return update_component(fs, new_add_b, Val(B))
end

"""
    ConstFunction(x)

Construct a function that returns `x` for all arguments.
"""
struct ConstFunction{T}
    value::T
end
(c::ConstFunction)(args...; kwargs...) = c.value

struct InteractionMatrixFunction{N,T,M<:SMatrix{N,N,T}}
    values::M
end
(imf::InteractionMatrixFunction)(i, j, args...) = imf.values[i, j]

function to_parameter_vector(address, t::Number)
    C = num_components(address)
    return SVector{C}(fill(float(t), C))
end
function to_parameter_vector(address, t)
    C = num_components(address)
    length(t) == C || throw(ArgumentError("`t` must be a $C-vector or a scalar"))
    return SVector{C}(float.(t))
end

struct ParameterColumnFunction{N,T}
    values::SVector{N,T}
end
function ParameterColumnFunction(address, t)
    return ParameterColumnFunction(to_parameter_vector(address, t))
end
(pcf::ParameterColumnFunction)(i, args...) = pcf.values[i]

struct KineticEnergyFunction{K,M}
    t::SVector{K,Float64}
    kes::SVector{M,Float64}
end

function KineticEnergyFunction(address, t, dispersion)
    M = num_modes(address)
    step = 2π/M
    if isodd(M)
        start = -π*(1+1/M) + step
    else
        start = -π + step
    end
    kr = range(start; step=step, length=M)
    ks = SVector{M}(kr)
    kes = SVector{M}(dispersion.(kr))
    return KineticEnergyFunction(to_parameter_vector(address, t), kes)
end
(ke::KineticEnergyFunction)(σ, p) = ke.kes[p] * ke.t[σ]


text:
include("test_abFunc.jl")
include("test_abGrad.jl")
include("test_bFunc.jl")
include("test_bGrad.jl")
include("test_exactdmd.jl")
include("test_projb.jl")
include("test_setup.jl")
include("test_solveDMD_withBFGS.jl")
include("test_solveDMD_withSVRG_huber.jl")
include("test_solveDMD_withSVRG.jl")


text:
import Base: ==, convert, promote_rule, length, reverse, show, div

struct Line{T <: Number}
    m::Matrix{T}
    function Line{T}(m::Matrix{T}) where {T <: Number}
        @assert size(m) == (2, 2) "Invalid values."
        new(m)
    end
    Line{T}(lx::T, ly::T, rx::T, ry::T) where {T <: Number} =
        new(Matrix([lx rx; ly ry]))
end

Line(m::Matrix{T}) where {T <: Number} = Line{T}(m)

start(l::Line) = l.m[:, 1]
endof(l::Line) = l.m[:, 2]

@inline sx(l::Line) = l.m[1, 1]
@inline sy(l::Line) = l.m[2, 1]
@inline ex(l::Line) = l.m[1, 2]
@inline ey(l::Line) = l.m[2, 2]

area(::Line{T}) where {T <:Number} = zero(T)

function Line(lx::Number, ly::Number, rx::Number, ry::Number)
    t = promote(lx, ly, rx, ry)
    return Line{typeof(t[1])}(t...)
end

convert(::Type{Line{T}}, r::Line{S}) where {T <: Number, S <: Number} =
    Line{T}(Matrix{T}(r.m))

promote_rule(::Type{Line{T}}, ::Type{Line{S}}) where {T <: Number, S <: Number} =
    Line{promote_type(T, S)}

show(io::IO, r::Line) =
    print(io, "Line:[$(r.m[1, 1]) $(r.m[2, 1]) $(r.m[1, 2]) $(r.m[2, 2])]")

==(l1::Line{T}, l2::Line{T}) where {T <: Number} = all(iszero.(l1.m - l2.m))
==(l1::Line, l2::Line) = ==(promote(l1, l2)...)

function reverse(l::Line)
    m = copy(l.m)
    m[:, 1], m[:, 2] = m[:, 2], m[:, 1]
    return Line(m)
end

axis_parallel(l::Line{T}; dir::Int=1) where {T <: Number} =
    iszero(l.m[dir, 1] - l.m[dir, 2])

"""
```
    isHorizontal(l::Line) -> Bool
    isVertcal(l::Line) -> Bool
```
If the `Line` is horizontal or vertical.
"""
isHorizontal(l::Line) = axis_parallel(l, dir=2)
isVertical(l::Line)   = axis_parallel(l, dir=1)

"""
```
    length(l::Line) -> Float64
```
The length of the line segment.
"""
length(l::Line) = (v = l.m[:, 1] - l.m[:, 2]; sqrt(dot(v, v)))

"""
```
    ratio(l1::Line{T}, p::Vector{T}) -> r::Real
```
If `p` is on `l1` it divides the line at ratio `r:(1-r)` else nothing.
"""
function ratio(l::Line{T}, p::Vector{T}) where {T <: Real}
    dv = l.m[:, 2] - l.m[:, 1]
    dp = p - l.m[:, 1]
    r, c = !iszero(dv[1]) ? (dp[1] / dv[1], 1) : (dp[2] / dv[2], 2)
    if c == 1
        tp = dv[2]*r + l.m[2, 1]
        iszero(tp - p[2]) && return r
    else
        iszero(dp[1]) && return r
    end
    return nothing
end

ratio(l::Line{T}, p::Vector{S}) where {T <: Number, S <: Number} =
    (ST = promote_type(S, T);
     ratio(convert(Line{ST}, l), convert(Vector{ST}, p)))

div(l::Line{T}, r::R) where {T <: Number, R <: Real} =
    l.m[:, 1]*(one(R) - r) + l.m[:, 2]*r

"""
```
    intersects(l1::Line{T}, l2::Line{T}) where {T <: Real} -> Bool
```
If `l1` and `l2` intersect each other.
"""
function intersects(l1::Line{T}, l2::Line{T}) where {T <: Real}
    l = [l1 l2; l2 l1]
    l1l21 = parallelogram_area(hcat(l1.m, l2.m[:, 1]))
    l1l22 = parallelogram_area(hcat(l1.m, l2.m[:, 2]))
    l2l11 = parallelogram_area(hcat(l2.m, l1.m[:, 1]))
    l2l12 = parallelogram_area(hcat(l2.m, l1.m[:, 2]))
    t = [l1l21 l1l22; l2l11 l2l12]
    for i = 1:2
        for j = 1:2
            if iszero(t[i, j])
                r = ratio(l[i, 1], l[i, 2].m[:, j])
                r === nothing && continue
                return zero(T) <= notvoid(r) <= one(T)
            end
        end
    end
    return t[1, 1]*t[1, 2] < zero(T) && t[2, 1]*t[2, 2] < zero(T)
end

intersects(l1::Line, l2::Line) = intersects(promote(l1, l2)...)

"""
```
    merge_axis_aligned(alines::Vector{Line{T}},
                       axis::Int=1,
                       order::Symbol=:increasing,
                       tol::T=pcTol(T)) -> Vector{Line{T}}
```
Given an array of axis aligned lines, if the line ends touch the lines are
merged into a larger segment. Lines which are not touching the other lines are
left intact.

`order` parameter can be in `:increasing` or `:decreasing` order in the direction
of the axis.

`axis` parameter can be `1` for horizontal lines and `2` for vertical lines.
"""
function merge_axis_aligned(alines::Vector{Line{T}},
                            axis::Int=1,
                            order::Symbol=:increasing,
                            tol::T=pcTol(T)) where {T}
    length(alines) == 0 && return Line{T}[]
    pl = alines[1]
    m = copy(pl.m)
    oaxis = axis == 1 ? 2 : 1
    vl = Vector{Line{T}}()
    for i = 2:length(alines)
        l = alines[i]
        if iszero(l.m[oaxis, 1] - pl.m[oaxis, 1])
            if order === :increasing && abs(l.m[axis, 1] - pl.m[axis, 2]) <= tol
                m[axis, 2] = l.m[axis, 2]
            elseif order === :decreasing &&
                abs(l.m[axis, 2] - pl.m[axis, 1]) <= tol
                m[axis, 1] = l.m[axis, 1]
            else
                push!(vl, Line{T}(m))
                m = copy(l.m)
            end
        else
            push!(vl, Line{T}(m))
            m = copy(l.m)
        end
        pl = l
    end
    push!(vl, Line{T}(m))
    return vl
end

function intersect_axis_aligned(hl::Line{T},
                                vl::Line{T}, tol::T) where T <: Number
    x, y  = sx(vl), sy(hl)
    if sx(hl) > ex(hl)
        hl = reverse(hl)
    end
    if sy(vl) > ey(vl)
        vl = reverse(vl)
    end
    if sx(hl) - tol <= x <= ex(hl) + tol && sy(vl) - tol <= y <= ey(vl) + tol
        return T[x, y]
    else
        return T[]
    end
end

function intersect_axis_aligned(hl::Line{T1},
                                vl::Line{T2},
                                tol::T=pcTol(T)) where {T1 <: Number,
                                                        T2 <: Number,
                                                        T <: Number}
    ST = promote_type(T1, T2, T)
    return intersect_axis_aligned(convert(Line{ST}, hl),
                                  convert(Line{ST}, vl),
                                  convert(ST, tol))
end


text:
# TODO:
# - Test both exact and inexact samples
# - Read from HDF5, CSV, Parquet, Arrow files, CSV datasets
# - Read from both S3 and Internet
# - Test with and without S3FS
# - Test shuffled, similar files, default, invalidated location but reused sample or reused location
# - Verify no error, length of returned sample, location files and their nrows,
# (eventually) ensure that rows come from the right files
@testset "$exact_or_inexact sample collected from $file_extension $(single_file ? "file" : "directory") on $on $with_or_without_s3fs S3FS $with_or_without_shuffled shuffled reusing $reusing" for exact_or_inexact in
                                                                                [
        "Exact",
        "Inexact",
    ],
    (file_extension, single_file, on, src_nrows) in [
        ("h5", true, "S3", 10),
        ("h5", true, "Internet", 10),
        ("csv", true, "S3", 150),
        ("parquet", true, "S3", 150),
        ("arrow", true, "S3", 150),
        ("csv", true, "Internet", 150),
        ("parquet", true, "Internet", 150),
        ("arrow", true, "Internet", 150),
        ("csv", false, "S3", 150 * 10),
        ("parquet", false, "S3", 150 * 10),
        ("arrow", false, "S3", 150 * 10),
    ],
    with_or_without_shuffled in ["with", "without"],
    reusing in ["nothing", "sample", "location", "sample and location"],
    with_or_without_s3fs in ["with", "without"]

    # Use job with appropriate sample collection configuration
    use_job_for_testing(
        sample_rate = 2,
        max_exact_sample_length = exact_or_inexact == "Exact" ? 1_024_000 : 0,
        with_s3fs = with_or_without_s3fs == "with",
    ) do

        # Use data to collect a sample from
        src_name = use_data(file_extension, on, single_file)

        # Construct location
        if reusing != "nothing"
            Remote(src_name, source_invalid = true, sample_invalid = true)
        end
        remote_source = Remote(
            src_name,
            source_invalid = (reusing == "nothing" || reusing == "sample"),
            sample_invalid = (reusing == "nothing" || reusing == "location"),
            shuffled = with_or_without_shuffled == "with",
        )

        # Verify the location
        if contains(src_name, "h5")
            @test remote_source.ndims == 2
            @test !contains(remote_source.path, "DS1")
            @test remote_source.subpath == "DS1"
            @test remote_source.size[1] == src_nrows
        else
            @test remote_source.nbytes > 0
            @test remote_source.nrows == src_nrows

            if contains(src_name, "dir")
                @test length(remote_source.files) == 10
                for f in remote_source.files
                    @test f["nrows"] == 150
                end
            else
                @test length(remote_source.files) == 1
            end
        end

        # TODO: Add these tests
        # TODO: Fix sample collection in the optimizations/reuse and nbytes

        # Verify the sample
        sample_nrows =
            contains(src_name, "h5") ? size(remote_source.sample.value, 1) :
            nrow(remote_source.sample.value)
        if exact_or_inexact == "Exact"
            @test sample_nrows == src_nrows
        else
            @test sample_nrows == cld(src_nrows, 2)
        end
    end
end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
using Flux
using Zygote
using MAT
#using Plots
using Statistics
using Distributions
using LinearAlgebra
using Random
using BSON
using JLD
#using PyPlot
using CUDA
#using Base.Iterators: partition

#gpu_id = 1
#CUDA.device!(gpu_id)

NVol = 36 # number of volumes
idx_valid = 1

# Unet define
channels_out=1

Nch_init = 16
kernelsize = (3, 3)

mf = matopen("../../Data/Data_processed_36Vol_NoCrop.mat")
Xtrain = read(mf, "Xtrain");
Ytrain = read(mf, "Ytrain");
close(mf)

for idx_subj = 1:size(Xtrain,5)
    Xtrain[:,:,:,:,idx_subj] = Xtrain[:,:,:,:,idx_subj]./maximum(Xtrain[:,:,:,:,idx_subj]);
    Ytrain[:,:,:,:,idx_subj] = Ytrain[:,:,:,:,idx_subj]./maximum(Ytrain[:,:,:,:,idx_subj]);
end

XValid = Xtrain[:,:,:,:,idx_valid];
YValid = Ytrain[:,:,:,:,idx_valid];

Xtrain = nothing;
Ytrain = nothing;

mf = matopen("../../Data/T1wBrain_processed_NoCrop.mat")
T1wBrain_resize_ds = read(mf, "T1wBrain_resize_ds");
close(mf)

for idx_subj = 1:size(T1wBrain_resize_ds,4)
    T1wBrain_resize_ds[:,:,:,idx_subj] = T1wBrain_resize_ds[:,:,:,idx_subj]./maximum(T1wBrain_resize_ds[:,:,:,idx_subj]);
end

Sall = size(XValid);
XValid = reshape(XValid,Sall[1],Sall[2],1,Sall[3]*Sall[4]);
YValid = reshape(YValid,Sall[1],Sall[2],1,Sall[3]*Sall[4]);

XValid_new = zeros(Float32,size(XValid,1),size(XValid,1),2,size(XValid,4));
XValid_new[:,:,1,:] = XValid;
XValid_new[:,:,2,:] = reshape(repeat(T1wBrain_resize_ds[:,:,:,idx_valid],1,1,NVol),size(XValid,1),size(XValid,1),1,size(XValid,4));

XValid = XValid_new;
XValid_new=nothing;

#XValid = JLD.load("Data/XYtrain_XYvalid_processed_36Vol_T1w.jld","XValid")
#YValid = JLD.load("Data/XYtrain_XYvalid_processed_36Vol_T1w.jld","YValid")

Nslc = 173
Nvol = 36

Nim = size(XValid,1);
#XValid = reshape(XValid,size(XValid,1),size(XValid,2),size(XValid,3),Nslc,Nvol)
XValid = reshape(XValid[:,:,:,:],size(XValid,1),size(XValid,2),2,Nslc,Nvol);
YValid = reshape(YValid,size(YValid,1),size(YValid,2),Nslc,Nvol);

#TestModelFile = "model_processed_36Vol_T1W_L2_SmallUNet30_k5_s14.bson"
#TestModelFile = "model_processed_36Vol_T1W_L2_SmallUNetGlobalSkip30_k5_s14.bson"
#TestModelFile = "model_EncoderDecoder_2345_1_noT1W.bson"
#TestModelFile = "model_EncoderDecoder_2345_1_T1W.bson"
#TestModelFile = "model_EncoderDecoderGlobalSkip_2345_1_T1W.bson"
#TestModelFile = "model_EncoderDecoderConvDownGlobalSkip_2345_1_T1WEpoch13.bson"
#TestModelFile = "model_EncoderDecoderConvDownReslearn_2345_1_T1WEpoch9.bson"
TestModelFile = "model_EncoderDecoderConvDownGlobalSkipMAE_2345_1_T1WEpoch24.bson"
BSON.@load string("Results/",TestModelFile) model
channels_in=2

XtestOutput4write = zeros(size(YValid));

Xtest = zeros(Float32,size(XValid,1),size(XValid,2),channels_in,1);
#Ytest = zeros(Float32,size(XValid,1),size(XValid,2),1,1)
for idx_dir_test in 1:size(XValid,5)
  for idx_slc_test in 1:size(XValid,4)
    Xtest[:,:,:,1] = XValid[:,:,1:channels_in,idx_slc_test,idx_dir_test];
    #Ytest[:,:,1,1] = YValid[:,:,idx_slc_test,idx_dir_test]
#    XtestOutput4write[:,:,idx_slc_test,idx_dir_test] = model(Xtest);
    XtestOutput4write[:,:,idx_slc_test,idx_dir_test] = model(Xtest).+Xtest[:,:,1]
  end
end

file = matopen(string("FromMSI/Results/CV/",TestModelFile[1:end-5],".mat"), "w")
write(file, "XtestOutput4write", XtestOutput4write)
close(file)


text:
function print_header(o)
    @printf "(%s_w,%s)-aCMA-ES (mu_w=%2.1f,w_1=%2d%%) in dimension %d (seed=%s, %s)\n" o.p.weights.μ o.p.λ o.p.weights.μeff round(Int, 100*o.p.weights.weights[1]) o.p.n o.p.seed now()
end
function print_state(o)
    l = o.logger
    push!(l.times, time())
    if o.stop.it == 1
        print_header(o)
        @printf "%6s %8s   %14s  %9s  %10s %9s\n" "iter" "fevals" "function value" "sigma" "axis ratio" "time[s]"
    end
    @printf "%6.d %8.d   %.8e   %.2e  %10.3e %9.3f\n" o.stop.it o.stop.it * o.p.λ + noisefevals(o.p.noise_handling) l.fmedian[end] sigma(o.p) maximum(o.p.cov.C.sqrtvalues)/minimum(o.p.cov.C.sqrtvalues) l.times[end] - o.stop.t0
end
function print_result(o)
    if o.stop.reason == :none
        return
    end
    println("  termination reason: $(o.stop.reason) = $(getproperty(o.stop, o.stop.reason)) ($(now()))")
    println("  lowest observed function value: $(fbest(o)) at $(xbest(o))")
    println("  population mean: $(population_mean(o))")
end

max_stars_repo_name: jbrea/CMAEvolutionStrategy

max_issues_repo_name: jbrea/CMAEvolutionStrategy


text:
using JuLIP.Chemistry: atomic_number
using JuLIP: Atoms
using StaticArrays: SVector

abstract type MSitePotential <: AbstractCalculator end

abstract type MPairPotential <: MSitePotential end

# Experimental Multi-component codes

# ----------------------------------------------------------------------
#    Managing a list of species
# ----------------------------------------------------------------------

abstract type AbstractZList end

"""
`ZList` and `SZList{NZ}` : simple data structures that store a list
of species and convert between atomic numbers and the index in the list.
Can be constructed via
* `ZList(zors)` : where `zors` is an Integer  or `Symbol` (single species)
* `ZList(zs1, zs2, ..., zsn)`
* `ZList([sz1, zs2, ..., zsn])`
* All of these take a kwarg `static = {true, false}`; if `true`, then `ZList`
will return a `SZList{NZ}` for (possibly) faster access.
"""
struct ZList <: AbstractZList
   list::Vector{Int16}
end

Base.length(zlist::AbstractZList) = length(zlist.list)

struct SZList{N} <: AbstractZList
   list::SVector{N, Int16}
end

ZList(zlist::AbstractVector{<: Integer}; static = false) = (
   static ? SZList(SVector( (Int16.(sort(zlist)))... ))
          :  ZList( convert(Vector{Int16}, sort(zlist)) ))

ZList(s::Symbol; kwargs...) =
      ZList( [ atomic_number(s) ]; kwargs... )

ZList(S::AbstractVector{Symbol}; kwargs...) =
      ZList( atomic_number.(S); kwargs... )

ZList(args...; kwargs) =
      ZList( [args...]; kwargs...)

i2z(Zs::AbstractZList, i::Integer) = Zs.list[i]

function z2i(Zs::AbstractZList, z::Integer)
   for j = 1:length(Zs.list)
      if Zs.list[j] == z
         return j
      end
   end
   error("z = $z not found in ZList $(Zs.list)")
end

i2z(V, i::Integer) = i2z(V.zlist, i)
z2i(V, z::Integer) = z2i(V.zlist, z)

Dict(zlist::ZList) = Dict("__id__" => "JuLIP_ZList",
                          "list" => zlist.list)
Base.convert(::Val{:JuLIP_ZList}, D::Dict) = ZList(D)
ZList(D::Dict) = ZList(D["list"])

Dict(zlist::SZList) = Dict("__id__" => "JuLIP_SZList",
                          "list" => zlist.list)
Base.convert(::Val{:JuLIP_SZList}, D::Dict) = SZList(D)
SZList(D::Dict) = SZList(SVector(Int16.(D["list"])...))

# ----------------------------------------------------------------------
#       Basic calculator extensions for multi-species
# ----------------------------------------------------------------------

"""
`neigsz!(tmp, nlist::PairList, at::Atoms, i::Integer) -> j, R Z`

requires a temporary storage array `tmp` with fields
`tmp.R, tmp.Z`.
"""
function neigsz!(tmp, nlist::PairList, at::Atoms, i::Integer)
   j, R = neigs!(tmp.R, nlist, i)
   Z = tmp.Z
   for n = 1:length(j)
      Z[n] = at.Z[j[n]]
   end
   return j, R, (@view Z[1:length(j)])
end

function neigsz(nlist::PairList, at::Atoms, i::Integer)
   j, R = NeighbourLists.neigs(nlist, i)
   return j, R, at.Z[j]
end

alloc_temp(V::MSitePotential, at::AbstractAtoms) =
      alloc_temp(V, maxneigs(neighbourlist(at, cutoff(V))))

alloc_temp(V::MSitePotential, N::Integer) =
      ( R = zeros(JVecF, N),
        Z = zeros(Int16, N), )

alloc_temp_d(V::MSitePotential, at::AbstractAtoms) =
      alloc_temp_d(V, maxneigs(neighbourlist(at, cutoff(V))))

alloc_temp_d(V::MSitePotential, N::Integer) =
      (dV = zeros(JVecF, N),
        R = zeros(JVecF, N),
        Z = zeros(Int16, N), )

energy(V::MSitePotential, at::AbstractAtoms; kwargs...) =
      energy!(alloc_temp(V, at), V, at; kwargs...)

virial(V::MSitePotential, at::AbstractAtoms; kwargs...) =
      virial!(alloc_temp_d(V, at), V, at; kwargs...)

forces(V::MSitePotential, at::AbstractAtoms{T}; kwargs...) where {T} =
      forces!(zeros(JVec{T}, length(at)), alloc_temp_d(V, at), V, at; kwargs...)

function energy!(tmp, calc::MSitePotential, at::Atoms{T};
                 domain=1:length(at)) where {T}
   E = zero(T)
   nlist = neighbourlist(at, cutoff(calc))
   for i in domain
      j, R, Z = neigsz!(tmp, nlist, at, i)
      E += evaluate!(tmp, calc, R, Z, Int16(at.Z[i]))
   end
   return E
end

function forces!(frc, tmp, calc::MSitePotential, at::Atoms{T};
                 domain=1:length(at), reset=true) where {T}
   if reset; fill!(frc, zero(JVec{T})); end
   nlist = neighbourlist(at, cutoff(calc))
   for i in domain
      j, R, Z = neigsz!(tmp, nlist, at, i)
      if length(j) > 0
         evaluate_d!(tmp.dV, tmp, calc, R, Z, Int16(at.Z[i]))
         for a = 1:length(j)
            frc[j[a]] -= tmp.dV[a]
            frc[i]    += tmp.dV[a]
         end
      end
   end
   return frc
end

function virial!(tmp, calc::MSitePotential, at::Atoms{T};
                 domain=1:length(at)) where {T}
   nlist = neighbourlist(at, cutoff(calc))
   vir = zero(JMat{T})
   for i in domain
      j, R, Z = neigsz!(tmp, nlist, at, i)
      if length(j) > 0
         evaluate_d!(tmp.dV, tmp, calc, R, Z, Int16(at.Z[i]))
         vir += site_virial(tmp.dV, R)
      end
   end
   return vir
end

evaluate(V::MSitePotential, R, Z, z) =
      evaluate!(alloc_temp(V, length(R)), V, R, Z, z)

evaluate_d(V::MSitePotential, R::AbstractVector{JVec{T}}, Z, z) where {T} =
      evaluate_d!(zeros(JVec{T}, length(R)),
                  alloc_temp_d(V, length(R)),
                  V, R, Z, z)

site_energy(V::MSitePotential, at::AbstractAtoms, i0::Integer) =
      energy(V, at; domain = (i0,))

site_energy_d(V::MSitePotential, at::AbstractAtoms, i0::Integer) =
      rmul!(forces(V, at; domain = (i0,)), -one(eltype(at)))

# -------- Dispatch for MPairPotentials ------------------------------------

evaluate!(tmp, V::MPairPotential,
          R::AbstractVector{<: JVec{T}}, Z::AbstractVector, z0) where {T} =
   length(R) == 0 ? zero(T) : T(0.5) * sum( evaluate!(tmp, V, norm(R[i]), Z[i], z0)
                                            for i = 1:length(R) )

function evaluate_d!(dV, tmp, V::MPairPotential, R, Z, z0)
   dV = tmp.dV
   for i = 1:length(R)
      r = norm(R[i])
      dV[i] = ((evaluate_d!(tmp, V, r, Z[i], z0)/r) * R[i]) / 2
   end
   return dV
end

evaluate(V::MPairPotential, r, z1, z0) =
      evaluate!(alloc_temp(V, 1), V, r, z1, z0)

evaluate_d(V::MPairPotential, r, z1, z0) =
      evaluate_d!(alloc_temp_d(V, 1), V, r, z1, z0)

# ------------------ Preconditioners for MPairPotentials ----------------------

function precon!(hEs, tmp, V::MPairPotential, R::AbstractVector{<: JVec},
                 Z::AbstractVector, z0, innerstab=T(0.0))
   n = length(R)
   for i = 1:n
      hEs[i,i] = precon!(tmp, V, norm(R[i]), R[i], Z[i], z0, innerstab)
   end
   return hEs
end

# # an FF preconditioner for pair potentials
# function precon!(tmp, V::PairPotential, r::T, R::JVec{T}, innerstab=T(0.1)
#                  ) where {T <: Number}
#    r = norm(R)
#    dV = evaluate_d!(tmp, V, r)
#    ddV = evaluate_dd!(tmp, V, r)
#    R̂ = R/r
#    return (1-innerstab) * (abs(ddV) * R̂ * R̂' + abs(dV / r) * (I - R̂ * R̂')) +
#              innerstab  * (abs(ddV) + abs(dV / r)) * I
# end

# struct MPairPotential <: MSitePotential
#    Z2V::Dict{Tuple{Int16,Int16}, WrappedPairPotential}
# end
#
# MPairPotential() =
#       MPairPotential(Dict{(Int16, Int16), WrappedPairPotential}())
#
# function MPairPotential(args...)
#    Vm = MPairPotential()
#    for p in args
#       @assert p isa Pair
#       z, V = _convert_multi_pp(sym::Symbol, p)
#       Vm.Z2V[z] = V
#    end
#    return Vm
# end
#
# F64fun(V::F64fun) = V
# _convert_multi_pp(sym::Tuple{Symbol, Symbol}, V) =
#       _convert_multi_pp(atomic_number.(sym), V)
# _convert_multi_pp(z::Tuple{<:Integer,<:Integer}, V) =
#       Int16.(z), V
#
# cutoff(V::MPairPotential) =
#    maximum( cutoff(V) for (key, V) in V.Z2V )


text:
@testset "201.bitwise-and-of-numbers-range.jl" begin
    @test range_bitwise_and(5, 7) == 4
    @test range_bitwise_and(0, 1) == 0
end


text:
"""
    dimlists = choose_dimensions(box::Box)

Choose combinations of dimensions that have been least recently used
for splitting in the vicinity of `box`. Passing these to
[`addpoint!`](@ref) will enable fresh estimates of the
corresponding coefficients of the polynomial.
"""
function choose_dimensions(box::Box{2})
    pairs = Union{Tuple{Int,Int},Tuple{Int}}[]
    n = ndims(box)
    if (isroot(box) && isleaf(box)) || n <= 2
        for i = 1:2:n
            if i < n
                push!(pairs, (i, i+1))
            else
                push!(pairs, (i,))
            end
        end
        return pairs
    end
    # Determine the least recently seen splits in the vicinity of
    # `box`.  We iterate over the splits until at least one dimension
    # has only one unseen pairing left. Allowing only unseen pairings,
    # we then construct a maximal matching
    # (https://en.wikipedia.org/wiki/Matching_(graph_theory))
    neven = n + isodd(n)
    agesentinel = typemax(Int)            # this value signals "never seen"
    age = fill(agesentinel, neven, neven) # recency of edge, younger is "worse"
    for i = 1:neven
        age[i,i] = 0   # a dimension can't be paired with itself
    end
    nfilled = fill(1, neven)
    agecounter = 1
    for split in splits(box)
        d1, d2 = split.dims
        if age[d1, d2] == agesentinel
            age[d1, d2] = age[d2, d1] = agecounter
            nfilled[d1] += 1
            nfilled[d2] += 1
        end
        agecounter += 1
        allfilled = true
        @inbounds @simd for i = 1:neven
            allfilled &= nfilled[i] >= neven-1
        end
        allfilled && break
    end
    # Build the weighted graph and solve the matching problem
    pairgraph = BlossomV.Matching(neven)
    for j = 1:neven, i = 1:j-1
        if age[i,j] > 0
            cost = agecounter - min(age[i,j], agecounter) # older has lower cost
            BlossomV.add_edge(pairgraph, i-1, j-1, cost)
        end
    end
    BlossomV.solve(pairgraph)
    seen = falses(n)
    extra = 0  # if a real node is paired with a fictitious node, add it last
    for i = 1:n
        seen[i] && continue
        j = BlossomV.get_match(pairgraph, i-1) + 1
        if j <= n
            if !seen[j]
                push!(pairs, (i, j))
            end
            seen[i] = seen[j] = true
        else
            extra = i
        end
    end
    if extra != 0
        push!(pairs, (extra,))
    end
    return pairs
end

"""
    c, g, Q, b = fit_quadratic(box)
    c, g, Q, b = fit_quadratic(f, box)

Return the coefficients of a quadratic model "centered" at `box`. The
model returns values

    m(x) = c + g'*Δx + (Δx'*Q*Δx)/2

where `Δx = x - b`.
"""
function fit_quadratic(f::Function, box::Box{2}; skip::Int=0)
    c, g, Q, b, y, prec, firstmatch = fit_quadratic_native(f, box; skip=skip)
    return canonicalize(c, g, Q, b, y, prec)
end
fit_quadratic(box::Box; skip::Int=0) = fit_quadratic(value, box; skip=skip)

"""
    c, g, Q, b, y, prec, firstmatch = fit_quadratic_native(box)
    c, g, Q, b, y, prec, firstmatch = fit_quadratic_native(f, box)

Return the coefficients of a CS2-native quadratic model "centered" at
`box`. The model returns values

    m(x) = c + g'*(x-b) + (1/2) ∑_{ij} (x[i] - b[i]) * Q[i,j] * (x[j] - χ[j,i])

where

             y_j          if i==j;
    χ[j,i] = b_j          if dimension i was split before or with j;
             2*y_j - b_j  if dimension i was split after j.

Precedence of splitting is determined with respect to `box`'s chain
(see [`CoordinateSplittingPTrees.chaintop`](@ref)).

See also [`CoordinateSplittingPTrees.fit_quadratic`](@ref).
"""
function fit_quadratic_native(f::Function, box::Box{2}; skip::Int=0)
    c, g, b, y, prec, firstmatch = fit_poly1(f, box; skip=skip)
    Q = coefficients_p(f, box; skip=skip)
    fit_poly2_diag!(f, Q, box, c, g, b, y, prec; skip=skip)
    return c, g, Q, b, y, prec, firstmatch
end
fit_quadratic_native(box::Box{2}; skip::Int=0) = fit_quadratic_native(value, box; skip=skip)

"""
    c, gc, Q, b = canonicalize(c, g, Q, b, y, prec)

Return the CS2-native quadratic model in canonical form,

    m(x) = c + gc'*Δx + (Δx'*Q*Δx)/2

where `Δx = x - b`.
"""
function canonicalize(c, g, Q, b, y, prec)
    gc = copy(g)
    n = length(g)
    for j = 1:n, i = 1:n
        χ = chi(j, i, b, y, prec)
        gc[i] += (b[j] - χ)*Q[i,j]/2
    end
    return c, gc, Q, b
end

function modelvalue(c, g, Q, b, y, prec::AbstractMatrix{Bool}, x)
    # For the CS2-native model
    n = length(x)
    val = c + g'*(x - b)
    for j = 1:n, i = j:n
        coef = CoordinateSplittingPTrees.Qcoef_value(i, j, x, b, y, prec)
        if coef != 0  # avoids using NaNs unless we need them
            val += Q[i,j]*coef*(1 + (i!=j))  # off-diagonals are done once
        end
    end
    return val
end

function modelvalue(c, g, Q, b, x)
    # For the canonical model
    Δx = x - b
    return c + g'*Δx + (Δx'*Q*Δx)/2
end

# Fit the first-order coefficients and set base positions
function fit_poly1(f::Function, box::Box{2,T}; skip::Int=0) where T
    box = getleaf(box)
    n = ndims(box)
    # Start at the top of the chain of splits as generated by `addpoint!`
    box, success = chaintop(box)
    success || error("no chain found at $box")

    b = position(box)
    y = similar(b)  # opposite point
    g = zeros(T, n)
    wassplit = falses(n)
    firstmatch = zeros(Int, n)
    # Rather than computing χ directly, store y and prec where
    # prec[i,j] is true if i ⪯ j. This avoids any risk of roundoff
    # error in the cancelation of Q coefficients (see Qcoef_value below).
    prec = fill(false, n, n)

    # Visit splits to cover each dimension once.
    # This suffices to determine the constant and gradient.
    c = f(box)
    npairs = ceil(Int, n/2)
    for paircounter = 1:npairs
        split = box.split
        d1, d2 = split.dims
        @assert(d1 != d2)
        @assert(!wassplit[d1] && (d2 > n || !wassplit[d2]))
        x1, x2 = split.xs
        f1, f2, f12 = f.(split.others.children)
        f0 = f(split.self)
        y[d1], wassplit[d1], firstmatch[d1] = x1, true, d2
        if d2 <= n
            y[d2], wassplit[d2], firstmatch[d2] = x2, true, d1
            prec[:, d2] = wassplit
        end
        prec[:, d1] = wassplit
        g[d1] = (f1 - f0)/(x1 - b[d1])
        if d2 <= n
            g[d2] = (f2 - f0)/(x2 - b[d2])
        end

        box = split.others.children[3]
    end
    return c, g, b, y, prec, firstmatch
end
fit_poly1(box) = fit_poly1(value, box)

function fit_poly2_diag!(f::Function, Q, box::Box{2}, c, g, b, y, prec; skip::Int=0)
    newinfo(x, b, y) = (x != b) & (x != y)
    function diag_equation(d, xnew, g, Q, x, b, y, prec)
        n = length(x)
        xd = x[d]
        coef = Qcoef_lineq(d, d, xnew-xd, xd, b[d], y[d], prec[d,d])
        consts = -g
        for k = 1:n
            k == d && continue  # already in coef
            xc = Qcoef_lineq(d, k, xnew-xd, x[k], b[k], y[k], prec[d,k])
            consts -= Q[d,k] * xc
        end
        return coef, consts
    end
    function process_split!(f, Q, checked, split, dimidx, x, g, b, y, prec)
        n = length(b)
        d = split.dims[dimidx]
        (d > n || checked[d]) && return 0
        bd = b[d]
        yd = y[d]
        xd = split.xs[dimidx]
        if newinfo(xd, bd, yd)
            coef, consts = diag_equation(d, xd, g[d], Q, x, b, y, prec)
            child = split.others.children[dimidx]
            ref = split.self
            Q[d,d] = (consts + (f(child)-f(ref))/(xd-x[d]))/coef
            checked[d] = true
            return 1
        end
        # Try it the other way too (we don't know which child we came from)
        if newinfo(x[d], bd, yd)
            xd, x[d] = x[d], xd
            d2 = split.dims[3-dimidx]
            xd2 = split.xs[3-dimidx]
            if d2 <= n
                xd2, x[d2] = x[d2], xd2
            end
            coef, consts = diag_equation(d, xd, g[d], Q, x, b, y, prec)
            child = split.others.children[3-dimidx]
            ref = split.others.children[3]
            Q[d,d] = (consts + (f(child)-f(ref))/(xd-x[d]))/coef
            checked[d] = true
            xd, x[d] = x[d], xd  # restore
            if d2 <= n
                xd2, x[d2] = x[d2], xd2
            end
            return 1
        end
        return 0
    end

    n = ndims(box)
    x = position(box)
    filled = similar(x, Bool)
    checked = falses(n)
    nremaining = n
    nsplits = 0
    for split in splits(box)
        nremaining <= 0 && break
        nsplits += 1
        nsplits <= skip && continue
        position!(x, filled, split.self)
        nremaining -= process_split!(f, Q, checked, split, 1, x, g, b, y, prec)
        nremaining -= process_split!(f, Q, checked, split, 2, x, g, b, y, prec)
    end
    return Q
end
fit_poly2_diag!(Q, box::Box{2}, c, g, b, y, prec) = fit_poly2_diag!(value, Q, box, c, g, b, y, prec)

function chi(j, i, b, y, prec)
    i == j && return y[j]
    return prec[i,j] ? b[j] : 2*y[j]-b[j]
end

function chi(b, y, prec)
    n = length(b)
    return [chi(j, i, b, y, prec) for j=1:n, i=1:n]
end

"""
    coef = Qcoef_value(i, j, x, b, y, prec)

Compute the symmetrized coefficient of `Q[i,j]` in the model value
expression.  The native expression is shown in the help for
[`CoordinateSplittingPTrees.fit_quadratic_native`](@ref). Note that
this expression is asymmetric. This function returns a symmetrized
version,

    ((x[i] - b[i])*(x[j] - χ[j,i]) + (x[j] - b[j])*(x[i] - χ[i,j]))/4

which is the mean of the coefficients for `Q[i,j]` and `Q[j,i]`. In
particular, it returns exactly 0 if the two should (analytically,
i.e. to infinite precision) cancel one another.
"""
function Qcoef_value(i, j, x, b, y, prec)
    i == j && return (x[i]-b[i])*(x[i]-y[i])/2
    bi, bj = b[i], b[j]
    yi, yj = y[i], y[j]
    xi, xj = x[i], x[j]
    pij, pji = prec[i, j], prec[j, i]
    (((pij & !pji) & (xi == yi)) | ((pji & !pij) & (xj == yj))) && return zero(xi*xj/2)
    χij = ifelse(pji, bi, 2yi-bi)
    χji = ifelse(pij, bj, 2yj-bj)
    return ((xi - bi)*(xj - χji) + (xj - bj)*(xi - χij))/4
end

function Qcoef_lineq(i, k, Δxi, xk::Real, bk::Real, yk::Real, precik::Bool)
    coef = zero(xk)
    if i == k
        (xk == bk && Δxi == yk-bk) || (xk == yk && Δxi == bk-yk) && return coef
        return Δxi/2 + xk - (bk+yk)/2
    end
    if precik
        xk == bk && return coef
        return coef + (xk - bk)
    end
    return coef + (xk - yk)
end

# An alternative to fit_poly1 and fit_poly2_diag! that accumulates 2x2
# equations. This does not depend on having chains, but it does mean
# that the gradient is dependent on all the off-diagonal elements of
# Q. (If some of those are missing, then g will be contaminated.)
function fit_quadratic_gdiag!(f::Function, Q, box::Box{2,T}, b=position(box); skip::Int=0) where T
    n = ndims(box)
    g = fill(T(NaN), n)
    qcoefs = fill(T(NaN), 2, n)
    rhs = fill(T(NaN), 2, n)
    filled = similar(g, Bool)
    x = position(box)
    nremaining = 2*n
    c = f(box)
    # Iterate until we fill all coefficients or exhaust the tree
    nsplits = 0
    for split in splits(box)
        nremaining == 0 && break
        nsplits += 1
        nsplits <= skip && continue
        for (i, d) in enumerate(split.dims)
            d <= n || continue
            if isnan(qcoefs[2,d])
                # There's at least one to fill
                position!(x, filled, split.self)
                xs = split.xs[i]
                xcoef = (x[d]+xs)/2 - b[d]
                xcoef == qcoefs[1,d] && continue  # not independent of the first
                idx = isnan(qcoefs[1,d]) ? 1 : 2
                Δx = xs - x[d]
                fd = (f(split.others.children[i]) - f(split.self))/Δx
                for j = 1:n
                    j == d && continue
                    fd -= Q[d,j]*(x[j] - b[j])
                end
                if isfinite(xcoef) && isfinite(fd)
                    qcoefs[idx,d] = xcoef
                    rhs[idx,d] = fd
                    nremaining -= 1
                end
            end
        end
    end
    # Solve the 2x2 equations
    for d = 1:n
        if !isnan(qcoefs[2,d])
            a12, a22 = qcoefs[1,d], qcoefs[2,d]  # a11 and a21 are both 1
            r1, r2 = rhs[1,d], rhs[2,d]
            det = a22 - a12    # manual inverse of 2x2 matrix
            g[d] = (a22*r1 - a12*r2)/det
            Q[d,d] = (-r1 + r2)/det
        end
    end
    return c, g, Q, b
end

fit_quadratic_gdiag!(Q, box::Box{2}, b=position(box); skip::Int=0) =
    fit_quadratic_gdiag!(value, Q, box, b; skip=skip)

### Lower-bound model
# Fit a model for which the observed function values are a lower bound for the model values

function lowerbound_model(f::Function, box::Box, optimizer, m::Integer=2*ndims(box))
    n = ndims(box)
    c = f(box)
    # Keep track of the points incorporated into the fit
    pq = PriorityQueue{typeof(box),typeof(c)}()
    function maybe_enqueue(box, val)
        leaf = getleaf(box)
        if !haskey(pq, leaf)
            enqueue!(pq, leaf, val)
        end
        nothing
    end
    # Allocate storage for the model coefficients and anything needed by `updater`
    Q = allocate_coefficients_p(n, box)
    g = similar(Q, n)
    x = position(box)
    xbx = similar(x)
    flag = similar(x, Bool)
    X = similar(x, m, n)
    XQX = similar(X)
    X2 = similar(X)
    Δf = Vector{typeof(c)}(undef, m)
    # Testing for independence
    ige = IGE{eltype(x)}(2n)
    xtmp = similar(x, (2n,))
    function isindependent(box)
        position!(x, flag, box)
        for i = 1:n
            xtmp[i] = x[i] - xbx[i]
            xtmp[i+n] = xtmp[i]^2
        end
        i = insert!(ige, xtmp, value(box))
        return i <= 2n && ige.rowtmp[i] != 0
    end
    function insertrow!(X, X2, Δf, i, box)
        for j = 1:n
            dx = x[j] - xbx[j]
            X[i,j] = dx
            X2[i,j] = dx^2/2
        end
        Δf[i] = value(box)
    end
    # Perform all box-specific setup (a box `bx` will be passed to this before solving the model)
    # The end results are `X`, which stores displacements from the base point in its columns,
    # and Δf, which stores the residual value that should ideally be fit by the gradient and the
    # diagonal of the Hessian.
    function updater(bx)
        empty!(pq)
        empty!(ige)
        c = f(bx)
        position!(xbx, flag, bx)
        # In case we run out of boxes, ensure these are zeroed
        fill!(X, 0)
        fill!(X2, 0)
        fill!(Δf, c)

        iter = splits(bx)
        state = skipsplits(iter, 0)
        coefficients_p!(f, Q, bx, iter, state, maybe_enqueue)
        # Unload the queue and test for independence
        irow = 0
        while !isempty(pq) && !iscomplete(ige)
            box = dequeue!(pq)
            if isindependent(box)
                insertrow!(X, X2, Δf, irow+=1, box)
            end
        end
        while !isempty(pq) && irow < m
            box = dequeue!(pq)
            position!(x, flag, box)
            insertrow!(X, X2, Δf, irow+=1, box)
        end
        # Continue traversing the tree until we get enough boxes
        while irow + length(pq) < m
            ret = iterate(iter, state)
            ret === nothing && break
            split, state = ret
            for ichild = 0:maxchildren(box)-1
                child = getchild(split, ichild)
                maybe_enqueue(child, f(child))
            end
        end
        while !isempty(pq) && irow < m
            box = dequeue!(pq)
            position!(x, flag, box)
            insertrow!(X, X2, Δf, irow+=1, box)
        end
        # Zero the NaNs so they don't poison Δf
        nanzero!(Q)
        # From Δf subtract off everything except the gradient term and the diagonal of Q
        mul!(XQX, X, Q)
        XQX .*= X
        Δf .-= c .+ vec(sum(XQX; dims=2)) ./ 2
        true  # return upon success
    end
    # Set up the solver
    model = Model(optimizer)
    g = [Variable(model) for _ = 1:n]
    d = [Variable(model) for _ = 1:n]
    Xparam = Parameter(model; val=X)
    X2param = Parameter(model; val=X2)
    Δfp = Parameter(model; val=Δf)
    pred = @expression Xparam*g + X2param*d
    @constraint(model, pred >= Δfp)
    @objective(model, Minimize, sum(pred))
    return updater, model, g, d, Q
end
lowerbound_model(box::Box, optimizer, m::Integer=2*ndims(box)) = lowerbound_model(value, box, optimizer, m)

function fit_quadratic_lowerbound!(Q, updater, model, g, d, box::Box{2}, b=position(box))
    updater(box)
    solve!(model)
    gv, dv = Parametron.value.(model, g), Parametron.value.(model, d)
    for i in axes(Q, 1)
        Q[i,i] = dv[i]
    end
    return value(box), gv, Q, b
end

### Incremental Gaussian Elimination utilities

function Base.insert!(ige::IGE{T}, x::AbstractVector, newrhs; canswap::Bool=true) where T
    @inline myapprox(x, y, rtol) = (x == y) | (abs(x-y) < rtol*(abs(x) + abs(y)))
    rtol = 1000*eps(T)

    coefs, rhs, newrow = ige.coefs, ige.rhs, ige.rowtmp
    copyto!(newrow, x)  # to avoid destroying x
    n = length(rhs)
    i = 0
    while i < n
        i += 1
        newrow[i]  == 0 && continue
        coefs[i,i] == 0 && break
        if canswap && abs(newrow[i]) > abs(coefs[i,i])
            # Swap (for numeric stability)
            for j = i:n
                newrow[j], coefs[i,j] = coefs[i,j], newrow[j]
            end
            newrhs, rhs[i] = rhs[i], newrhs
        end
        # Cancel the ith value
        c = newrow[i]/coefs[i,i]
        @simd for j = i+1:n
            sub = c * coefs[i,j]
            newrow[j] = ifelse(myapprox(newrow[j], sub, rtol), zero(sub), newrow[j] - sub)
        end
        newrow[i] = 0  # rather than subtracting, just set it (to avoid roundoff error)
        newrhs -= c*rhs[i]
    end
    if i <= n && newrow[i] != 0
        # Insert new row
        for j = i:n
            coefs[i,j] = newrow[j]
        end
        rhs[i] = newrhs
    end
    return i
end

function iscomplete(ige::IGE)
    coefs = ige.coefs
    n = size(coefs, 1)
    for i = 1:n
        coefs[i,i] == 0 && return false
    end
    return true
end

function solve(ige::IGE)
    return ige.coefs \ ige.rhs
end

function Base.empty!(ige::IGE)
    fill!(ige.coefs, 0)
    fill!(ige.rhs, 0)
    return ige
end


text:
x = 5
sqrt(x)


text:
# BNF from https://semver.org/

#=
<valid semver> ::= <version core>
                 | <version core> "-" <pre-release>
                 | <version core> "+" <build>
                 | <version core> "-" <pre-release> "+" <build>

<version core> ::= <major> "." <minor> "." <patch>

<major> ::= <numeric identifier>

<minor> ::= <numeric identifier>

<patch> ::= <numeric identifier>

<pre-release> ::= <dot-separated pre-release identifiers>

<dot-separated pre-release identifiers> ::= <pre-release identifier>
                                          | <pre-release identifier> "." <dot-separated pre-release identifiers>

<build> ::= <dot-separated build  identifiers>

<dot-separated build identifiers> ::= <build identifier>
                                    | <build identifier> "." <dot-separated build identifiers>

<pre-release identifier> ::= <alphanumeric identifier>
                           | <numeric identifier>

<build identifier> ::= <alphanumeric identifier>
                     | <digits>

<alphanumeric identifier> ::= <non-digit>
                            | <non-digit> <identifier characters>
                            | <identifier characters> <non-digit>
                            | <identifier characters> <non-digit> <identifier characters>

<numeric identifier> ::= "0"
                       | <positive digit>
                       | <positive digit> <digits>

<identifier characters> ::= <identifier character>
                          | <identifier character> <identifier characters>

<identifier character> ::= <digit>
                         | <non-digit>

<non-digit> ::= <letter>
              | "-"

<digits> ::= <digit>
           | <digit> <digits>

<digit> ::= "0"
          | <positive digit>

<positive digit> ::= "1" | "2" | "3" | "4" | "5" | "6" | "7" | "8" | "9"

<letter> ::= "A" | "B" | "C" | "D" | "E" | "F" | "G" | "H" | "I" | "J"
           | "K" | "L" | "M" | "N" | "O" | "P" | "Q" | "R" | "S" | "T"
           | "U" | "V" | "W" | "X" | "Y" | "Z" | "a" | "b" | "c" | "d"
           | "e" | "f" | "g" | "h" | "i" | "j" | "k" | "l" | "m" | "n"
           | "o" | "p" | "q" | "r" | "s" | "t" | "u" | "v" | "w" | "x"
           | "y" | "z"
=#

# To simplify how we construct the intermediate results, we have
# inserted some intermediate nodes to do the construction.  These are
# tagged with an asterisk after teh opening < sign in the name.  They
# are used to distinguish the entry points of recursive rules from the
# recursion steps.

SemVerGrammar = BNFGrammar(:SemVer)

#=
macro ref(name)
    Expr(:call, :BNFRef, :SemVerGrammar, name)
end
=#

function gref(name)
    BNFRef(SemVerGrammar, name)
end

str2int(s) = parse(Int, s)

function undot(x)
    result = []
    function ud(x)
        if x == '.'
            return
        end
        if x isa AbstractVector
            for x1 in x
                ud(x1)
            end
        else
            push!(result, x)
        end
    end
    ud(x)
    return Tuple(result)
end

DerivationRule(
    SemVerGrammar, "<valid semver>",
    Alternatives(Constructor(gref("<version core>"),
                               v -> VersionNumber(v...)),
                  Constructor(Sequence(gref("<version core>"),
                                         CharacterLiteral('-'),
                                         gref("<pre-release>")),
                               function(v)
                                   (vc, dash, pre) = v
                                   VersionNumber(vc..., pre)
                               end),
                  Constructor(Sequence(gref("<version core>"),
                                         CharacterLiteral('+'),
                                         gref("<build>")),
                               function(v)
                                   (vc, plus, build) = v
                                   VersionNumber(vc..., (), build)
                               end),
                  Constructor(Sequence(gref("<version core>"),
                                         CharacterLiteral('-'),
                                         gref("<pre-release>"),
                                         CharacterLiteral('+'),
                                         gref("<build>")),
                               function(v)
                                   (vc, dash, pre, plus, build) = v
                                   VersionNumber(vc..., pre, build)
                               end)))

DerivationRule(
    SemVerGrammar, "<version core>",
    Constructor(Sequence(gref("<major>"), CharacterLiteral('.'),
                           gref("<minor>"), CharacterLiteral('.'),
                           gref("<patch>")),
                 undot))

DerivationRule(
    SemVerGrammar, "<major>",
    gref("<numeric identifier>"))

DerivationRule(
    SemVerGrammar, "<minor>",
    gref("<numeric identifier>"))

DerivationRule(
    SemVerGrammar, "<patch>",
    gref("<numeric identifier>"))

DerivationRule(
    SemVerGrammar, "<pre-release>",
    Constructor(
        gref("<dot-separated pre-release identifiers>"),
        undot))

DerivationRule(
    SemVerGrammar, "<dot-separated pre-release identifiers>",
    Alternatives(
        gref("<pre-release identifier>"),
        Sequence(gref("<pre-release identifier>"),
                  CharacterLiteral('.'),
                  gref("<dot-separated pre-release identifiers>"))))

DerivationRule(
    SemVerGrammar, "<build>",
    gref("<dot-separated build identifiers>"))

DerivationRule(
    SemVerGrammar, "<dot-separated build identifiers>",
    Constructor(
        Alternatives(
            gref("<build identifier>"),
            Sequence(gref("<build identifier>"),
                      CharacterLiteral('.'),
                      gref("<dot-separated build identifiers>"))),
        undot))

DerivationRule(
    SemVerGrammar, "<pre-release identifier>",
    Alternatives(gref("<alphanumeric identifier>"),
                  gref("<numeric identifier>")))

DerivationRule(
    SemVerGrammar, "<build identifier>",
    Alternatives(gref("<alphanumeric identifier>"),
                  Constructor(gref("<digits>"),
                               str2int)))

DerivationRule(
    SemVerGrammar, "<alphanumeric identifier>",
    StringCollector(
        Alternatives(gref("<non-digit>"),
                      Sequence(gref("<non-digit>"),
                                gref("<identifier characters>")),
                      Sequence(gref("<identifier characters>"),
                                gref("<non-digit>")),
                      Sequence(gref("<identifier characters>"),
                                gref("<non-digit>"),
                                gref("<identifier characters>")))))

DerivationRule(
    SemVerGrammar, "<numeric identifier>",
    Constructor(
        StringCollector(
            Alternatives(CharacterLiteral('0'),
                          gref("<positive digit>"),
                          Sequence(gref("<positive digit>"),
                                    gref("<digits>")))),
        str2int))

DerivationRule(
    SemVerGrammar, "<identifier characters>",
    StringCollector(
        gref("<*identifier characters>")))

DerivationRule(
    SemVerGrammar, "<*identifier characters>",
    Alternatives(gref("<identifier character>"),
                  Sequence(gref("<identifier character>"),
                            gref("<*identifier characters>"))))

DerivationRule(
    SemVerGrammar, "<identifier character>",
    Alternatives(gref("<digit>"),
                  gref("<non-digit>")))

DerivationRule(
    SemVerGrammar, "<non-digit>",
    Alternatives(gref("<letter>"),
                  CharacterLiteral('-')))

DerivationRule(
    SemVerGrammar, "<digits>",
    StringCollector(
        gref("<*digits>")))

DerivationRule(
    SemVerGrammar, "<*digits>",
    Alternatives(
        # Reordered the alternatives to reduce calls to <digit>.
        Sequence(gref("<digit>"),
                  gref("<*digits>")),
        gref("<digit>")))

DerivationRule(
    SemVerGrammar, "<digit>",
    Alternatives(CharacterLiteral('0'),
                  gref("<positive digit>")))

DerivationRule(
    SemVerGrammar, "<positive digit>",
    Alternatives([CharacterLiteral(c) for c in '1':'9']...))

DerivationRule(
    SemVerGrammar, "<letter>",
    Alternatives([CharacterLiteral(c) for c in 'A':'Z']...,
                  [CharacterLiteral(c) for c in 'a':'z']...))

@testset "SemVer grammar" begin
    @test undot("foo") == ("foo",)
    @test undot(["foo", '.', ["bar"]]) == ("foo", "bar")
    @test undot(["foo", '.', ["bar", '.', ["baz"]]]) == ("foo", "bar", "baz")

    @test recognize(CharacterLiteral('a'), "abc") == (true, 'a', 2)
    @test recognize(CharacterLiteral('a'), "b") == (false, nothing, 1)
    @test recognize(CharacterLiteral('a'), "") == (false, nothing, 1)
    @test recognize(SemVerGrammar["<letter>"], "abc") == (true, 'a', 2)
    @test recognize(SemVerGrammar["<digits>"], "1234 ") == (true, "1234", 5)
    @test recognize(gref("<positive digit>"), "0") == (false, nothing, 1)
    @test recognize(gref("<positive digit>"), "1") == (true, '1', 2)
    @test recognize(SemVerGrammar["<numeric identifier>"], "1") == (true, 1, 2)
    @test recognize(SemVerGrammar["<major>"],"2") == (true, 2, 2)
    @test recognize(SemVerGrammar["<minor>"],"21") == (true, 21, 3)
    @test recognize(SemVerGrammar["<version core>"], "3.2.11") == (true, (3, 2, 11), 7)
    @test recognize(SemVerGrammar["<valid semver>"], "1.2.3") == (true, VersionNumber(1, 2,3), 6)
    @test recognize(SemVerGrammar["<alphanumeric identifier>"], "dev") == (true, "dev", 4)
    @test recognize(SemVerGrammar["<pre-release identifier>"], "dev") == (true, "dev", 4)
    @test recognize(SemVerGrammar["<pre-release identifier>"], "15") == (true, 15, 3)
    @test recognize(SemVerGrammar["<pre-release>"], "dev") == (true, ("dev",), 4)
    @test recognize(SemVerGrammar["<pre-release>"],
                    "foo.3.bar") == (true, ("foo", 3, "bar"), 10)
    @test recognize(SemVerGrammar["<pre-release>"],"3") == (true, (3,), 2)
    @test recognize(SemVerGrammar["<pre-release>"], "dev") == (true, ("dev",), 4)
    @test recognize(SemVerGrammar["<pre-release>"],"dev.naha.2") == (true, ("dev", "naha", 2), 11)
    @test recognize(SemVerGrammar["<pre-release>"],"dev.naha-2") == (true, ("dev", "naha-2"), 11)
    @test recognize(SemVerGrammar["<valid semver>"], "1.2.3-dev.3") ==
        (true, VersionNumber(1, 2, 3, ("dev", 3)), 12)
end


text:
using CUDArt: CudaArray
import Base.reshape
reshape(a::CudaArray, dims::Tuple) = CudaArray(a.ptr, dims, a.dev)
reshape(a::CudaArray, dims::Integer...) = reshape(a, dims)

import CUDArt.to_host
to_host{T <: AbstractFloat}(a::Array{T}) = copy(a)


text:
using Documenter, RGBDSegmentation

makedocs(
    format = :html,
    sitename = "RGBDSegmentatione",
    pages = [
        "index.md",
        "Algorithms" => [
       		"JCSA"     => "JCSA.md",
                "DepthWMM" => "DepthWMM.md",
                "Felzenszwalb" => "felzenszwalb.md"
        ]
    ],
    modules = [RGBDSegmentation]
)


text:
###################################################
## Create city graph from lonlat polygon
## Automatically download it from OSM overpass API
###################################################

"""
    Create the road network inside the defined box.
    query the osm database using the overpass API if no osm file is given
"""
function queryOsmBox{T<:AbstractFloat}(left::T, right::T, bottom::T, top::T; osmfile="")
    if isempty(osmfile)
        url = "http://overpass-api.de/api/map?bbox=$(left),$(bottom),$(right),$(top)"
        osmfile = "$(Pkg.dir("RoutingNetworks"))/.cache/data.osm"
        println("Downloading OSM file (may take time)...")
        download(url, osmfile)
    end
    println("Parsing OSM file...")
    n = osm2network(osmfile)
    println("Creating routing Network")
    n = removeNodes(n,singleNodes(n))
    n = subsetNetwork(n,inPolygon(n,[(left,bottom),(left,top),(right,top),(right,bottom)]))
    n = removeNodes(n,singleNodes(n))

    return n
end

"""
    Query a polygon
    right now, query the whole bounding-box and subset the polygon: not optimal
"""
function queryOsmPolygon{T<:AbstractFloat}(poly::Vector{Tuple{T,T}}; osmfile="")
    (minLon,maxLon,minLat,maxLat) = boundingBox(poly)
    n = queryOsmBox(minLon,maxLon,minLat,maxLat,osmfile=osmfile)
    n = subsetNetwork(n,inPolygon(n,poly))
    n = removeNodes(n,singleNodes(n))
end

max_stars_repo_name: andrew-zhm/RoutingNetworks

max_issues_repo_name: andrew-zhm/RoutingNetworks

max_forks_repo_name: andrew-zhm/RoutingNetworks


text:
const CURRENT_NAME = [""]

################################################################################
## Makes a plot with all relevant transforms to visualize the data
################################################################################

function plot_transformed_data(df::DataFrame;
                               svg_path::String = "checkTransforms.svg",
                               overwrite::Bool = false)

    if !overwrite && isfile(svg_path)
        throw(ArgumentError("file '$svg_path' already exists. Either specify " *
                    "'overwrite=true' or 'svg_path=<some other file name>'."))
    end

    nms = names(df)

    n, p = size(df)
    relevant_columns = findall(x -> eltype(x) <: Real, [df[!, i] for i = 1:p])

    p_rel = length(relevant_columns)
    plots = Matrix{Plot}(undef, p_rel, 3)

    p_df = DataFrame(trait = names(df)[relevant_columns],
            untransformed = fill(NaN, p_rel),
            log = fill(NaN, p_rel),
            logit = fill(NaN, p_rel))

    for row in 1:length(relevant_columns)
        i = relevant_columns[row]
        col = df[!, i]
        if !(eltype(col) <: Real)
            error("all columns (except the first) must have elements of type Real")
        end

        name = nms[i]
        CURRENT_NAME[1] = name

        col = col[findall(!isnan, col)]

        pval = pvalue(JarqueBeraTest(col))
        p_df.untransformed[row] = pval

        p_original = plot(x=col,
                        Geom.histogram,
                        Guide.title(plot_title(row, name, "untransformed", pval)))
        plots[row, 1] = p_original

        if can_transform_log(col)
            logged_col = log_transform(col)
            pval = pvalue(JarqueBeraTest(logged_col))
            p_df.log[row] = pval
            plots[row, 2] = plot(x = logged_col,
                                Geom.histogram,
                                Guide.title(plot_title(row, name, "log", pval)))
        else
            plots[row, 2] = blank_plot("NA")
        end
        if can_transform_logit(col)
            logit_col = logit_transform(col)
            pval = pvalue(JarqueBeraTest(logit_col))
            p_df.logit[row] = pval
            plots[row, 3] = plot(x = logit_col,
                                Geom.histogram,
                                Guide.title(plot_title(row, name, "logit", pval)))
        else
            plots[row, 3] = blank_plot("NA")
        end

    end

    plt = gridstack(plots)
    draw(SVG(svg_path, 10inch, p * 2inch), plt)
    return p_df
end

function plot_title(i::Int, name::String, transform::String, p::Float64)
    p = round(p, sigdigits=3)
    return "$i -> $name $transform (p = $p)"
end

function blank_plot(title::String)
    return plot(x=Float64[], Geom.histogram, Guide.title(title))
end

################################################################################
## General tools for transforming data
################################################################################

function can_transform_log(x::Vector{<:Real})
    neg_ind = findfirst(y -> y < 0.0, x)
    return isnothing(neg_ind)
end

function log_transform(x::Vector{<:Real})
    n = length(x)
    pos_inds = findall(y -> y > 0.0, x)
    log_x = fill(NaN, n)

    if length(pos_inds) != n
        min = minimum(@view x[pos_inds])
        half_min = 0.5 * min
        log_half_min = log(half_min)
        zero_inds = findall(isequal(0.0), x)
        log_x[zero_inds] .= log_half_min

        @warn "0's in log-transform of $(CURRENT_NAME[1])." *
                " Converting all 0's to $half_min (half of the minium " *
                "value $min) prior to log-transforming."
    end

    for i in pos_inds
        log_x[i] = log(x[i])
    end
    return log_x
end

function can_transform_logit(x::Vector{<:Real})
    bad_ind = findfirst(y -> y > 1.0 || y < 0.0, x)
    return isnothing(bad_ind)
end

function logit_transform(x::Vector{<:Real})
    n = length(x)
    good_inds = findall(y -> y > 0.0 && y < 1.0, x)

    logit_x = fill(NaN, n)

    if length(good_inds) != n
        good_x = @view x[good_inds]
        half_min = 0.5 * minimum(good_x)
        logit_half_min = logit(half_min)

        half_max = 0.5 * (1.0 + maximum(good_x))
        logit_half_max = logit(half_max)

        zero_inds = findall(isequal(0.0), x)
        logit_x[zero_inds] .= logit_half_min

        one_inds = findall(isequal(1.0), x)
        logit_x[one_inds] .= logit_half_max

        @warn "0's and/or 1's in logit-transform of $(CURRENT_NAME[1])." *
                " Converting all 0's to $half_min (half of the non-zero minium " *
                "value $min) and all 1's to $half_max (halfways between the " *
                "non-one maximum $half_max and 1.0) prior to logit-transforming."
    end

    for i in good_inds
        logit_x[i] = logit(x[i])
    end
    return logit_x
end

function logit(x::Float64)
    return log(x / (1 - x))
end

################################################################################
## transforming the dataframe
################################################################################

function transform_data!(df::DataFrame, transform::Function;
                        inds::AbstractVector{Int} = 1:(size(df, 2) - 1))
    nms = names(df)
    for ind in inds
        CURRENT_NAME[1] = nms[ind + 1]
        col = df[!, ind + 1]
        col .= transform(col)
    end
end

function log_transform_data!(df::DataFrame;
                             inds::AbstractVector{Int} = 1:(size(df, 2) - 1))
    transform_data!(df, log_transform, inds = inds)
end

function standardize_data!(df::DataFrame)
    for i = 2:size(df, 2)
        col = df[findall(!isnan, df[!, i]), i]
        μ = mean(col)
        σ = std(col, mean=μ, corrected=false)
        df[!, i] .-= μ
        df[!, i] ./= σ
    end
end