text:
"""
    InputShapeVertex <: AbstractVertex
    InputShapeVertex(v::AbstractVertex, t::FluxLayer)

Input type vertex which also has information about what type of layer the input is shaped for.
"""
struct InputShapeVertex{V<:AbstractVertex, L<:FluxLayer} <: AbstractVertex
    base::V
    t::L
end

@functor InputShapeVertex

const inputshapemotivation = """
Providing the input type is not strictly necessary for the package to work and in many cases a normal `inputvertex`
will do.

One example of when it is useful is the [`concat`](@ref) function which needs to know the input type to
automatically determine which dimension to concatenate.
"""

"""
    inputvertex(name, size, type::FluxLayer)

Return an immutable input type vertex with the given `name` and `size` and a `type` which can be used to
indicate what type of input is expected.

$inputshapemotivation
"""
NaiveNASlib.inputvertex(name, size, type::FluxLayer) = InputShapeVertex(inputvertex(name, size), type)

"""
    convinputvertex(name, nchannel, ndim)

Return an input type vertex with the given `name` which promises convolution shaped input
with `nchannel` channels and `ndim` number of dimensions for feature maps (e.g. 2 for images)
suitable for `Flux`s convolution layers.

$inputshapemotivation
"""
convinputvertex(name, nchannel, ndim) = inputvertex(name, nchannel, GenericFluxConvolutional{ndim}())

doc_convninputvertex(n) = """
    conv$(n)dinputvertex(name, nchannel)

Return an input type vertex with the given `name` which promises convolution shaped input
with `nchannel` channels suitable for `Flux`s convolution layers.

Equivalent to [`convinputvertex(name, nchannel, $(n))`](@ref).

$inputshapemotivation
"""

@doc doc_convninputvertex(1)
conv1dinputvertex(name, nchannel) = convinputvertex(name, nchannel, 1)
@doc doc_convninputvertex(2)
conv2dinputvertex(name, nchannel) = convinputvertex(name, nchannel, 2)
@doc doc_convninputvertex(3)
conv3dinputvertex(name, nchannel) = convinputvertex(name, nchannel, 3)

"""
    denseinputvertex(name, size)

Return an input type vertex with the given `name` which promises 2D shaped input
with `size` number of features suitable for e.g. `Flux`s `Dense` layer.

$inputshapemotivation
"""
denseinputvertex(name, size) = inputvertex(name, size, GenericFlux2D())

"""
    rnninputvertex(name, size)

Return an input type vertex with the given `name` which promises 2D shaped input
with `size` number of features suitable for `Flux`s recurrent layers.

$inputshapemotivation
"""
rnninputvertex(name, size) = inputvertex(name, size, GenericFluxRecurrent())

layertype(v::InputShapeVertex) = v.t
layer(v::InputShapeVertex) = LayerTypeWrapper(v.t) # so that layertype(layer(v)) works
NaiveNASlib.base(v::InputShapeVertex) = v.base
NaiveNASlib.name(v::InputShapeVertex) = name(base(v))
NaiveNASlib.nout(v::InputShapeVertex) = nout(base(v))
NaiveNASlib.nin(v::InputShapeVertex) = nin(base(v))
NaiveNASlib.outputs(v::InputShapeVertex) = outputs(base(v))
NaiveNASlib.inputs(::InputShapeVertex) = []

# Only to prevent stack overflow above
struct LayerTypeWrapper{L}
    t::L
end
layertype(l::LayerTypeWrapper) = l.t

"""
    SizeNinNoutConnected <: NaiveNASlib.DecoratingTrait
    SizeNinNoutConnected(t)

Trait for computations for which a change in output size results in a change in input size but which
is not fully `SizeTransparent`.

Example of this is DepthWiseConv where output size must be an integer multiple of the input size.

Does not create any constraints or objectives, only signals that vertices after a
`SizeNinNoutConnected` might need to change size if the size of the `SizeNinNoutConnected` vertex changes.
"""
struct SizeNinNoutConnected{T <: NaiveNASlib.MutationTrait} <: NaiveNASlib.DecoratingTrait
    base::T
end
NaiveNASlib.base(t::SizeNinNoutConnected) = t.base

NaiveNASlib.all_in_Δsize_graph(mode, ::SizeNinNoutConnected, args...) = NaiveNASlib.all_in_Δsize_graph(mode, SizeInvariant(), args...)

Flux.trainable(v::CompVertex) = Flux.trainable(v.computation)
Flux.trainable(g::CompGraph) = Flux.trainable(vertices(g))

const doc_layerfun_and_traitfun = """
Keyword argument `layerfun` can be used to wrap the computation, e.g. in an [`ActivationContribution`](@ref).

Keyword argument `traitfun` can be used to wrap the `MutationTrait` of the vertex in a `DecoratingTrait`
"""

"""
    fluxvertex(l, in::AbstractVertex; layerfun=LazyMutable, traitfun=validated())

Return a vertex which wraps the layer `l` and has input vertex `in`.

$doc_layerfun_and_traitfun
"""
fluxvertex(l, in::AbstractVertex; layerfun=LazyMutable, traitfun=validated()) = fluxvertex(layertype(l), l, in, layerfun, traitfun)

"""
    fluxvertex(name::AbstractString, l, in::AbstractVertex; layerfun=LazyMutable, traitfun=validated())

Return a vertex with name `name` which wraps the layer `l` and has input vertex `in`.

Name is only used when displaying or logging and does not have to be unique (although it probably is a good idea).

$doc_layerfun_and_traitfun
"""
fluxvertex(name::AbstractString, l, in::AbstractVertex; layerfun=LazyMutable, traitfun=validated()) = fluxvertex(layertype(l), l, in, layerfun, traitfun ∘ named(name))

fluxvertex(::FluxParLayer, l, in::AbstractVertex, layerfun, traitfun) = absorbvertex(layerfun(MutableLayer(l)), in, traitdecoration = traitfun)

fluxvertex(::FluxDepthwiseConv, l, in::AbstractVertex, layerfun, traitfun) = absorbvertex(layerfun(MutableLayer(l)), in; traitdecoration=traitfun ∘ SizeNinNoutConnected)

fluxvertex(::FluxParInvLayer, l, in::AbstractVertex, layerfun, traitfun) = invariantvertex(layerfun(MutableLayer(l)), in, traitdecoration=traitfun ∘ FixedSizeTrait)

fluxvertex(::FluxNoParLayer, l, in::AbstractVertex, layerfun, traitfun) = invariantvertex(layerfun(NoParams(l)), in, traitdecoration=traitfun)

"""
    concat(v::AbstractVertex, vs::AbstractVertex...; traitfun=identity, layerfun=identity)

Return a vertex which concatenates input along the activation (e.g. channel if convolution, first dimension if dense) dimension.

Inputs must have compatible activation shapes or an exception will be thrown.

$doc_layerfun_and_traitfun

See also `NaiveNASlib.conc`.
"""
function concat(v::AbstractVertex, vs::AbstractVertex...; traitfun=identity, layerfun=identity)
    allactdims = unique(mapreduce(actdim, vcat, [v, vs...]))
    if length(allactdims) != 1
        throw(DimensionMismatch("Can not concatenate activations with different shapes! Got: $(join(allactdims, ", ", " and "))"))
    end

    allactranks = unique(mapreduce(actrank, vcat, [v, vs...]))
    if length(allactranks) != 1
        throw(DimensionMismatch("Can not concatenate activations with different shapes! Got:  $(join(allactranks, ", ", " and "))"))
    end

    conc(v, vs...; dims=allactdims[], traitdecoration=traitfun, outwrap=layerfun)
end

"""
    concat(name::AbstractString, v::AbstractVertex, vs::AbstractVertex...; traitfun=identity, layerfun=identity)

Return a vertex with name `name` which concatenates input along the activation (e.g. channel if convolution, first dimension if dense) dimension.

Name is only used when displaying or logging and does not have to be unique (although it probably is a good idea).

Inputs must have compatible activation shapes or an exception will be thrown.

$doc_layerfun_and_traitfun

See also `NaiveNASlib.conc`.
"""
concat(name::AbstractString, v::AbstractVertex, vs::AbstractVertex...; traitfun = identity, layerfun=identity) = concat(v, vs..., traitfun=traitfun ∘ named(name), layerfun=layerfun)

"""
    layer(v)

Return the computation wrapped inside `v` and inside any mutable wrappers.

# Examples
```jldoctest
julia> using NaiveNASflux, Flux

julia> layer(fluxvertex(Dense(2,3), inputvertex("in", 2)))
Dense(2, 3)         # 9 parameters
```
"""
layer(v::AbstractVertex) = layer(base(v))
layer(v::CompVertex) = layer(v.computation)
layer(::InputVertex) = nothing

layertype(v::AbstractVertex) = layertype(base(v))
layertype(v::CompVertex) = layertype(v.computation)
layertype(::InputVertex) = nothing

actdim(v::AbstractVertex) = actdim.(layer.(NaiveNASlib.findterminating(v, inputs)))
actrank(v::AbstractVertex) = actrank.(layer.(NaiveNASlib.findterminating(v, inputs)))

mutate_weights(v::AbstractVertex, w) = mutate_weights(base(v), w)
mutate_weights(v::CompVertex, w) = mutate_weights(v.computation, w)

"""
    setlayer!(x, propval)

Set the properties `propval` to the layer wrapped in `x` where `propval` is a named tuple with fieldname->value pairs.

This typically means create a new layer with the given values and set the wrapped layer to it.

### Examples

```julia-repl
julia> v = fluxvertex(Dense(3, 4, relu), inputvertex("in", 3));

julia> layer(v)
Dense(3, 4, relu)   # 16 parameters

julia> NaiveNASflux.setlayer!(v, (;σ=tanh));

julia> layer(v)
Dense(3, 4, tanh)   # 16 parameters
```
"""
function setlayer!(x, propval) end
setlayer!(v::AbstractVertex, propval) = setlayer!(base(v), propval)
setlayer!(v::CompVertex, propval) = setlayer!(v.computation, propval)
setlayer!(m::AbstractMutableComp, propval) = setlayer!(wrapped(m), propval)
setlayer!(m::ResetLazyMutable, propval) = setlayer!(m.wrapped, propval)
setlayer!(m::MutationTriggered, propval) = setlayer!(m.wrapped, propval)
function setlayer!(m::MutableLayer, propval)
    m.layer = setproperties(m.layer, propval)
end

ext: jl


text:
module SPDEPrecisionMatrices

using TriangleMesh
using Interpolations
using SparseArrays
using LinearAlgebra
using SpecialFunctions
using Statistics
using NearestNeighbors
using ConcaveHull

export component_matrices,
    unscaled_precision_matrix,
    calculate_κ,
    calculate_τ,
    precision_matrix,
    make_border_points,
    observation_matrix

include("mesh_utils.jl")
include("spde_precision_matrices.jl")

end # module


text:
#!/bin/bash
#=
JULIA="${JULIA:-julia}"
JULIA_CMD="${JULIA_CMD:-$JULIA --color=yes --startup-file=no}"
# below gets the directory name of the script, even if there is a symlink involved
# from https://stackoverflow.com/questions/59895/get-the-source-directory-of-a-bash-script-from-within-the-script-itself
SOURCE="${BASH_SOURCE[0]}"
while [ -h "$SOURCE" ]; do # resolve $SOURCE until the file is no longer a symlink
  DIR="$( cd -P "$( dirname "$SOURCE" )" >/dev/null 2>&1 && pwd )"
  SOURCE="$(readlink "$SOURCE")"
  [[ $SOURCE != /* ]] && SOURCE="$DIR/$SOURCE" # if $SOURCE was a relative symlink, we need to resolve it relative to the path where the symlink file was located
done
DIR="$( cd -P "$( dirname "$SOURCE" )" >/dev/null 2>&1 && pwd )"
export JULIA_PROJECT=$DIR/..
export JULIA_LOAD_PATH=@:@stdlib  # exclude default environment
exec $JULIA_CMD -e 'include(popfirst!(ARGS))' "$SOURCE" "$@"
=#
using ArgParse
s = ArgParseSettings(description="Coallate ljh files and convert to off with given model file\n"*
"example usage:\n"*
"./ljh2off.jl 20181205_B/20181205_B_model.hdf5 20181205_ --endings=B C\n"
)
@add_arg_table s begin
    "model_file"
        help = "name of the hdf5 with basis to use"
        required = true
        arg_type = String
    "ljh_file"
        help = "name of the pulse files containing ljh file to use to make basis, will process all channels"
        required = true
        arg_type = String
    "--endings"
        nargs = '+'
        help = "a series of endings appended to ljh_file, eg A B C"
        default = String[""]
        arg_type = String
    "--maxchannels"
        default = 10000
        help = "process at most this many channels, in order from lowest channel number"
        arg_type = Int
    "--outdir"
        arg_type = String
        help ="path to the output directory, otherwise calculate automatically"
    "--replaceoutput", "-r"
        help = "delete and replace any existing output files"
        action = :store_true
end
parsed_args = parse_args(ARGS, s)
if parsed_args["outdir"] == nothing
    parts = splitpath(parsed_args["ljh_file"])
    parts[end] = "off"
    ljhpathPlusoff = joinpath(parts...)
    parsed_args["outdir"] = ljhpathPlusoff*prod(parsed_args["endings"])
end
display(parsed_args);println()
using Pope
using Pope.LJH
using HDF5
using JSON
import Statistics

function off_header_json(ljh::LJH.LJHFile, z::Pope.SVDBasisWithCreationInfo)
    offVersion = "0.2.0" # version 0.2.0 has projectors as binary after the header
    offHeader = Dict(
        "FileFormatVersion" => offVersion,
        "FramePeriodSeconds" => LJH.frameperiod(ljh),
        "NumberOfBases" => size(z.svdbasis.projectors)[1],
        "FileFormat" => "OFF"
    )
    offHeader["ModelInfo"]=Dict(
        "Projectors" =>Dict(
            "Rows" =>size(z.svdbasis.projectors)[1],
            "Cols" =>size(z.svdbasis.projectors)[2],
            "SavedAs" => "float64 binary data after header and before records. projectors first then basis, nbytes = rows*cols*8 for each projectors and basis"
    ),
        "Basis" =>Dict(
            "Rows" =>size(z.svdbasis.basis)[1],
            "Cols" =>size(z.svdbasis.basis)[2],
            "SavedAs" => "float64 binary data after header and before records. projectors first then basis, nbytes = rows*cols*8 for each projectors and basis"
    ),
        "NoiseStandardDeviation" => z.noise_std_dev,
        "NoiseModelFile" => z.noise_model_file,
        "PulseFile" => z.pulse_file,
        "ModelFile" => abspath(parsed_args["ljh_file"])
    )
    offHeader["ReadoutInfo"] = Dict(
        "ColumnNum" => LJH.column(ljh),
        "RowNum" => LJH.row(ljh),
        "NumberOfColumns" => ljh.headerdict["Number of columns"],
        "NumberOfRows" => ljh.headerdict["Number of rows"],
    )
    offHeader["CreationInfo"]=Dict(
        "SourceName" => "ljh2off.jl"
    )
    return JSON.json(offHeader,4) # includes "\n" termination
end

function write_off_header(io, ljh, z)
    nheader = write(io, off_header_json(ljh, z)) #  includes "\n" termination
    # write the projectors as float64 row-major (julia arrays are column major) binary after the header
    nprojectors = write(io, transpose(Float64.(z.svdbasis.projectors)))
    nbasis = write(io, transpose(Float64.(z.svdbasis.basis)))
end

ljh_files = [parsed_args["ljh_file"]*ending for ending in parsed_args["endings"]]
@show ljh_files
@show parsed_args["endings"]

earliest_timestamp_usec = [typemax(Int64) for i in ljh_files]
latest_timestamp_usec = [typemin(Int64) for i in ljh_files]
# make sure the directory exists
outputDir = parsed_args["outdir"]
offPrefix = splitpath(outputDir)[end]
if isdir(outputDir)
    if parsed_args["replaceoutput"]
        rm(outputDir, force=true, recursive=true)
    else
        println("outdir $(outputDir) exists, pass --replaceoutput to overwrite")
        println("exiting")
        exit(1)
    end
end
mkdir(outputDir)
HDF5.h5open(parsed_args["model_file"],"r") do h5
    channels_processed = 0
    channums = sort(parse.(Int,names(h5))) # parse to int early so we can sort
    for channum in channums
        channels_processed >= parsed_args["maxchannels"] && break
        offPrefix = splitpath(outputDir)[end]
        offFilename = joinpath(outputDir, offPrefix*"_chan$(channum).off")
        z = Pope.hdf5load(Pope.SVDBasisWithCreationInfo,h5["$channum"])
        open(offFilename,"w") do f
            for (i,ljhbase) in enumerate(ljh_files)
                dir, base, ext = LJH.dir_base_ext(ljhbase)
                ljh_filename = joinpath(dir,base)*"_chan$(channum)$(ext)"
                ljh = try
                    LJH.LJHFile(ljh_filename)
                catch ex
                    println("Channel $channum ljh file failed to open: $ljh_filename")
                    println(ex)
                    continue
                end
                println("$ljh_filename => $offFilename")
                if i==1 #write header when processing first ljh file
                    write_off_header(f,ljh,z)
                    header_written = true
                end
                earliest_timestamp_usec[i] = min(ljh[1].timestamp_usec, earliest_timestamp_usec[i])
                latest_timestamp_usec[i] = max(ljh[end].timestamp_usec, latest_timestamp_usec[i])
                channels_processed += 1
                for record in ljh
                    dp = Pope.record2dataproduct(z.svdbasis,record)
                    write(f,Int32(dp.nsamples), Int32(dp.first_rising_sample-4), Int64(dp.frame1index),
                    Int64(dp.timestamp_usec*1000), Float32(Statistics.mean(record.data[1:dp.first_rising_sample-3])),
                    Float32(dp.residual_std), Float32.(dp.reduced))
                end # end of loop over ljh
            end # end of loop over parsed_args["ljh_files"]
        end # end of open do block
    end # end of loop over channels
end # end of h5open do block

experimental_state_filename = joinpath(outputDir, offPrefix*"_experiment_state.txt")
state_names = parsed_args["endings"]
@show state_names
@show experimental_state_filename
open(experimental_state_filename,"w") do f
    write(f,"# unix time in nanoseconds, state label\n")
    write(f,string(earliest_timestamp_usec[1]*1000),", ", "START\n")
    for (i, timestamp_usec) in enumerate(earliest_timestamp_usec)
        write(f, string(timestamp_usec*1000),", ",state_names[i],"\n")
    end
    write(f,string(latest_timestamp_usec[end]*1000),", ", "END\n")
end


text:
import Base.vcat, Base.hcat, Base.hvcat, Base.start
export vcat, hcat, hvcat

function vcat(args::AffineExpr...)
  num_cols = args[1].size[2]
  num_rows = 0
  # Calculate numbers of rows after vertical concatenation
  for arg in args
    if arg.size[2] != num_cols
      error("All arguments must have same number of columns")
    end
    num_rows += arg.size[1]
  end

  vars_to_coeffs_map = Dict{UInt64, Constant}()
  vec_sz = num_rows * num_cols
  constant = Constant(spzeros(vec_sz, 1))
  index_start = 0

  for i = 1:length(args)
    for (v, c) in args[i].vars_to_coeffs_map
      if !haskey(vars_to_coeffs_map, v)
        vars_to_coeffs_map[v] = Constant(spzeros(vec_sz, c.size[2]))
      end
      for j in 1 : num_cols
        index = index_start + (j - 1) * num_rows
        vars_to_coeffs_map[v].value[index + 1 : index + args[i].size[1], :] =
            c.value[(j - 1) * args[i].size[1] + 1 : j * args[i].size[1], :]
      end
    end

    for j in 1 : num_cols
      index = index_start + (j - 1) * num_rows
      constant.value[index + 1 : index + args[i].size[1]] =
          args[i].constant.value[(j - 1) * args[i].size[1] + 1 : j * args[i].size[1]]
    end

    index_start += args[i].size[1]
  end

  this = AffineExpr(:vcat, tuple(args...), vars_to_coeffs_map, constant, (num_rows, num_cols))
  this.evaluate = ()->Base.vcat([arg.evaluate() for arg in args]...)
  return this
end

function hcat(args::AffineExpr...)
  num_rows = args[1].size[1]
  num_cols = 0

  for arg in args
    if arg.size[1] != num_rows
      error("All arguments must have same number of rows")
    end
    num_cols += arg.size[2]
  end

  vars_to_coeffs_map = Dict{UInt64, Constant}()
  vec_sz = num_rows * num_cols
  constant = Constant(spzeros(vec_sz, 1))
  index = 0

  for i = 1:length(args)
    arg_vec_sz = args[i].size[1] * args[i].size[2]

    for (v, c) in args[i].vars_to_coeffs_map
      if !haskey(vars_to_coeffs_map, v)
        vars_to_coeffs_map[v] = Constant(spzeros(vec_sz, c.size[2]))
      end
      vars_to_coeffs_map[v].value[index + 1 : index + arg_vec_sz, :] = c.value
    end
    constant.value[index + 1 : index + arg_vec_sz] = args[i].constant.value

    index += arg_vec_sz
  end

  this = AffineExpr(:hcat, tuple(args...), vars_to_coeffs_map, constant, (num_rows, num_cols))
  this.evaluate = ()->Base.hcat([arg.evaluate() for arg in args]...)
  return this
end

function hvcat(rows::Tuple{Vararg{Int64}}, args::AffineExpr...)
  row_exprs = AffineExpr[]
  index = 0
  for row_size in rows
    push!(row_exprs, hcat(args[index + 1 : index + row_size]...))
    index += row_size
  end
  return vcat(row_exprs...)
end


text:
include("header.jl")

#=
These Jupyter notebooks and associated code serve as a companion to the course.

These are not static notebooks. **You are encouraged to play with them: run with different parameter values, create other plots, etc.!** The whole point is to enable an understanding of fluid mechanics through interactiveness.

*The code is released under the [MIT license](https://opensource.org/licenses/MIT).*
<hr />
=#

#=
## Table of Notebooks

- [Introduction](0.1-WelcomeToJulia.ipynb)

#### Chapter 1 - Fluid properties, flow quantities, and dimensional analysis

- [1.0 Units and quantities](1.0-UnitsAndQuantities.ipynb)
- [1.1 Velocity profiles](1.1-VelocityProfiles.ipynb)
- [1.2 Dimensionless quantities](1.2-DimensionlessQuantities.ipynb)
- [1.3 Dimensional analysis](1.3-DimensionalAnalysis.ipynb)

#### Chapter 2 - Fluid statics

- [2.0 Basics of fluid statics](2.0-FluidStatics.ipynb)
- [2.1 Manometry](2.1-Manometry.ipynb)
- [2.2 Submerged surfaces](2.2-SubmergedSurfaces.ipynb)

#### Chapter 3 - Flow kinematics

- [3.0 Streamlines](3.0-Streamlines.ipynb)

#### Chapter 4 - Control volume analysis

-

#### Chapter 5 - Pipe flow systems

- [5.0 Energy equation example](5.0-EnergyEquationExample.ipynb)

#### Chapter 6 - Compressible flow

- [6.0 Basic gas dynamics](6.0-BasicGasDynamics.ipynb)
- [6.1 Isentropic gas dynamics](6.1-IsentropicGasDynamics.ipynb)
- [6.2 Converging-diverging nozzles](6.2-ConvergingDivergingNozzle.ipynb)

=#


text:
# --------------------------------------------------
# Helper functions

function closest_index(arr::AbstractArray, val)
    arr .- val .|> abs |> argmin
end

function cart_index_to_array(idx::Array{CartesianIndex{3}, 1})
    a = zeros(Float32, (3, length(idx)));
    Threads.@threads for i=1:length(idx)
        a[:, i] .= idx[i].I
    end
    return a
end

function cart_index_to_array(idx::Array{CartesianIndex{2}, 1})
    a = zeros(Float32, (2, length(idx)));
    Threads.@threads for i=1:length(idx)
        a[:, i] .= idx[i].I
    end
    return a
end

function closest_index(arr::AbstractArray, val)
    arr .- val .|> abs |> argmin
end

function to_kps_val(x)
    uconvert(u"km / s", x).val
end


text:
@safetestset "SVG API support in Renderer.Html" begin

  @safetestset "" begin
    using Genie
    using Genie.Renderer.Html

    @test_throws UndefVarError svg()
  end;

  @safetestset "Loading SVG support makes SVG API available" begin
    using Genie
    using Genie.Renderer.Html

    Html.register_svg_elements()
    @test svg() == "<svg></svg>"

    @test_throws UndefVarError clippath()
    @test clipPath() == "<clipPath></clipPath>"
  end;

end;


text:
using SparseMatrix
using LinearAlgebra
using SparseArrays
using Test

@testset "nnz1" begin
    for i=1:100
        m = rand(1:20)
        n = rand(1:20)
        p = rand()
        A=Matrix(sprandn(m,n,p))
        x1=sparse(A)
        x2=SparseCSR(A)
        x3=SparseCSC(A)
        x4=SparseCOO(A)
        @test nnz(x1) == nnz(x2) == nnz(x3) == nnz(x4)
    end
end

@testset "tocsc1" begin
    for i=1:100
        m = rand(1:20)
        n = rand(1:20)
        p = rand()
        A=Matrix(sprandn(m,n,p))
        x0 = sparse(A)
        x1=SparseCSR(A)
        x2=SparseCSC(A)
        x3=SparseCOO(A)
        y1=sparse(x1)
        y2=sparse(x2)
        y3=sparse(x3)
        if nnz(x1) != 0
            x1.val[1] = 0.0
            x2.val[1] = 0.0
            x3.val[1] = 0.0
        end
        @test x0 == y1 == y2 == y3
    end
end

@testset "multiple1" begin
    for i=1:100
        m = rand(1:20)
        n = rand(1:20)
        p = rand()
        A=Matrix(sprandn(m,n,p))
        x=randn(n)
        C1 = A * x
        C2 = SparseCSR(A) * x
        C3 = SparseCSC(A) * x
        C4 = SparseCOO(A) * x
        @test C1 ≈ C2 ≈ C3 ≈ C4
    end
end

@testset "multiple2" begin
    for i=1:100
        m = rand(1:20)
        n = rand(1:20)
        p = rand()
        A=Matrix(sprandn(n,m,p))
        x=randn(n)
        C1 = A' * x
        C2 = SparseCSR(A)' * x
        C3 = SparseCSC(A)' * x
        C4 = SparseCOO(A)' * x
        @test C1 ≈ C2 ≈ C3 ≈ C4
    end
end

@testset "multiple3" begin
    for i=1:100
        m = rand(1:20)
        n = rand(1:20)
        k = rand(1:20)
        p = rand()
        A=Matrix(sprandn(m,k,p))
        B=randn(k,n)
        C1 = A * B
        C2 = SparseCSR(A) * B
        C3 = SparseCSC(A) * B
        C4 = SparseCOO(A) * B
        @test C1 ≈ C2 ≈ C3 ≈ C4
    end
end

@testset "multiple4" begin
    for i=1:100
        m = rand(1:20)
        n = rand(1:20)
        k = rand(1:20)
        p = rand()
        A=Matrix(sprandn(k,m,p))
        B=randn(k,n)
        C1 = A' * B
        C2 = SparseCSR(A)' * B
        C3 = SparseCSC(A)' * B
        C4 = SparseCOO(A)' * B
        @test C1 ≈ C2 ≈ C3 ≈ C4
    end
end

@testset "multiple5" begin
    for i=1:100
        m = rand(1:20)
        n = rand(1:20)
        k = rand(1:20)
        p = rand()
        A=Matrix(sprandn(m,k,p))
        B=randn(n,k)
        C1 = A * B'
        C2 = SparseCSR(A) * B'
        C3 = SparseCSC(A) * B'
        C4 = SparseCOO(A) * B'
        @test C1 ≈ C2 ≈ C3 ≈ C4
    end
end

@testset "multiple6" begin
    for i=1:100
        m = rand(1:20)
        n = rand(1:20)
        k = rand(1:20)
        p = rand()
        A=Matrix(sprandn(k,m,p))
        B=randn(n,k)
        C1 = A' * B'
        C2 = SparseCSR(A)' * B'
        C3 = SparseCSC(A)' * B'
        C4 = SparseCOO(A)' * B'
        @test C1 ≈ C2 ≈ C3 ≈ C4
    end
end

@testset "multiple7" begin
    for i=1:100
        m = rand(1:20)
        n = rand(1:20)
        k = rand(1:20)
        p = rand()
        A=randn(m,k)
        B=Matrix(sprandn(k,n,p))
        C1 = A * B
        C2 = A * SparseCSR(B)
        C3 = A * SparseCSC(B)
        C4 = A * SparseCOO(B)
        @test C1 ≈ C2 ≈ C3 ≈ C4
    end
end

@testset "multiple8" begin
    for i=1:100
        m = rand(1:20)
        n = rand(1:20)
        k = rand(1:20)
        p = rand()
        A=randn(k,m)
        B=Matrix(sprandn(k,n,p))
        C1 = A' * B
        C2 = A' * SparseCSR(B)
        C3 = A' * SparseCSC(B)
        C4 = A' * SparseCOO(B)
        @test C1 ≈ C2 ≈ C3 ≈ C4
    end
end

@testset "multiple9" begin
    for i=1:100
        m = rand(1:20)
        n = rand(1:20)
        k = rand(1:20)
        p = rand()
        A=randn(m,k)
        B=Matrix(sprandn(n,k,p))
        C1 = A * B'
        C2 = A * SparseCSR(B)'
        C3 = A * SparseCSC(B)'
        C4 = A * SparseCOO(B)'
        @test C1 ≈ C2 ≈ C3 ≈ C4
    end
end

@testset "multiple10" begin
    for i=1:100
        m = rand(1:20)
        n = rand(1:20)
        k = rand(1:20)
        p = rand()
        A=randn(k,m)
        B=Matrix(sprandn(n,k,p))
        C1 = A' * B'
        C2 = A' * SparseCSR(B)'
        C3 = A' * SparseCSC(B)'
        C4 = A' * SparseCOO(B)'
        @test C1 ≈ C2 ≈ C3 ≈ C4
    end
end

@testset "diag1" begin
    for i = 1:100
        m = rand(1:20)
        n = rand(1:20)
        p = rand()
        A = Matrix(sprandn(m,n,p))
        d1 = spdiag(A)
        d2 = convert(Vector{Float64}, spdiag(SparseCSR(A)))
        d3 = convert(Vector{Float64}, spdiag(SparseCSC(A)))
        d4 = convert(Vector{Float64}, spdiag(SparseCOO(A)))
        @test d1 ≈ d2 ≈ d3 ≈ d4
    end
end

@testset "diag2" begin
    for i = 1:100
        m = rand(1:20)
        n = rand(1:20)
        p = rand()
        A = Matrix(sprandn(m,n,p))
        d1 = spdiag(A)
        d2 = spdiag(SparseCSR(A))
        d3 = spdiag(SparseCSC(A))
        d4 = spdiag(SparseCOO(A))
        for k = 1:length(d1)
            @test d1[k] ≈ d2[k] ≈ d3[k] ≈ d4[k]
        end
    end
end

@testset "diag3" begin
    for i = 1:100
        m = rand(1:20)
        n = rand(1:20)
        p = rand()
        A = Matrix(sprandn(m,n,p))
        d1 = abs.(spdiag(A))
        d2 = abs.(spdiag(SparseCSR(A)))
        d3 = abs.(spdiag(SparseCSC(A)))
        d4 = abs.(spdiag(SparseCOO(A)))
        @test d1 ≈ d2 ≈ d3 ≈ d4
    end
end

@testset "diag4" begin
    for i = 1:100
        m = rand(1:20)
        n = rand(1:20)
        p = rand()
        A = Matrix(sprandn(m,n,p))
        d1 = maximum(spdiag(A))
        d2 = maximum(spdiag(SparseCSR(A)))
        d3 = maximum(spdiag(SparseCSC(A)))
        d4 = maximum(spdiag(SparseCOO(A)))
        @test d1 ≈ d2 ≈ d3 ≈ d4
    end
end

@testset "copy" begin
    for i = 1:100
        m = rand(1:20)
        n = rand(1:20)
        p = rand()
        A = Matrix(sprandn(m,n,p))
        d1 = Matrix(copy(SparseCSR(A)))
        d2 = Matrix(copy(SparseCSC(A)))
        d3 = Matrix(copy(SparseCOO(A)))
        @test A ≈ d1 ≈ d2 ≈ d3
    end
end

@testset "zero" begin
    for i = 1:100
        m = rand(1:20)
        n = rand(1:20)
        p = rand()
        A = Matrix(sprandn(m,n,p))
        z = zero(A)
        d1 = Matrix(zero(SparseCSR(A)))
        d2 = Matrix(zero(SparseCSC(A)))
        d3 = Matrix(zero(SparseCOO(A)))
        d4 = Matrix(zero(SparseArrays.SparseMatrixCSC(A)))
        @test z ≈ d1 ≈ d2 ≈ d3 ≈ d4
    end
end

@testset "linearlize" begin
    for i = 1:100
        m = rand(1:20)
        n = rand(1:20)
        p = rand()
        A = SparseCSC(sprandn(m,n,p))
        B = copy(A)
        for i = 1:length(B)
            B[i] *= 100.0
        end
        @test A.val * 100.0 ≈ B.val
        A = SparseCSC(sprandn(m,n,p))
        B = copy(A)
        for i = eachindex(B)
            B[i] *= 100.0
        end
        @test A.val * 100.0 ≈ B.val
    end
end

@testset "block 1" begin
    for i = 1:100
        m = rand(1:20)
        n = rand(1:20)
        p = rand()
        A1 = SparseCSC(sprandn(m,n,p))
        x1 = rand(n,1)
        m = rand(1:20)
        n = rand(1:20)
        p = rand()
        A2 = SparseCSC(sprandn(m,n,p))
        x2 = rand(n,1)

        y0 = [A1 * x1; A2 * x2]

        X = BlockCOO(2,2, [(1,1,A1), (2,2,A2)])
        y = SparseCOO(X) * [x1; x2]

        @test Matrix(SparseCOO(X)) == Matrix(SparseCSR(X)) == Matrix(SparseCSC(X)) == Matrix(sparse(X))
        @test y0 ≈ y
    end
end

@testset "block 2" begin
    for i = 1:100
        m = rand(1:20)
        n = rand(1:20)
        p = rand()
        A1 = SparseCSR(sprandn(m,n,p))
        x1 = rand(n,1)
        m = rand(1:20)
        n = rand(1:20)
        p = rand()
        A2 = SparseCOO(sprandn(m,n,p))
        x2 = rand(n,1)

        y0 = [A1 * x1; A2 * x2]

        X = BlockCOO(2,2, [(1,1,A1), (2,2,A2)])
        y = SparseCOO(X) * [x1; x2]

        @test Matrix(SparseCOO(X)) == Matrix(SparseCSR(X)) == Matrix(SparseCSC(X)) == Matrix(sparse(X))
        @test y0 ≈ y
    end
end

@testset "block 3" begin
    for i = 1:100
        m = rand(1:20)
        n = rand(1:20)
        p = rand()
        A1 = Matrix(sprandn(m,n,p))
        x1 = rand(n,1)
        m = rand(1:20)
        n = rand(1:20)
        p = rand()
        A2 = SparseCSC(sprandn(m,n,p))
        x2 = rand(n,1)

        y0 = [A1 * x1; A2 * x2]

        X = BlockCOO(2,2, [(1,1,A1), (2,2,A2)])
        y = SparseCOO(X) * [x1; x2]

        @test Matrix(SparseCOO(X)) == Matrix(SparseCSR(X)) == Matrix(SparseCSC(X)) == Matrix(sparse(X))
        @test y0 ≈ y
    end
end

@testset "spdiag 1" begin
    for ite = 1:100
        m = rand(1:20)
        n = rand(1:20)
        p = rand()
        A = Matrix(sprandn(m,n,p))
        d0 = spdiag(A)
        d1 = spdiag(SparseCSR(A))
        d2 = spdiag(SparseCSC(A))
        d3 = spdiag(SparseCOO(A))
        d4 = spdiag(SparseMatrixCSC(A))
        v0 = [d0.val[i] for i = d0.index if !iszero(d0.val[i])]
        v1 = [d1.val[i] for i = d1.index if !iszero(i)]
        v2 = [d2.val[i] for i = d2.index if !iszero(i)]
        v3 = [d3.val[i] for i = d3.index if !iszero(i)]
        v4 = [d4.val[i] for i = d4.index if !iszero(i)]
        @test v0 == v1 == v2 == v3 == v4
    end
end

@testset "block 4" begin
    for i = 1:100
        m = rand(1:20)
        n = rand(1:20)
        p = rand()
        A1 = Matrix(sprandn(m,n,p))
        m = rand(1:20)
        n = rand(1:20)
        p = rand()
        A2 = SparseCSC(sprandn(m,n,p))
        X = spzeros(AbstractMatrix{Float64}, 2, 2)
        X[1,1] = A1
        X[2,2] = A2
        if sum(A1 .!= 0) != 0 && SparseArrays.nnz(A2) != 0
            AA = Matrix(sparse(block(X)))
            A1 = Matrix(A1)
            A2 = Matrix(A2)
            BB = zeros(size(A1) .+ size(A2))
            BB[1:size(A1)[1],1:size(A1)[2]] = A1
            BB[size(A1)[1]+1:end,size(A1)[2]+1:end] = A2
            @test AA == BB
        end
    end
end


text:
@doc raw"""
    truncated_conjugate_gradient_descent(M, F, gradF, x, η, HessF, trust_region_radius)

solve the trust-region subproblem

```math
\operatorname*{arg\,min}_{η ∈ T_xM}
m_x(η) \quad\text{where}
m_x(η) = F(x) + ⟨\operatorname{grad}F(x),η⟩_x + \frac{1}{2}⟨\operatorname{Hess}F(x)[η],η⟩_x,
```

```math
\text{such that}\quad ⟨η,η⟩_x ≤ Δ^2
```

with the [`truncated_conjugate_gradient_descent`](@ref).
For a description of the algorithm and theorems offering convergence guarantees,
see the reference:

* P.-A. Absil, C.G. Baker, K.A. Gallivan,
    Trust-region methods on Riemannian manifolds, FoCM, 2007.
    doi: [10.1007/s10208-005-0179-9](https://doi.org/10.1007/s10208-005-0179-9)
* A. R. Conn, N. I. M. Gould, P. L. Toint, Trust-region methods, SIAM,
    MPS, 2000. doi: [10.1137/1.9780898719857](https://doi.org/10.1137/1.9780898719857)

# Input

* `M` – a manifold ``\mathcal M``
* `F` – a cost function ``F: \mathcal M → ℝ`` to minimize
* `gradF` – the gradient ``\operatorname{grad}F: \mathcal M → T\mathcal M`` of `F`
* `HessF` – the hessian ``\operatorname{Hess}F(x): T_x\mathcal M → T_x\mathcal M``, ``X ↦ \operatoname{Hess}F(x)[X] = ∇_ξ\operatorname{grad}f(x)``
* `x` – a point on the manifold ``x ∈ \mathcal M``
* `η` – an update tangential vector ``η ∈ T_x\mathcal M``
* `trust_region_radius` – a trust-region radius

# Optional

* `evaluation` – ([`AllocatingEvaluation`](@ref)) specify whether the gradient and hessian work by
   allocation (default) or [`MutatingEvaluation`](@ref) in place
* `preconditioner` – a preconditioner for the hessian H
* `θ` – 1+θ is the superlinear convergence target rate. The algorithm will
    terminate early if the residual was reduced by a power of 1+theta.
* `κ` – the linear convergence target rate: algorithm will terminate
    early if the residual was reduced by a factor of kappa.
* `randomize` – set to true if the trust-region solve is to be initiated with a
    random tangent vector. If set to true, no preconditioner will be
    used. This option is set to true in some scenarios to escape saddle
    points, but is otherwise seldom activated.
* `project_vector!` : (`copyto!`) specify a projection operation for tangent vectors
    for numerical stability. A function `(M, Y, p, X) -> ...` working in place of `Y`.
    per default, no projection is perfomed, set it to `project!` to activate projection.
* `stopping_criterion` – ([`StopWhenAny`](@ref), [`StopAfterIteration`](@ref),
    [`StopIfResidualIsReducedByFactor`](@ref), [`StopIfResidualIsReducedByPower`](@ref),
    [`StopWhenCurvatureIsNegative`](@ref), [`StopWhenTrustRegionIsExceeded`](@ref) )
    a functor inheriting from [`StoppingCriterion`](@ref) indicating when to stop,
    where for the default, the maximal number of iterations is set to the dimension of the
    manifold, the power factor is `θ`, the reduction factor is `κ`.
    .
* `return_options` – (`false`) – if activated, the extended result, i.e. the
    complete [`Options`](@ref) re returned. This can be used to access recorded values.
    If set to false (default) just the optimal value `x_opt` is returned

and the ones that are passed to [`decorate_options`](@ref) for decorators.

# Output

* `η` – an approximate solution of the trust-region subproblem in ``T_{x}\mathcal M``.
OR
* `options` - the options returned by the solver (see `return_options`)

# see also
[`trust_regions`](@ref)
"""
function truncated_conjugate_gradient_descent(
    M::AbstractManifold,
    F::TF,
    gradF::TG,
    x,
    η,
    H::TH,
    trust_region_radius::Float64;
    kwargs...,
) where {TF,TG,TH}
    x_res = allocate(x)
    copyto!(M, x_res, x)
    return truncated_conjugate_gradient_descent!(
        M, F, gradF, x_res, η, H, trust_region_radius; kwargs...
    )
end
@doc raw"""
    truncated_conjugate_gradient_descent!(M, F, gradF, x, η, HessF, trust_region_radius; kwargs...)

solve the trust-region subproblem in place of `x`.

# Input
# Input
* `M` – a manifold ``\mathcal M``
* `F` – a cost function ``F: \mathcal M → ℝ`` to minimize
* `gradF` – the gradient ``\operatorname{grad}F: \mathcal M → T\mathcal M`` of `F`
* `HessF` – the hessian ``\operatorname{Hess}F(x): T_x\mathcal M → T_x\mathcal M``, ``X ↦ \operatoname{Hess}F(x)[X] = ∇_ξ\operatorname{grad}f(x)``
* `x` – a point on the manifold ``x ∈ \mathcal M``
* `η` – an update tangential vector ``η ∈ T_x\mathcal M``
* `trust_region_radius` – a trust-region radius

For more details and all optional arguments, see [`truncated_conjugate_gradient_descent`](@ref).
"""
function truncated_conjugate_gradient_descent!(
    M::AbstractManifold,
    F::TF,
    gradF::TG,
    x,
    η,
    H::TH,
    trust_region_radius::Float64;
    evaluation=AllocatingEvaluation(),
    preconditioner::Tprec=(M, x, ξ) -> ξ,
    θ::Float64=1.0,
    κ::Float64=0.1,
    randomize::Bool=false,
    stopping_criterion::StoppingCriterion=StopWhenAny(
        StopAfterIteration(manifold_dimension(M)),
        StopWhenAll(StopIfResidualIsReducedByPower(θ), StopIfResidualIsReducedByFactor(κ)),
        StopWhenTrustRegionIsExceeded(),
        StopWhenCurvatureIsNegative(),
        StopWhenModelIncreased(),
    ),
    project_vector!::Proj=copyto!,
    return_options=false,
    kwargs..., #collect rest
) where {TF,TG,TH,Tprec,Proj}
    p = HessianProblem(M, F, gradF, H, preconditioner; evaluation=evaluation)
    o = TruncatedConjugateGradientOptions(
        p,
        x,
        η,
        trust_region_radius,
        randomize;
        θ=θ,
        κ=κ,
        stop=stopping_criterion,
        (project_vector!)=project_vector!,
    )
    o = decorate_options(o; kwargs...)
    resultO = solve(p, o)
    if return_options
        resultO
    else
        return get_solver_result(resultO)
    end
end
function initialize_solver!(p::HessianProblem, o::TruncatedConjugateGradientOptions)
    (o.randomize) || zero_vector!(p.M, o.η, o.x)
    o.Hη = o.randomize ? get_hessian(p, o.x, o.η) : zero_vector(p.M, o.x)
    o.gradient = get_gradient(p, o.x)
    o.residual = o.randomize ? o.gradient + o.Hη : o.gradient
    o.z = o.randomize ? o.residual : get_preconditioner(p, o.x, o.residual)
    o.δ = -deepcopy(o.z)
    o.Hδ = zero_vector(p.M, o.x)
    o.δHδ = inner(p.M, o.x, o.δ, o.Hδ)
    o.ηPδ = o.randomize ? inner(p.M, o.x, o.η, o.δ) : zero(o.δHδ)
    o.δPδ = inner(p.M, o.x, o.residual, o.z)
    o.ηPη = o.randomize ? inner(p.M, o.x, o.η, o.η) : zero(o.δHδ)
    if o.randomize
        o.model_value = inner(p.M, o.x, o.η, o.gradient) + 0.5 * inner(p.M, o.x, o.η, o.Hη)
    else
        o.model_value = 0
    end
    o.z_r = inner(p.M, o.x, o.z, o.residual)
    o.initialResidualNorm = sqrt(inner(p.M, o.x, o.residual, o.residual))
    return o
end
function step_solver!(
    p::P, o::O, ::Int
) where {P<:HessianProblem,O<:TruncatedConjugateGradientOptions}
    # Updates
    get_hessian!(p, o.Hδ, o.x, o.δ)
    o.δHδ = inner(p.M, o.x, o.δ, o.Hδ)
    α = o.z_r / o.δHδ
    ηPη_new = o.ηPη + 2 * α * o.ηPδ + α^2 * o.δPδ
    # Check against negative curvature and trust-region radius violation.
    if o.δHδ <= 0 || ηPη_new >= o.trust_region_radius^2
        τ = (-o.ηPδ + sqrt(o.ηPδ^2 + o.δPδ * (o.trust_region_radius^2 - o.ηPη))) / o.δPδ
        o.η = o.η + τ * o.δ
        o.Hη = o.Hη + τ * o.Hδ
        o.ηPη = ηPη_new
        return o
    end
    o.ηPη = ηPη_new
    new_η = o.η + α * o.δ
    new_Hη = o.Hη + α * o.Hδ
    # No negative curvature and o.η - α * (o.δ) inside TR: accept it.
    o.new_model_value =
        inner(p.M, o.x, new_η, o.gradient) + 0.5 * inner(p.M, o.x, new_η, new_Hη)
    o.new_model_value >= o.model_value && return o
    copyto!(p.M, o.η, o.x, new_η)
    o.model_value = o.new_model_value
    copyto!(p.M, o.Hη, o.x, new_Hη)
    o.residual = o.residual + α * o.Hδ
    # Precondition the residual.
    o.z = o.randomize ? o.residual : get_preconditioner(p, o.x, o.residual)
    zr = inner(p.M, o.x, o.z, o.residual)
    # Compute new search direction.
    β = zr / o.z_r
    o.z_r = zr
    o.δ = -o.z + β * o.δ
    o.project_vector!(p.M, o.δ, o.x, o.δ)
    o.ηPδ = β * (α * o.δPδ + o.ηPδ)
    o.δPδ = o.z_r + β^2 * o.δPδ
    return o
end
get_solver_result(o::O) where {O<:TruncatedConjugateGradientOptions} = o.η


text:
@Base.kwdef struct BaseAPI
	token::String = ""
	endpoint::String = ""
	proxy::Union{Nothing,String} = nothing
	# verify TLS certificate?
	verify::Bool = true
end

function send_request(api::BaseAPI, url::String, data::IO;)
	url = api.endpoint * url

	body = Dict(
		:file => data,
		:name => "file", # what does this do?
	) |> HTTP.Form

	HTTP.request("POST", url,
		verbose=2,
		require_ssl_verification = api.verify,
		proxy = api.proxy,
		headers = [
			"Authorization" => api.token,
			# "Content-Type" => "application/x-www-form-urlencoded"
		],
		body=body
	)
end

function send_request(api::BaseAPI, url::String;
			query=Dict(),
			body=HTTP.Form(Dict()),
			verb="GET",
		)

	url = api.endpoint * url

	HTTP.request(
		verb,
		url;
		body=body,
		query=query,
		verbose=2,
		require_ssl_verification = api.verify,
		proxy = api.proxy,
		headers = [
			"Authorization" => api.token,
			"Content-Type" => "application/x-www-form-urlencoded"
		],
	)
end

to_json(r::HTTP.Messages.Response) = r.body |> String |> JSON3.read

function parse_response(r::HTTP.Messages.Response, ::Type{T}) where T
	JSON3.read(r.body, T;
		parsequoted=true,
		dateformat=dateformat"yyyy-mm-dd HH:MM:ss"
	)
end

function struct_to_form_data(s::T) where T
	# HTTP.Form() seems to be stupid and can't handle DateTime and Nothing types.
	# we filter these out.
	dict = Dict(key => getfield(s, key) for key in fieldnames(T) if !isnothing(getfield(s, key)) && key ∉ [:lastchange])
	HTTP.Form(dict)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "x264Builder"
version = v"2019.05.25"

# Collection of sources required to build x264Builder
sources = [
    "https://download.videolan.org/pub/videolan/x264/snapshots/x264-snapshot-20190525-2245-stable.tar.bz2" =>
    "638581a18bff8e9375211955422eff145011c8ccfd0994d43bd194cd82984f7a",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd x264-snapshot-20190525-2245-stable/
if [[ "${target}" == x86_64* ]] || [[ "${target}" == i686* ]]; then
    apk add nasm
    export AS=nasm
else
    export AS="${CC}"
fi
./configure --prefix=$prefix --host=$target --enable-shared --disable-cli
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    # Windows
    Windows(:i686),
    Windows(:x86_64),

    # linux
    Linux(:i686, :glibc),
    Linux(:x86_64, :glibc),
    Linux(:aarch64, :glibc),
    Linux(:armv7l, :glibc),
    Linux(:powerpc64le, :glibc),

    # musl
    Linux(:i686, :musl),
    Linux(:x86_64, :musl),
    Linux(:aarch64, :musl),
    Linux(:armv7l, :musl),

    # The BSD's
    FreeBSD(:x86_64),
    MacOS(:x86_64),
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libx264", :libx264)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
# Test equivalence of functions and prox mappings by means of calculus rules

stuff = [
  Dict( "funcs"  => (IndBallLinf(), Conjugate(NormL1())),
        "args"   => ( randn(10), ),
        "gammas" => ( 1.0, )
      ),

  Dict( "funcs"  => (lambda -> (NormL1(lambda), Conjugate(IndBallLinf(lambda))))(0.1+10.0*rand()),
        "args"   => ( 5.0*sign.(randn(10)) + 5.0*randn(10),
                      5.0*sign.(randn(20)) + 5.0*randn(20) ),
        "gammas" => ( 0.5+rand(), 0.5+rand() )
      ),

  Dict( "funcs"  => (lambda -> (IndBallLinf(lambda), Conjugate(NormL1(lambda))))(0.1+10.0*rand()),
        "args"   => ( 5.0*sign.(randn(10)) + 5.0*randn(10),
                      5.0*sign.(randn(20)) + 5.0*randn(20) ),
        "gammas" => ( 0.5+rand(), 0.5+rand() )
      ),

  Dict( "funcs"  => (lambda -> (NormL1(lambda), Conjugate(IndBox(-lambda,lambda))))(0.1+10.0*rand(30)),
        "args"   => ( 5.0*sign.(randn(30)) + 5.0*randn(30), ),
        "gammas" => ( 0.5+rand(), 0.5+rand() )
      ),

  Dict( "funcs"  => (lambda -> (IndBox(-lambda,lambda), Conjugate(NormL1(lambda))))(0.1+10.0*rand(30)),
        "args"   => ( 5.0*sign.(randn(30)) + 5.0*randn(30), ),
        "gammas" => ( 0.5+rand(), 0.5+rand() )
      ),

  Dict( "funcs"  => (lambda -> (NormL2(lambda), Conjugate(IndBallL2(lambda))))(0.1+10.0*rand()),
        "args"   => ( 5.0*sign.(randn(10)) + 5.0*randn(10),
                      5.0*sign.(randn(20)) + 5.0*randn(20) ),
        "gammas" => ( 0.5+rand(), 0.5+rand() )
      ),

  Dict( "funcs"  => (lambda -> (IndBallL2(lambda), Conjugate(NormL2(lambda))))(0.1+10.0*rand()),
        "args"   => ( 5.0*sign.(randn(10)) + 5.0*randn(10),
                      5.0*sign.(randn(20)) + 5.0*randn(20) ),
        "gammas" => ( 0.5+rand(), 0.5+rand() )
      ),

  Dict( "funcs"  => ((a, b, mu) -> (LogBarrier(a, b, mu), Postcompose(Precompose(LogBarrier(), a, b), mu)))(2.0, 0.5, 1.0),
        "args"   => ( rand(10), rand(10) ),
        "gammas" => ( 0.5+rand(), 0.5+rand() )
      ),

  Dict( "funcs"  => (p -> (IndPoint(p), IndBox(p, p)))(randn(50)),
        "args"   => ( randn(50), randn(50), randn(50) ),
        "gammas" => ( 1.0, rand(), 5.0*rand() )
      ),

  Dict( "funcs"  => (IndZero(), IndBox(0, 0)),
        "args"   => ( randn(50), randn(50), randn(50) ),
        "gammas" => ( 1.0, rand(), 5.0*rand() )
      ),

  Dict( "funcs"  => (IndFree(), IndBox(-Inf, +Inf)) ,
        "args"   => ( randn(50), randn(50), randn(50) ),
        "gammas" => ( 1.0, rand(), 5.0*rand() )
      ),

  Dict( "funcs"  => (IndNonnegative(), IndBox(0.0, Inf)),
        "args"   => ( randn(50), randn(50), randn(50) ),
        "gammas" => ( 1.0, rand(), 5.0*rand() )
      ),

  Dict( "funcs"  => (IndNonpositive(), IndBox(-Inf, 0.0)),
        "args"   => ( randn(50), randn(50), randn(50) ),
        "gammas" => ( 1.0, rand(), 5.0*rand() )
      ),

  Dict( "funcs"  => (lambda -> (SqrNormL2(lambda), Conjugate(SqrNormL2(1.0/lambda))))(0.1+5.0*rand()),
        "args"   => ( randn(50), randn(50), randn(50) ),
        "gammas" => ( 1.0, rand(), 5.0*rand() )
      ),

  Dict( "funcs"  => ((A, b) -> (LeastSquares(A, b), Tilt(LeastSquares(A, zeros(size(A, 1))), -A'*b, 0.5*vecdot(b, b))))(randn(10,20), randn(10)),
        "args"   => ( randn(20), randn(20), randn(20) ),
        "gammas" => ( 1.0, rand(), 5.0*rand() )
      ),

  Dict( "funcs"  => ((lambda, rho) -> (ElasticNet(lambda,rho), Regularize(NormL1(lambda),rho)))(rand(), rand()),
        "args"   => ( randn(20), randn(20), randn(20) ),
        "gammas" => ( 1.0, rand(), 5.0*rand() )
      ),

  Dict( "funcs"  => ((b, mu) -> (HingeLoss(b, mu), Postcompose(Precompose(SumPositive(), -b, 1.0), mu)))([0.5+rand(10); -0.5-rand(10)], 0.5+rand()),
        "args"   => ( randn(20), randn(20), randn(20) ),
        "gammas" => ( 1.0, rand(), 5.0*rand() )
      )
]

for i = 1:length(stuff)
  println("----------------------------------------------------------")
  f = stuff[i]["funcs"][1]
  g = stuff[i]["funcs"][2]
  println(string(rpad(typeof(f), 25, " "), " VS ", typeof(g)))

  for j = 1:length(stuff[i]["args"])
    x = stuff[i]["args"][j]
    gamma = stuff[i]["gammas"][j]

    # compare the three versions (for f)
    yf, fy = prox_test(f, x, gamma)

    # compare the three versions (for g)
    yg, gy = prox_test(g, x, gamma)

    # compare results of f and g
    @test vecnorm(yf-yg, Inf)/(1+norm(yf, Inf)) <= TOL_ASSERT
    @test fy == gy || abs(fy-gy)/(1+abs(fy)) <= TOL_ASSERT
  end

end


text:
@generated function apodization!(p::NFFTPlan{D,DIM,T}, f::AbstractArray{U,D},
                                 g::StridedArray{Complex{T},D}) where {D,DIM,T,U}
    quote
        offset = round(Int, p.n[$DIM] - p.N[$DIM]/2) - 1

        @nloops $D l f d->begin
            # preexpr
            if d == $DIM
                gidx_d = rem(l_d+offset, p.n[d]) + 1
                winidx = l_d
            else
                gidx_d = l_d
            end
        end begin
            # bodyexpr
            (@nref $D g gidx) = (@nref $D f l) * p.windowHatInvLUT[1][winidx]
        end
    end
end

@generated function apodization_adjoint!(p::NFFTPlan{D,DIM,T},
           g::AbstractArray{Complex{T},D}, f::StridedArray{U,D}) where {D,DIM,T,U}
    quote
        offset = round(Int, p.n[$DIM] - p.N[$DIM]/2) - 1

        @nloops $D l f d->begin
            # preexpr
            if d == $DIM
                gidx_d = rem(l_d+offset, p.n[d]) + 1
                winidx = l_d
            else
                gidx_d = l_d
            end
        end begin
            # bodyexpr
            (@nref $D f l) = (@nref $D g gidx) * p.windowHatInvLUT[1][winidx]
        end
    end
end

@generated function convolve!(p::NFFTPlan{D,DIM,T}, g::AbstractArray{Complex{T},D},
                              fHat::StridedArray{U,D}) where {D,DIM,T,U}
    quote
        fill!(fHat, zero(T))
        scale = 1.0 / p.m * (p.K-1)

        for k in 1:p.M
            xscale = p.x[k] * p.n[$DIM]
            c = floor(Int, xscale)
            @nloops $D l d->begin
                # rangeexpr
                if d == $DIM
                    (c-p.m):(c+p.m)
                else
                    1:size(g,d)
                end
            end d->begin
                # preexpr
                if d == $DIM
                    gidx_d = rem(l_d+p.n[d], p.n[d]) + 1
                    idx = abs( (xscale - l_d)*scale ) + 1
                    idxL = floor(idx)
                    idxInt = Int(idxL)
                    tmpWin = p.windowLUT[1][idxInt] + ( idx-idxL ) * (p.windowLUT[1][idxInt+1] - p.windowLUT[1][idxInt])
                    fidx_d = k
                else
                    gidx_d = l_d
                    fidx_d = l_d
                end
            end begin
                # bodyexpr
                (@nref $D fHat fidx) += (@nref $D g gidx) * tmpWin
            end
        end
    end
end

@generated function convolve_adjoint!(p::NFFTPlan{D,DIM,T}, fHat::AbstractArray{U,D},
                                      g::StridedArray{Complex{T},D}) where {D,DIM,T,U}
    quote
        fill!(g, zero(T))
        scale = 1.0 / p.m * (p.K-1)

        for k in 1:p.M
            xscale = p.x[k] * p.n[$DIM]
            c = floor(Int, xscale)
            @nloops $D l d->begin
                # rangeexpr
                if d == $DIM
                    (c-p.m):(c+p.m)
                else
                    1:size(g,d)
                end
            end d->begin
                # preexpr
                if d == $DIM
                    gidx_d = rem(l_d+p.n[d], p.n[d]) + 1
                    idx = abs( (xscale - l_d)*scale ) + 1
                    idxL = floor(idx)
                    idxInt = Int(idxL)
                    tmpWin = p.windowLUT[1][idxInt] + ( idx-idxL ) * (p.windowLUT[1][idxInt+1] - p.windowLUT[1][idxInt])
                    fidx_d = k
                else
                    gidx_d = l_d
                    fidx_d = l_d
                end
            end begin
                # bodyexpr
                (@nref $D g gidx) += (@nref $D fHat fidx) * tmpWin
            end
        end
    end
end


text:
module HermiteDualContouring

using RegionTrees
using ConstructiveSolidGeometry: Surface, distance_field, normal_field, unitize
using Roots
using LinearAlgebra

import RegionTrees: needs_refinement, refine_data

export HermiteDualContour

struct DualContourRefinery{S<:Surface} <: AbstractRefinery
    surface_def::S
    atol::Float64
    rtol::Float64
    surfcellmax::Float64
end

function find_surface_intersect(v1, v2, contour::T) where {T<:Surface}
    edge_func = (α)->v1+α*(v2-v1)
    line_func = (α)->distance_field(contour, edge_func(α))
    #Because of how edge_func is parameterized, the range is always [0,1]
    α_final = find_zero(line_func,(0,1))
    return edge_func(α_final)
end

function needs_refinement(refinery::DualContourRefinery, cell::Cell)
    minimum(cell.boundary.widths) > refinery.atol || return false

    if !isempty(cell.data.residual) && !isnan(cell.data.residual)
        if minimum(cell.boundary.widths)>refinery.surfcellmax
            return true
        elseif !isapprox(cell.data.residual, 0.0, rtol=refinery.rtol, atol=refinery.atol)
            return true
        else
            return false
        end
    else # Don't check for zero or odd normals, just break up the cell if the QEF fails. Is this robust/efficient?
        if isempty(cell.data.residual)
            # check to see if the distances from the corners could indicate a surface within the cell
            max_width = maximum(cell.boundary.widths)
            for vertex in vertices(cell)
                if abs(distance_field(refinery.surface_def,vertex)) <= max_width/sqrt(2)
                    return true #still creates a lot of fine granularity in empty space, but much better
                end
            end
            return false #line only reached if no vertices meet distance condition
        else
            if minimum(cell.boundary.widths)>refinery.surfcellmax
                return true
            else
                return false #isnan case
            end
        end
    end
end

function refine_data(refinery::DualContourRefinery, cell::Cell, indices)
    refine_data(refinery, child_boundary(cell, indices))
end

struct CellData
    edges
    p
    n
    qef_min
    residual
    function CellData(
        edges=[],
        p=[],
        n=[],
        qef_min=[],
        residual=[])
        new(edges,p,n,qef_min,residual)
    end
end

function isOutside(vert_array, c, rtol, atol)
    cat_verts = hcat(collect(vert_array)...)
    maxX = maximum(cat_verts[1,:])
    minX = minimum(cat_verts[1,:])

    maxY = maximum(cat_verts[2,:])
    minY = minimum(cat_verts[2,:])

    if !isapprox(maxX,c[1], rtol=rtol, atol=atol) && maxX<c[1]
        return true
    elseif !isapprox(minX,c[1], rtol=rtol, atol=atol) && minX>c[1]
        return true
    elseif !isapprox(maxY,c[2], rtol=rtol, atol=atol) && maxY<c[2]
        return true
    elseif !isapprox(minY,c[2], rtol=rtol, atol=atol) && minY>c[2]
        return true
    else
        return false
    end
end

function refine_data(refinery::DualContourRefinery, boundary::HyperRectangle)
    p=[]
    n=[]
    edges=[]
    vert_array=vertices(boundary) #ideally this would be examining edges not vertices
    v_order = [1,2,4,3,1] #2D for now, edges would not require the list

    #logic to check if cell is entirely in/outside of surface
    val_array = map((x)->distance_field(refinery.surface_def, x), vert_array)#for now, 2D
    sign_set = Set(sign.(val_array))
    is_on_plane = length(val_array)-count(!iszero,val_array)>1 && -1.0 in sign_set
    is_not_on_surface = length(sign_set)==1 && !(0.0 in sign_set || -0.0 in sign_set) #what to do if all verts are on surface?
    has_single_corner = (length(sign_set)<3 && (0.0 in sign_set || -0.0 in sign_set)) #because floats are terrible
    if  !is_on_plane && (is_not_on_surface || has_single_corner)
        #check to see if the cell contains a surface or part of one within it
        return CellData()
    end

    for v_idx = 1:length(vert_array)
        current_vert = vert_array[v_order[v_idx]]
        next_vert = vert_array[v_order[v_idx+1]]

        current_val = val_array[v_order[v_idx]]
        next_val = val_array[v_order[v_idx+1]]

        current_sign = sign(current_val)
        next_sign = sign(next_val)

        sign_change = current_sign!=next_sign
        current_sign_zero = current_sign==0.0 || current_sign==-0.0
        next_sign_zero = next_sign==0.0 || next_sign == -0.0

        if sign_change || (current_sign_zero && next_sign_zero)
            if current_val==zero(current_val)
                p_local = current_vert
            elseif next_val==zero(next_val)
                continue
            else
                p_local = find_surface_intersect(current_vert, next_vert, refinery.surface_def)
            end
            n_local = normal_field(refinery.surface_def, p_local) #for now, 2D
            push!(p,p_local)
            push!(n,n_local) #for now, 2D
            push!(edges, [current_vert, next_vert])
        end
    end

    A = hcat(n...)'
    simA=similar(A,size(A)...)
    map!(x->x,simA,A) #convert out of SMatrix so that A\b works correctly
    b = dot.(p,n)
    qef_min=0#remember how to initialize as undefined
    residual=0#remember how to initialize as undefined
    if length(p)==1 #may be useful for debugging where flat surfaces are not involved
        qef_min = Inf
        residual = Inf
    else
        try
            qef_min = simA\b
            residual = distance_field(refinery.surface_def, qef_min)
        catch
            qef_min = NaN
            residual = NaN
        end
        if !isnan(residual)
            #occurs when a cell corner aligns with a surface corner, needs to be adjusted for dimensionality
            corner_idx = findall(i->isapprox(i,qef_min), p)
            if !isempty(corner_idx) && length(p)>2
                deleteat!(edges, corner_idx)
                deleteat!(p, corner_idx)
                deleteat!(n, corner_idx)
            end
            if isOutside(vert_array,qef_min, refinery.rtol, refinery.atol)
                qef_min = Inf
                residual = Inf
            end
        end
    end

    return CellData(
        edges,
        p,
        n,
        qef_min,
        residual)
end

struct HermiteDualContour{C<:Cell} <: Surface
    root::C
end

function HermiteDualContour(surface_def::Surface, origin::AbstractArray,
        widths::AbstractArray,
        rtol = 1e-2,
        atol = 1e-2,
        surfcellmax = 1e-1)
    refinery = DualContourRefinery(surface_def, atol, rtol, surfcellmax)
    boundary = HyperRectangle(origin, widths)
    root = Cell(boundary, refine_data(refinery, boundary))
    adaptivesampling!(root, refinery)
    HermiteDualContour(root)
end

# end

end # module

max_stars_repo_name: CrashBurnRepeat/HermiteDualContouring

max_issues_repo_name: CrashBurnRepeat/HermiteDualContouring

max_forks_repo_name: CrashBurnRepeat/HermiteDualContouring


text:
#emul TIMES
function frule(
    (_, ΔA, ΔB, _),
    ::typeof(emul),
    A::GBArray,
    B::GBArray,
    ::typeof(BinaryOps.TIMES)
)
    Ω = emul(A, B, BinaryOps.TIMES)
    ∂Ω = emul(unthunk(ΔA), B, BinaryOps.TIMES) + emul(unthunk(ΔB), A, BinaryOps.TIMES)
    return Ω, ∂Ω
end
function frule((_, ΔA, ΔB), ::typeof(emul), A::GBArray, B::GBArray)
    return frule((nothing, ΔA, ΔB, nothing), emul, A, B, BinaryOps.TIMES)
end

function rrule(::typeof(emul), A::GBArray, B::GBArray, ::typeof(BinaryOps.TIMES))
    function timespullback(ΔΩ)
        ∂A = emul(unthunk(ΔΩ), B)
        ∂B = emul(unthunk(ΔΩ), A)
        return NoTangent(), ∂A, ∂B, NoTangent()
    end
    return emul(A, B, BinaryOps.TIMES), timespullback
end

function rrule(::typeof(emul), A::GBArray, B::GBArray)
    Ω, fullpb = rrule(emul, A, B, BinaryOps.TIMES)
    emulpb(ΔΩ) = fullpb(ΔΩ)[1:3]
    return Ω, emulpb
end

############
# eadd rules
############

# PLUS
######

function frule(
    (_, ΔA, ΔB, _),
    ::typeof(eadd),
    A::GBArray,
    B::GBArray,
    ::typeof(BinaryOps.PLUS)
)
    Ω = eadd(A, B, BinaryOps.PLUS)
    ∂Ω = eadd(unthunk(ΔA), unthunk(ΔB), BinaryOps.PLUS)
    return Ω, ∂Ω
end
function frule((_, ΔA, ΔB), ::typeof(eadd), A::GBArray, B::GBArray)
    return frule((nothing, ΔA, ΔB, nothing), eadd, A, B, BinaryOps.PLUS)
end

function rrule(::typeof(eadd), A::GBArray, B::GBArray, ::typeof(BinaryOps.PLUS))
    function pluspullback(ΔΩ)
        return (
            NoTangent(),
            mask(unthunk(ΔΩ), A; structural = true),
            mask(unthunk(ΔΩ), B; structural = true),
            NoTangent()
        )
    end
    return eadd(A, B, BinaryOps.PLUS), pluspullback
end

function rrule(::typeof(eadd), A::GBArray, B::GBArray)
    Ω, fullpb = rrule(eadd, A, B, BinaryOps.PLUS)
    eaddpb(ΔΩ) = fullpb(ΔΩ)[1:3]
    return Ω, eaddpb
end


text:
# ---
# title: 893. Groups of Special-Equivalent Strings
# id: problem893
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/groups-of-special-equivalent-strings/description/>
# hidden: true
# ---
#
# You are given an array `A` of strings.
#
# A _move  onto `S`_ consists of swapping any two even indexed characters of
# `S`, or any two odd indexed characters of `S`.
#
# Two strings `S` and `T` are  _special-equivalent_  if after any number of
# _moves onto`S`_, `S == T`.
#
# For example, `S = "zzxy"` and `T = "xyzz"` are special-equivalent because we
# may make the moves `"zzxy" -> "xzzy" -> "xyzz"` that swap `S[0]` and `S[2]`,
# then `S[1]` and `S[3]`.
#
# Now, a _group of special-equivalent strings from`A`_ is a non-empty subset of
# A such that:
#
#   1. Every pair of strings in the group are special equivalent, and;
#   2. The group is the largest size possible (ie., there isn't a string S not in the group such that S is special equivalent to every string in the group)
#
# Return the number of groups of special-equivalent strings from `A`.
#
#
#
# **Example 1:**
#
#
#
#     Input: ["abcd","cdab","cbad","xyzz","zzxy","zzyx"]
#     Output: 3
#     Explanation:
#     One group is ["abcd", "cdab", "cbad"], since they are all pairwise special equivalent, and none of the other strings are all pairwise special equivalent to these.
#
#     The other two groups are ["xyzz", "zzxy"] and ["zzyx"].  Note that in particular, "zzxy" is not special equivalent to "zzyx".
#
#
# **Example 2:**
#
#
#
#     Input: ["abc","acb","bac","bca","cab","cba"]
#     Output: 3
#
#
#
# **Note:**
#
#   * `1 <= A.length <= 1000`
#   * `1 <= A[i].length <= 20`
#   * All `A[i]` have the same length.
#   * All `A[i]` consist of only lowercase letters.
#
#
## @lc code=start
using LeetCode

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


text:
module DAslabTests

using Test
using ForwardDiff
using LightPropagation: fluence_DA_semiinf_CW, fluence_DA_semiinf_TD
using LightPropagation: fluence_DA_slab_CW, fluence_DA_slab_TD
using LightPropagation: flux_DA_slab_CW, flux_DA_slab_TD, flux_DA_semiinf_TD, flux_DA_semiinf_CW

# Test slabs against semi-inf for sufficiently thick slabs

# CW
@test fluence_DA_slab_CW(1.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, s = 20.0, z = 0.0) == 0.024035998288332954
@test fluence_DA_slab_CW(2.4, 0.6, 18.0, n_ext = 1.3, n_med = 1.2, s = 15.0, z = 1.0) == 1.77784095787581e-7
@test fluence_DA_semiinf_CW(1.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, z = 0.0) ≈ fluence_DA_slab_CW(1.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, z = 0.0, s = 20.0)
@test fluence_DA_semiinf_CW(4.0, 0.01, 50.0, n_ext = 1.0, n_med = 1.0, z = 0.0) ≈ fluence_DA_slab_CW(4.0, 0.01, 50.0,n_ext = 1.0, n_med = 1.0, z = 0.0, s = 20.0)

# TD
t = 0.1:0.5:8
@test fluence_DA_semiinf_TD(t, 1.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, z = 0.0) ≈ fluence_DA_slab_TD(t, 1.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, z = 0.0, s = 40.0)
@test fluence_DA_semiinf_TD(t, 1.0, 0.01, 70.0, n_ext = 1.0, n_med = 1.0, z = 0.0) ≈ fluence_DA_slab_TD(t, 1.0, 0.01, 70.0, n_ext = 1.0, n_med = 1.0, z = 0.0, s = 40.0)

@test fluence_DA_semiinf_TD(t, 3.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, z = 0.0) ≈ fluence_DA_slab_TD(t, 3.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, z = 0.0, s = 40.0)
@test fluence_DA_semiinf_TD(t, 3.0, 0.01, 70.0, n_ext = 1.0, n_med = 1.0, z = 0.0) ≈ fluence_DA_slab_TD(t, 3.0, 0.01, 70.0, n_ext = 1.0, n_med = 1.0, z = 0.0, s = 40.0)
@test fluence_DA_semiinf_TD(t, 3.0, 0.01, 70.0, n_ext = 1.4, n_med = 1.2, z = 0.1) ≈ fluence_DA_slab_TD(t, 3.0, 0.01, 70.0, n_ext = 1.4, n_med = 1.2, z = 0.1, s = 40.0)

## check flux

# reflectance

D = 1 / (3 * 10.0)
@test (D * ForwardDiff.derivative(dz -> fluence_DA_slab_TD(1.0, 1.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, s = 1.0, z = dz), 0.0)) ≈ flux_DA_slab_TD(1.0, 1.0, 0.1, 10.0,  n_ext = 1.0, n_med = 1.0, s = 1.0, z = 0.0)

D = 1 / (3 * 20.0)
@test (D * ForwardDiff.derivative(dz -> fluence_DA_slab_TD(1.0, 1.3, 0.3, 20.0, n_ext = 1.0, n_med = 1.0, s = 1.0, z = dz), 0.0)) ≈ flux_DA_slab_TD(1.0, 1.3, 0.3, 20.0,  n_ext = 1.0, n_med = 1.0, s = 1.0, z = 0.0)

# transmittance
D = 1 / (3 * 20.0)
@test (-D * ForwardDiff.derivative(dz -> fluence_DA_slab_TD(1.0, 1.3, 0.3, 20.0, n_ext = 1.0, n_med = 1.0, s = 1.0, z = dz), 1.0)) ≈ flux_DA_slab_TD(1.0, 1.3, 0.3, 20.0,  n_ext = 1.0, n_med = 1.0, s = 1.0, z = 1.0)

## CW
@test flux_DA_slab_CW(1.0, 0.1, 10.0; n_ext = 1.0, n_med = 1.0, s = 40.0, z = 0.0, xs = 15) ≈ flux_DA_semiinf_CW(1.0, 0.1, 10.0; n_ext = 1.0, n_med = 1.0)

end

max_issues_repo_name: heltonmc/LightPropagation


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule FastJet_Julia_Wrapper_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("FastJet_Julia_Wrapper")
JLLWrappers.@generate_main_file("FastJet_Julia_Wrapper", UUID("33ff6b16-4e4c-5b96-9396-f34816859298"))
end  # module FastJet_Julia_Wrapper_jll


text:
@testset "dijkstra" begin

    @testset "small_network" begin
        file = "data/networks/distance/weighted_elist_small.txt"
        graph = read_graph(file, "distance")
        distances, pred = dijkstra(graph, 1)
        @test distances[5] == 36.0
        @test distances[10] == 19.0
        dist_5, path_5 = dijkstra(graph, 1, 5)
        @test dist_5 == 36.0
        @test path_5 == [1, 12, 5]
        dist_10, path_10 = dijkstra(graph, 1, 10)
        @test dist_10 == 19.0
        @test path_10 == [1, 10]
    end

    @testset "large_network" begin
        file = "data/networks/distance/weighted_elist_large.txt"
        graph = read_graph(file, "distance")
        distances, pred = dijkstra(graph, 1)
        @test distances[130] == 49.0
        @test distances[180] == 25.0
        @test distances[200] == 46.0
        dist_130, path_130 = dijkstra(graph, 1, 130)
        @test dist_130 == 49.0
        @test path_130 == [1, 96, 30, 130]
        dist_200, path_200 = dijkstra(graph, 1, 200)
        @test dist_200 == 46.0
        @test path_200 == [1, 96, 140, 14, 200]
    end

    @testset "vlarge_network" begin
        file = "data/networks/distance/weighted_elist_vlarge.txt"
        graph = read_graph(file, "distance")
        distances, pred = dijkstra(graph, 1)
        @test distances[250] == 19.0
        @test distances[1000] == 20.0
        @test distances[1500] == 16.0
        dist_250, path_250 = dijkstra(graph, 1, 250)
        @test dist_250 == 19.0
        @test path_250 == [1, 293, 680, 139, 1439, 1109, 250]
        dist_1500, path_1500 = dijkstra(graph, 1, 1500)
        @test dist_1500 == 16.0
        @test path_1500 == [1, 293, 680, 139, 1371, 980, 1500]
    end

end


text:
"""
An abstract type for which all concrete types expose a database-like
interface.

**Common methods**

An AbstractDataFrame is a two-dimensional table with Symbols for
column names. An AbstractDataFrame is also similar to an Associative
type in that it allows indexing by a key (the columns).

The following are normally implemented for AbstractDataFrames:

* [`describe`](@ref) : summarize columns
* [`dump`](@ref) : show structure
* `hcat` : horizontal concatenation
* `vcat` : vertical concatenation
* `names` : columns names
* [`names!`](@ref) : set columns names
* [`rename!`](@ref) : rename columns names based on keyword arguments
* [`eltypes`](@ref) : `eltype` of each column
* `length` : number of columns
* `size` : (nrows, ncols)
* [`head`](@ref) : first `n` rows
* [`tail`](@ref) : last `n` rows
* `convert` : convert to an array
* [`completecases`](@ref) : boolean vector of complete cases (rows with no missings)
* [`dropmissing`](@ref) : remove rows with missing values
* [`dropmissing!`](@ref) : remove rows with missing values in-place
* [`nonunique`](@ref) : indexes of duplicate rows
* [`unique!`](@ref) : remove duplicate rows
* `similar` : a DataFrame with similar columns as `d`

**Indexing**

Table columns are accessed (`getindex`) by a single index that can be
a symbol identifier, an integer, or a vector of each. If a single
column is selected, just the column object is returned. If multiple
columns are selected, some AbstractDataFrame is returned.

```julia
d[:colA]
d[3]
d[[:colA, :colB]]
d[[1:3; 5]]
```

Rows and columns can be indexed like a `Matrix` with the added feature
of indexing columns by name.

```julia
d[1:3, :colA]
d[3,3]
d[3,:]
d[3,[:colA, :colB]]
d[:, [:colA, :colB]]
d[[1:3; 5], :]
```

`setindex` works similarly.
"""
abstract type AbstractDataFrame end

##############################################################################
##
## Interface (not final)
##
##############################################################################

# index(df) => AbstractIndex
# nrow(df) => Int
# ncol(df) => Int
# getindex(...)
# setindex!(...) exclusive of methods that add new columns

##############################################################################
##
## Basic properties of a DataFrame
##
##############################################################################

struct Cols{T <: AbstractDataFrame} <: AbstractVector{Any}
    df::T
end
Base.start(::Cols) = 1
Base.done(itr::Cols, st) = st > length(itr.df)
Base.next(itr::Cols, st) = (itr.df[st], st + 1)
Base.length(itr::Cols) = length(itr.df)
Base.size(itr::Cols, ix) = ix==1 ? length(itr) : throw(ArgumentError("Incorrect dimension"))
Base.size(itr::Cols) = (length(itr.df),)
Base.IndexStyle(::Type{<:Cols}) = IndexLinear()
Base.getindex(itr::Cols, inds...) = getindex(itr.df, inds...)

# N.B. where stored as a vector, 'columns(x) = x.vector' is a bit cheaper
columns(df::T) where {T <: AbstractDataFrame} = Cols{T}(df)

Base.names(df::AbstractDataFrame) = names(index(df))
_names(df::AbstractDataFrame) = _names(index(df))

# TODO: remove allow_duplicates after deprecation period
"""
Set column names

```julia
names!(df::AbstractDataFrame, vals)
```

**Arguments**

* `df` : the AbstractDataFrame
* `vals` : column names, normally a Vector{Symbol} the same length as
  the number of columns in `df`
* `makeunique` : if `false` (the default), an error will be raised
  if duplicate names are found; if `true`, duplicate names will be suffixed
  with `_i` (`i` starting at 1 for the first duplicate).

**Result**

* `::AbstractDataFrame` : the updated result

**Examples**

```julia
df = DataFrame(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
names!(df, [:a, :b, :c])
names!(df, [:a, :b, :a])  # throws ArgumentError
names!(df, [:a, :b, :a], makeunique=true)  # renames second :a to :a_1
```

"""
function names!(df::AbstractDataFrame, vals; allow_duplicates=false, makeunique::Bool=false)
    if allow_duplicates
        Base.depwarn("Keyword argument allow_duplicates is deprecated. Use makeunique.", :names!)
    end
    names!(index(df), vals, allow_duplicates=allow_duplicates, makeunique=makeunique)
    return df
end

function rename!(df::AbstractDataFrame, args...)
    rename!(index(df), args...)
    return df
end
function rename!(f::Function, df::AbstractDataFrame)
    rename!(f, index(df))
    return df
end

rename(df::AbstractDataFrame, args...) = rename!(copy(df), args...)
rename(f::Function, df::AbstractDataFrame) = rename!(f, copy(df))

"""
Rename columns

```julia
rename!(df::AbstractDataFrame, (from => to)::Pair{Symbol, Symbol}...)
rename!(df::AbstractDataFrame, d::AbstractDict{Symbol,Symbol})
rename!(df::AbstractDataFrame, d::AbstractArray{Pair{Symbol,Symbol}})
rename!(f::Function, df::AbstractDataFrame)
rename(df::AbstractDataFrame, (from => to)::Pair{Symbol, Symbol}...)
rename(df::AbstractDataFrame, d::AbstractDict{Symbol,Symbol})
rename(df::AbstractDataFrame, d::AbstractArray{Pair{Symbol,Symbol}})
rename(f::Function, df::AbstractDataFrame)
```

**Arguments**

* `df` : the AbstractDataFrame
* `d` : an Associative type or an AbstractArray of pairs that maps
  the original names to new names
* `f` : a function which for each column takes the old name (a Symbol)
  and returns the new name (a Symbol)

**Result**

* `::AbstractDataFrame` : the updated result

New names are processed sequentially. A new name must not already exist in the `DataFrame`
at the moment an attempt to rename a column is performed.

**Examples**

```julia
df = DataFrame(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
rename(df, :i => :A, :x => :X)
rename(df, [:i => :A, :x => :X])
rename(df, Dict(:i => :A, :x => :X))
rename(x -> Symbol(uppercase(string(x))), df)
rename!(df, Dict(:i =>: A, :x => :X))
```

"""
(rename!, rename)

"""
Return element types of columns

```julia
eltypes(df::AbstractDataFrame)
```

**Arguments**

* `df` : the AbstractDataFrame

**Result**

* `::Vector{Type}` : the element type of each column

**Examples**

```julia
df = DataFrame(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
eltypes(df)
```

"""
eltypes(df::AbstractDataFrame) = map!(eltype, Vector{Type}(undef, size(df,2)), columns(df))

Base.size(df::AbstractDataFrame) = (nrow(df), ncol(df))
function Base.size(df::AbstractDataFrame, i::Integer)
    if i == 1
        nrow(df)
    elseif i == 2
        ncol(df)
    else
        throw(ArgumentError("DataFrames only have two dimensions"))
    end
end

Base.length(df::AbstractDataFrame) = ncol(df)
Compat.lastindex(df::AbstractDataFrame) = ncol(df)
Compat.lastindex(df::AbstractDataFrame, i) = last(axes(df, i))
Compat.axes(df, i) = axes(df)[i]

Base.ndims(::AbstractDataFrame) = 2

##############################################################################
##
## Similar
##
##############################################################################

"""
    similar(df::DataFrame[, rows::Integer])

Create a new `DataFrame` with the same column names and column element types
as `df`. An optional second argument can be provided to request a number of rows
that is different than the number of rows present in `df`.
"""
function Base.similar(df::AbstractDataFrame, rows::Integer = size(df, 1))
    rows < 0 && throw(ArgumentError("the number of rows must be positive"))
    DataFrame(Any[similar(x, rows) for x in columns(df)], copy(index(df)))
end

##############################################################################
##
## Equality
##
##############################################################################

function Base.:(==)(df1::AbstractDataFrame, df2::AbstractDataFrame)
    size(df1, 2) == size(df2, 2) || return false
    isequal(index(df1), index(df2)) || return false
    eq = true
    for idx in 1:size(df1, 2)
        coleq = df1[idx] == df2[idx]
        # coleq could be missing
        !isequal(coleq, false) || return false
        eq &= coleq
    end
    return eq
end

function Base.isequal(df1::AbstractDataFrame, df2::AbstractDataFrame)
    size(df1, 2) == size(df2, 2) || return false
    isequal(index(df1), index(df2)) || return false
    for idx in 1:size(df1, 2)
        isequal(df1[idx], df2[idx]) || return false
    end
    return true
end

##############################################################################
##
## Associative methods
##
##############################################################################

Base.haskey(df::AbstractDataFrame, key::Any) = haskey(index(df), key)
Base.get(df::AbstractDataFrame, key::Any, default::Any) = haskey(df, key) ? df[key] : default
Base.isempty(df::AbstractDataFrame) = size(df, 1) == 0 || size(df, 2) == 0

##############################################################################
##
## Description
##
##############################################################################

head(df::AbstractDataFrame, r::Int) = df[1:min(r,nrow(df)), :]
head(df::AbstractDataFrame) = head(df, 6)
tail(df::AbstractDataFrame, r::Int) = df[max(1,nrow(df)-r+1):nrow(df), :]
tail(df::AbstractDataFrame) = tail(df, 6)

"""
Show the first or last part of an AbstractDataFrame

```julia
head(df::AbstractDataFrame, r::Int = 6)
tail(df::AbstractDataFrame, r::Int = 6)
```

**Arguments**

* `df` : the AbstractDataFrame
* `r` : the number of rows to show

**Result**

* `::AbstractDataFrame` : the first or last part of `df`

**Examples**

```julia
df = DataFrame(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
head(df)
tail(df)
```

"""
(head, tail)

# get the structure of a df
"""
Show the structure of an AbstractDataFrame, in a tree-like format

```julia
dump(df::AbstractDataFrame, n::Int = 5)
dump(io::IO, df::AbstractDataFrame, n::Int = 5)
```

**Arguments**

* `df` : the AbstractDataFrame
* `n` : the number of levels to show
* `io` : optional output descriptor

**Result**

* nothing

**Examples**

```julia
df = DataFrame(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
dump(df)
```

"""
function Base.dump(io::IO, df::AbstractDataFrame, n::Int, indent)
    println(io, typeof(df), "  $(nrow(df)) observations of $(ncol(df)) variables")
    if n > 0
        for (name, col) in eachcol(df)
            print(io, indent, "  ", name, ": ")
            dump(io, col, n - 1, string(indent, "  "))
        end
    end
end

# summarize the columns of a df
# TODO: clever layout in rows
"""
Summarize the columns of an AbstractDataFrame

```julia
describe(df::AbstractDataFrame)
describe(io, df::AbstractDataFrame)
```

**Arguments**

* `df` : the AbstractDataFrame
* `io` : optional output descriptor

**Result**

* nothing

**Details**

If the column's base type derives from Number, compute the minimum, first
quantile, median, mean, third quantile, and maximum. Missings are filtered and
reported separately.

For boolean columns, report trues, falses, and missings.

For other types, show column characteristics and number of missings.

**Examples**

```julia
df = DataFrame(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
describe(df)
```

"""
StatsBase.describe(df::AbstractDataFrame) = describe(stdout, df)
function StatsBase.describe(io, df::AbstractDataFrame)
    for (name, col) in eachcol(df)
        println(io, name)
        describe(io, col)
        println(io, )
    end
end

function StatsBase.describe(io::IO, X::AbstractVector{Union{T, Missing}}) where T
    missingcount = count(ismissing, X)
    pmissing = 100 * missingcount/length(X)
    if pmissing != 100 && T <: Real
        show(io, StatsBase.summarystats(collect(skipmissing(X))))
    else
        println(io, "Summary Stats:")
    end
    println(io, "Length:         $(length(X))")
    println(io, "Type:           $(eltype(X))")
    !(T <: Real) && println(io, "Number Unique:  $(length(unique(X)))")
    println(io, "Number Missing: $(missingcount)")
    @printf(io, "%% Missing:      %.6f\n", pmissing)
    return
end

##############################################################################
##
## Miscellaneous
##
##############################################################################

function _nonmissing!(res, col)
    # workaround until JuliaLang/julia#21256 is fixed
    eltype(col) >: Missing || return

    @inbounds for (i, el) in enumerate(col)
        res[i] &= !ismissing(el)
    end
end

function _nonmissing!(res, col::CategoricalArray{>: Missing})
    for (i, el) in enumerate(col.refs)
        res[i] &= el > 0
    end
end

"""
Indexes of complete cases (rows without missing values)

```julia
completecases(df::AbstractDataFrame)
```

**Arguments**

* `df` : the AbstractDataFrame

**Result**

* `::Vector{Bool}` : indexes of complete cases

See also [`dropmissing`](@ref) and [`dropmissing!`](@ref).

**Examples**

```julia
df = DataFrame(i = 1:10,
               x = Vector{Union{Missing, Float64}}(rand(10)),
               y = Vector{Union{Missing, String}}(rand(["a", "b", "c"], 10)))
df[[1,4,5], :x] = missing
df[[9,10], :y] = missing
completecases(df)
```

"""
function completecases(df::AbstractDataFrame)
    res = trues(size(df, 1))
    for i in 1:size(df, 2)
        _nonmissing!(res, df[i])
    end
    res
end

"""
Remove rows with missing values.

```julia
dropmissing(df::AbstractDataFrame)
```

**Arguments**

* `df` : the AbstractDataFrame

**Result**

* `::AbstractDataFrame` : the updated copy

See also [`completecases`](@ref) and [`dropmissing!`](@ref).

**Examples**

```julia
df = DataFrame(i = 1:10,
               x = Vector{Union{Missing, Float64}}(rand(10)),
               y = Vector{Union{Missing, String}}(rand(["a", "b", "c"], 10)))
df[[1,4,5], :x] = missing
df[[9,10], :y] = missing
dropmissing(df)
```

"""
dropmissing(df::AbstractDataFrame) = deleterows!(copy(df), findall(!, completecases(df)))

"""
Remove rows with missing values in-place.

```julia
dropmissing!(df::AbstractDataFrame)
```

**Arguments**

* `df` : the AbstractDataFrame

**Result**

* `::AbstractDataFrame` : the updated version

See also [`dropmissing`](@ref) and [`completecases`](@ref).

**Examples**

```julia
df = DataFrame(i = 1:10,
               x = Vector{Union{Missing, Float64}}(rand(10)),
               y = Vector{Union{Missing, String}}(rand(["a", "b", "c"], 10)))
df[[1,4,5], :x] = missing
df[[9,10], :y] = missing
dropmissing!(df)
```

"""
dropmissing!(df::AbstractDataFrame) = deleterows!(df, findall(!, completecases(df)))

"""
    filter(function, df::AbstractDataFrame)

Return a copy of data frame `df` containing only rows for which `function`
returns `true`. The function is passed a `DataFrameRow` as its only argument.

# Examples
```
julia> df = DataFrame(x = [3, 1, 2, 1], y = ["b", "c", "a", "b"])
4×2 DataFrames.DataFrame
│ Row │ x │ y │
├─────┼───┼───┤
│ 1   │ 3 │ b │
│ 2   │ 1 │ c │
│ 3   │ 2 │ a │
│ 4   │ 1 │ b │

julia> filter(row -> row[:x] > 1, df)
2×2 DataFrames.DataFrame
│ Row │ x │ y │
├─────┼───┼───┤
│ 1   │ 3 │ b │
│ 2   │ 2 │ a │
```
"""
Base.filter(f, df::AbstractDataFrame) = df[collect(f(r)::Bool for r in eachrow(df)), :]

"""
    filter!(function, df::AbstractDataFrame)

Remove rows from data frame `df` for which `function` returns `false`.
The function is passed a `DataFrameRow` as its only argument.

# Examples
```
julia> df = DataFrame(x = [3, 1, 2, 1], y = ["b", "c", "a", "b"])
4×2 DataFrames.DataFrame
│ Row │ x │ y │
├─────┼───┼───┤
│ 1   │ 3 │ b │
│ 2   │ 1 │ c │
│ 3   │ 2 │ a │
│ 4   │ 1 │ b │

julia> filter!(row -> row[:x] > 1, df);

julia> df
2×2 DataFrames.DataFrame
│ Row │ x │ y │
├─────┼───┼───┤
│ 1   │ 3 │ b │
│ 2   │ 2 │ a │
```
"""
Base.filter!(f, df::AbstractDataFrame) =
    deleterows!(df, findall(collect(!f(r)::Bool for r in eachrow(df))))

function Base.convert(::Type{Array}, df::AbstractDataFrame)
    convert(Matrix, df)
end
function Base.convert(::Type{Matrix}, df::AbstractDataFrame)
    T = reduce(promote_type, eltypes(df))
    convert(Matrix{T}, df)
end
function Base.convert(::Type{Array{T}}, df::AbstractDataFrame) where T
    convert(Matrix{T}, df)
end
function Base.convert(::Type{Matrix{T}}, df::AbstractDataFrame) where T
    n, p = size(df)
    res = Matrix{T}(undef, n, p)
    idx = 1
    for (name, col) in zip(names(df), columns(df))
        !(T >: Missing) && any(ismissing, col) && error("cannot convert a DataFrame containing missing values to array (found for column $name)")
        copyto!(res, idx, convert(Vector{T}, col))
        idx += n
    end
    return res
end

"""
Indexes of duplicate rows (a row that is a duplicate of a prior row)

```julia
nonunique(df::AbstractDataFrame)
nonunique(df::AbstractDataFrame, cols)
```

**Arguments**

* `df` : the AbstractDataFrame
* `cols` : a column indicator (Symbol, Int, Vector{Symbol}, etc.) specifying the column(s) to compare

**Result**

* `::Vector{Bool}` : indicates whether the row is a duplicate of some
  prior row

See also [`unique`](@ref) and [`unique!`](@ref).

**Examples**

```julia
df = DataFrame(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
df = vcat(df, df)
nonunique(df)
nonunique(df, 1)
```

"""
function nonunique(df::AbstractDataFrame)
    gslots = row_group_slots(df)[3]
    # unique rows are the first encountered group representatives,
    # nonunique are everything else
    res = fill(true, nrow(df))
    @inbounds for g_row in gslots
        (g_row > 0) && (res[g_row] = false)
    end
    return res
end

nonunique(df::AbstractDataFrame, cols::Union{Real, Symbol}) = nonunique(df[[cols]])
nonunique(df::AbstractDataFrame, cols::Any) = nonunique(df[cols])

if isdefined(Base, :unique!)
    import Base.unique!
end

unique!(df::AbstractDataFrame) = deleterows!(df, findall(nonunique(df)))
unique!(df::AbstractDataFrame, cols::Any) = deleterows!(df, findall(nonunique(df, cols)))

# Unique rows of an AbstractDataFrame.
Base.unique(df::AbstractDataFrame) = df[(!).(nonunique(df)), :]
Base.unique(df::AbstractDataFrame, cols::Any) = df[(!).(nonunique(df, cols)), :]

"""
Delete duplicate rows

```julia
unique(df::AbstractDataFrame)
unique(df::AbstractDataFrame, cols)
unique!(df::AbstractDataFrame)
unique!(df::AbstractDataFrame, cols)
```

**Arguments**

* `df` : the AbstractDataFrame
* `cols` :  column indicator (Symbol, Int, Vector{Symbol}, etc.)
specifying the column(s) to compare.

**Result**

* `::AbstractDataFrame` : the updated version of `df` with unique rows.
When `cols` is specified, the return DataFrame contains complete rows,
retaining in each case the first instance for which `df[cols]` is unique.

See also [`nonunique`](@ref).

**Examples**

```julia
df = DataFrame(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
df = vcat(df, df)
unique(df)   # doesn't modify df
unique(df, 1)
unique!(df)  # modifies df
```

"""
(unique, unique!)

# Count the number of missing values in every column of an AbstractDataFrame.
function colmissing(df::AbstractDataFrame) # -> Vector{Int}
    nrows, ncols = size(df)
    missing = zeros(Int, ncols)
    for j in 1:ncols
        missing[j] = countmissing(df[j])
    end
    return missing
end

function without(df::AbstractDataFrame, icols::Vector{Int})
    newcols = setdiff(1:ncol(df), icols)
    df[newcols]
end
without(df::AbstractDataFrame, i::Int) = without(df, [i])
without(df::AbstractDataFrame, c::Any) = without(df, index(df)[c])

##############################################################################
##
## Hcat / vcat
##
##############################################################################

# hcat's first argument must be an AbstractDataFrame
# or AbstractVector if the second argument is AbstractDataFrame
# Trailing arguments (currently) may also be vectors.

# hcat! is defined in DataFrames/DataFrames.jl
# Its first argument (currently) must be a DataFrame.

# catch-all to cover cases where indexing returns a DataFrame and copy doesn't

Base.hcat(df::AbstractDataFrame, x; makeunique::Bool=false) =
    hcat!(df[:, :], x, makeunique=makeunique)
Base.hcat(x, df::AbstractDataFrame; makeunique::Bool=false) =
    hcat!(x, df[:, :], makeunique=makeunique)
Base.hcat(df1::AbstractDataFrame, df2::AbstractDataFrame; makeunique::Bool=false) =
    hcat!(df1[:, :], df2, makeunique=makeunique)
Base.hcat(df::AbstractDataFrame, x, y...; makeunique::Bool=false) =
    hcat!(hcat(df, x, makeunique=makeunique), y..., makeunique=makeunique)
Base.hcat(df1::AbstractDataFrame, df2::AbstractDataFrame, dfn::AbstractDataFrame...;
          makeunique::Bool=false) =
    hcat!(hcat(df1, df2, makeunique=makeunique), dfn..., makeunique=makeunique)

@generated function promote_col_type(cols::AbstractVector...)
    T = mapreduce(x -> Missings.T(eltype(x)), promote_type, cols)
    if CategoricalArrays.iscatvalue(T)
        T = CategoricalArrays.leveltype(T)
    end
    if any(col -> eltype(col) >: Missing, cols)
        if any(col -> col <: AbstractCategoricalArray, cols)
            return :(CategoricalVector{Union{$T, Missing}})
        else
            return :(Vector{Union{$T, Missing}})
        end
    else
        if any(col -> col <: AbstractCategoricalArray, cols)
            return :(CategoricalVector{$T})
        else
            return :(Vector{$T})
        end
    end
end

"""
    vcat(dfs::AbstractDataFrame...)

Vertically concatenate `AbstractDataFrames`.

Column names in all passed data frames must be the same, but they can have
different order. In such cases the order of names in the first passed
`DataFrame` is used.

# Example
```jldoctest
julia> df1 = DataFrame(A=1:3, B=1:3);
julia> df2 = DataFrame(A=4:6, B=4:6);
julia> vcat(df1, df2)
6×2 DataFrames.DataFrame
│ Row │ A │ B │
├─────┼───┼───┤
│ 1   │ 1 │ 1 │
│ 2   │ 2 │ 2 │
│ 3   │ 3 │ 3 │
│ 4   │ 4 │ 4 │
│ 5   │ 5 │ 5 │
│ 6   │ 6 │ 6 │

julia> df1 = DataFrame(A=1:3, B=1:3);
julia> df2 = DataFrame(B=4:6, A=4:6);
julia> vcat(df1, df2)
6×2 DataFrames.DataFrame
│ Row │ A │ B │
├─────┼───┼───┤
│ 1   │ 1 │ 1 │
│ 2   │ 2 │ 2 │
│ 3   │ 3 │ 3 │
│ 4   │ 4 │ 4 │
│ 5   │ 5 │ 5 │
│ 6   │ 6 │ 6 │
```
"""
Base.vcat(df::AbstractDataFrame) = df
Base.vcat(dfs::AbstractDataFrame...) = _vcat(collect(dfs))
function _vcat(dfs::AbstractVector{<:AbstractDataFrame})
    isempty(dfs) && return DataFrame()
    allheaders = map(names, dfs)
    uniqueheaders = unique(allheaders)
    unionunique = union(uniqueheaders...)
    intersectunique = intersect(uniqueheaders...)
    coldiff = setdiff(unionunique, intersectunique)

    if !isempty(coldiff)
        # if any DataFrames are a full superset of names, skip them
        filter!(u -> Set(u) != Set(unionunique), uniqueheaders)
        estrings = Vector{String}(undef, length(uniqueheaders))
        for (i, u) in enumerate(uniqueheaders)
            matching = findall(h -> u == h, allheaders)
            headerdiff = setdiff(coldiff, u)
            cols = join(headerdiff, ", ", " and ")
            args = join(matching, ", ", " and ")
            estrings[i] = "column(s) $cols are missing from argument(s) $args"
        end
        throw(ArgumentError(join(estrings, ", ", ", and ")))
    end

    header = allheaders[1]
    length(header) == 0 && return DataFrame()
    cols = Vector{Any}(undef, length(header))
    for (i, name) in enumerate(header)
        data = [df[name] for df in dfs]
        lens = map(length, data)
        cols[i] = promote_col_type(data...)(undef, sum(lens))
        offset = 1
        for j in 1:length(data)
            copyto!(cols[i], offset, data[j])
            offset += lens[j]
        end
    end
    return DataFrame(cols, header)
end

##############################################################################
##
## Hashing
##
##############################################################################

const hashdf_seed = UInt == UInt32 ? 0xfd8bb02e : 0x6215bada8c8c46de

function Base.hash(df::AbstractDataFrame, h::UInt)
    h += hashdf_seed
    h += hash(size(df))
    for i in 1:size(df, 2)
        h = hash(df[i], h)
    end
    return h
end

## Documentation for methods defined elsewhere

# nrow, ncol
"""
Number of rows or columns in an AbstractDataFrame

```julia
nrow(df::AbstractDataFrame)
ncol(df::AbstractDataFrame)
```

**Arguments**

* `df` : the AbstractDataFrame

**Result**

* `::AbstractDataFrame` : the updated version

See also [`size`](@ref).

NOTE: these functions may be depreciated for `size`.

**Examples**

```julia
df = DataFrame(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
size(df)
nrow(df)
ncol(df)
```

"""


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using Test

Test.@testset "parameters" begin
    include("test_spatialvolume.jl")
    include("test_parambounds.jl")
    include("test_density_sample.jl")
    include("test_parameter_mapping.jl")
end


text:
baremodule IRDumps

function tofiles end
function dump end
function parse end
function parsefile end

module Internal

using ..IRDumps: IRDumps

import InteractiveUtils
using ArgCheck: @argcheck

const HAS_PRINT_LLVM = try
    InteractiveUtils.print_llvm(devnull, "")
    true
catch
    false
end

if HAS_PRINT_LLVM
    using InteractiveUtils: print_llvm
else
    const print_llvm = print
end

include("collection.jl")
include("parse.jl")
include("split.jl")
include("dump.jl")

end  # module Internal

end  # baremodule IRDumps


text:
##Testing the Generation of Kronecker Graph
include("KronGraph500NoPerm.jl")

SCALE = 12
EdgesPerVertex = 16
Nmax = 2.^SCALE
M = EdgesPerVertex .* Nmax

v1,v2 = KronGraph500NoPerm(SCALE,EdgesPerVertex)
A = sparse(v1,v2,1)

#There is a spy in PyPlot
#=
using PyPlot

figure
spy(A)
xlabel("end vertex")
ylabel("start vertex")

figure
loglog(full(OutDegree(A)),'o')
xlabel("degree")
ylabel("count")

spy(A)

=#


text:
# Testing part 2
using Test, ReferenceTests

# for testing, comment the function call: TM_2D(; do_visu = true) at the end of TM_2D_perf.jl
include("../scripts-part2/TM_2D_perf.jl")

# ------------------test macros defined in TM_2D_perf.jl---------------------
# Initialize arbitrary values for varialbes necessary for the macros
λb, θp, θv, θT = 1.0, 1.0, 1.0, 1.0
dx, dy = 1.0, 1.0
nx, ny = 5, 5                               # including ghost cells, choose small values for test speed
Vx, Vy = rand(nx-1,ny), rand(nx  ,ny-1)
ηc, ηv = rand(nx  ,ny), rand(nx-1,ny-1)
P , T  = rand(nx  ,ny), rand(nx  ,ny  )
# Calculate variables using macros
∇V   = [  @∇V(ix,iy) for ix=2:nx-1, iy=2:ny-1 ]  # inner cells, size (nx-2)(ny-2)
ϵxx  = [ @ϵxx(ix,iy) for ix=2:nx-1, iy=2:ny-1 ]  # inner cells
ϵyy  = [ @ϵyy(ix,iy) for ix=2:nx-1, iy=2:ny-1 ]  # inner cells
ϵxy  = [ @ϵxy(ix,iy) for ix=1:nx-1, iy=1:ny-1 ]  # cell vertices, size (nx-1)(ny-1)
τxx  = [ @τxx(ix,iy) for ix=2:nx-1, iy=2:ny-1 ]  # inner cells
τyy  = [ @τyy(ix,iy) for ix=2:nx-1, iy=2:ny-1 ]  # inner cells
τxy  = [ @τxy(ix,iy) for ix=1:nx-1, iy=1:ny-1 ]  # cell vertices, size (nx-1)(ny-1)
dτP  = [ @dτP(ix,iy) for ix=2:nx-1, iy=2:ny-1 ]  # inner cells
dτVx = [@dτVx(ix,iy) for ix=2:nx-2, iy=2:ny-1 ]  # inner Vx, size (nx-3)(ny-2)
dτVy = [@dτVy(ix,iy) for ix=2:nx-1, iy=2:ny-2 ]  # inner Vy, size (nx-2)(ny-3)
dτT  = [ @dτT(ix,iy) for ix=2:nx-1, iy=2:ny-1 ]  # inner cells
# calculate variables using definitions (as references)
 ∇V_ref = diff(Vx[:,2:end-1],dims=1)./dx .+ diff(Vy[2:end-1,:],dims=2)./dy
ϵxx_ref = diff(Vx[:,2:end-1],dims=1)./dx .- 1/2 .* ∇V_ref
ϵyy_ref = diff(Vy[2:end-1,:],dims=2)./dy .- 1/2 .* ∇V_ref
ϵxy_ref = 0.5.*(diff(Vx,dims=2)./dy .+ diff(Vy,dims=1)./dx)
τxx_ref = 2 .* ηc[2:end-1,2:end-1] .*(ϵxx_ref .+ λb.*∇V_ref) .- P[2:end-1,2:end-1]
τyy_ref = 2 .* ηc[2:end-1,2:end-1] .*(ϵyy_ref .+ λb.*∇V_ref) .- P[2:end-1,2:end-1]
τxy_ref = 2 .* ηv .* ϵxy_ref

 dτP_ref = θp*4.1/max(nx-2,ny-2).*  ηc[2:end-1,2:end-1] .* (1.0+λb)
dτVx_ref = θv/4.1*(min(dx,dy)^2 ./( 0.5.*(ηc[2:end-2,2:end-1] + ηc[3:end-1,2:end-1]) ))./(1+λb)
dτVy_ref = θv/4.1*(min(dx,dy)^2 ./( 0.5.*(ηc[2:end-1,2:end-2] + ηc[2:end-1,3:end-1]) ))./(1+λb)
 dτT_ref = [θT/4.1* min(dx,dy)^2 for ix=2:nx-1, iy=2:ny-1]
# Test macro calculations against reference
tol = 1e-10;
@testset "Test strain rate macros" begin
    @test isapprox(∇V ,  ∇V_ref, atol=tol)
    @test isapprox(ϵxx, ϵxx_ref, atol=tol)
    @test isapprox(ϵyy, ϵyy_ref, atol=tol)
    @test isapprox(ϵxy, ϵxy_ref, atol=tol)
end

@testset "Test stress macros" begin
    @test isapprox(τxx, τxx_ref, atol=tol)
    @test isapprox(τyy, τyy_ref, atol=tol)
    @test isapprox(τxy, τxy_ref, atol=tol)
end

@testset "Test pseudo time step macros" begin
    @test isapprox(dτP ,  dτP_ref, atol=tol)
    @test isapprox(dτVx, dτVx_ref, atol=tol)
    @test isapprox(dτVy, dτVy_ref, atol=tol)
    @test isapprox(dτT ,  dτT_ref, atol=tol)
end


text:
if !(@isdefined FramePointCache)
    @computed struct FramePointCache{n,S}
        link_cache::fulltype(LinkFrameCache{n,S})
        point::Point3D{SVector{3,S}}
        world_position_jacobian::CacheElement{Vector{SMatrix{3,12,S,36}}}
        world_position_jacobian_dot::CacheElement{Vector{SMatrix{3,12,S,36}}}
    end
end

function FramePointCache(link_cache::LinkFrameCache{n,S}, position_on_link) where {n,S}
    point = Point3D(default_frame(link_cache.link), position_on_link)
    world_position_jacobian = CacheElement([@SMatrix zeros(3,12)])
    world_position_jacobian_dot = CacheElement([@SMatrix zeros(3,12)])
    FramePointCache{n,S}(link_cache, point, world_position_jacobian, world_position_jacobian_dot)
end

function setdirty!(point_cache::FramePointCache)
    CustomCollections.setdirty!(point_cache.world_position_jacobian)
    CustomCollections.setdirty!(point_cache.world_position_jacobian_dot)
    nothing
end

# Frame transform and position in frame to world frame position
function point_world_position(Rt, point_frame_position)
    R = SMatrix{3,3,eltype(Rt)}(Rt[static(1):static(9)])
    t = Rt[static(10):static(12)]
    R*point_frame_position + t
end

function point_world_position_jacobian(Rt, point_frame_position)
    f = (Rt) -> point_world_position(Rt, point_frame_position)
    ForwardDiff.jacobian(f, Rt)
end

function point_world_position_jacobian_dot(Rt, Rt_dot, point_frame_position)
    m, n = 12, 3
    ∇xf = SMatrix{n*m,m,eltype(Rt)}(ForwardDiff.jacobian(Rt -> reshape(
        point_world_position_jacobian(Rt, point_frame_position), m*n), Rt))
    Jfdot = reshape(∇xf*Rt_dot, Size(n,m))
end

@inline function update_point_world_position_jacobian!(world_position_jacobian::CacheElement, Rt,
        point_frame_position)
    isdirty(world_position_jacobian) && _update_point_world_position_jacobian!(world_position_jacobian, Rt, point_frame_position)
    nothing
end

function _update_point_world_position_jacobian!(world_position_jacobian::CacheElement, Rt,
        point_frame_position)
    f = (Rt) -> point_world_position(Rt, point_frame_position)
    world_position_jacobian.data[1] = ForwardDiff.jacobian(f, Rt)
    world_position_jacobian.dirty = false
    nothing
end

@inline function update_point_world_position_jacobian_dot!(
    world_position_jacobian_dot::CacheElement, Rt, Rt_dot, point_frame_position)
    isdirty(world_position_jacobian_dot) && _update_point_world_position_jacobian_dot!(world_position_jacobian_dot, Rt, Rt_dot, point_frame_position)
    nothing
end

function _update_point_world_position_jacobian_dot!(world_position_jacobian_dot::CacheElement,
        Rt, Rt_dot, point_frame_position)
    n, m = size(world_position_jacobian_dot.data[1])
    f = (Rt) -> point_world_position(Rt, point_frame_position)
    Jf = (Rt) -> ForwardDiff.jacobian(f, Rt)
    ∇xf = ForwardDiff.jacobian(Rt -> reshape(Jf(Rt), m*n), Rt)
    Jfdot = SMatrix{n,m,eltype(Rt)}(reshape(∇xf*Rt, n, m))::SArray{Tuple{n,m},eltype(Rt),2,m*n}

    world_position_jacobian_dot.data[1] = Jfdot
    world_position_jacobian_dot.dirty = false
    nothing
end


text:
import RDatasets
iris = RDatasets.dataset("datasets", "iris")
using StatsPlots, Interact
using Blink
w = Window()
body!(w, dataviewer(iris))