text:
## COV_EXCL_START

# TODO
# - serial version for lower latency
# - group-stride loop to delay need for second kernel launch

# Reduce a value across a group, using local memory for communication
@inline function reduce_group(op, val::T, neutral, ::Val{maxitems}) where {T, maxitems}
    items = workgroupDim().x
    item = workitemIdx().x

    # shared mem for a complete reduction
    shared = ROCDeviceArray((2*maxitems,), alloc_special(Val(:reduce_block), T, Val(AS.Local), Val(2*maxitems)))
    @inbounds shared[item] = val

    # perform a reduction
    d = items>>1
    while d > 0
        sync_workgroup()
        if item <= d
            shared[item] = op(shared[item], shared[item+d])
        end
        d >>= 1
    end

    # load the final value on the first item
    if item == 1
        val = @inbounds shared[item]
    end

    return val
end

Base.@propagate_inbounds _map_getindex(args::Tuple, I) = ((args[1][I]), _map_getindex(Base.tail(args), I)...)
Base.@propagate_inbounds _map_getindex(args::Tuple{Any}, I) = ((args[1][I]),)
Base.@propagate_inbounds _map_getindex(args::Tuple{}, I) = ()

# Reduce an array across the grid. All elements to be processed can be addressed by the
# product of the two iterators `Rreduce` and `Rother`, where the latter iterator will have
# singleton entries for the dimensions that should be reduced (and vice versa).
function partial_mapreduce_device(f, op, neutral, maxitems, Rreduce, Rother, R, As...)
    # decompose the 1D hardware indices into separate ones for reduction (across items
    # and possibly groups if it doesn't fit) and other elements (remaining groups)
    localIdx_reduce = workitemIdx().x
    localDim_reduce = workgroupDim().x
    groupIdx_reduce, groupIdx_other = fldmod1(workgroupIdx().x, length(Rother))
    groupDim_reduce = gridDimWG().x ÷ length(Rother)

    # group-based indexing into the values outside of the reduction dimension
    # (that means we can safely synchronize items within this group)
    iother = groupIdx_other
    @inbounds if iother <= length(Rother)
        Iother = Rother[iother]

        # load the neutral value
        Iout = CartesianIndex(Tuple(Iother)..., groupIdx_reduce)
        neutral = if neutral === nothing
            R[Iout]
        else
            neutral
        end

        val = op(neutral, neutral)

        # reduce serially across chunks of input vector that don't fit in a group
        ireduce = localIdx_reduce + (groupIdx_reduce - 1) * localDim_reduce
        while ireduce <= length(Rreduce)
            Ireduce = Rreduce[ireduce]
            J = Base.max(Iother, Ireduce)
            val = op(val, f(_map_getindex(As, J)...))
            ireduce += localDim_reduce * groupDim_reduce
        end

        val = reduce_group(op, val, neutral, maxitems)

        # write back to memory
        if localIdx_reduce == 1
            R[Iout] = val
        end
    end

    return
end

## COV_EXCL_STOP

function GPUArrays.mapreducedim!(f::F, op::OP, R::ROCArray{T},
                                 A::Union{AbstractArray,Broadcast.Broadcasted};
                                 init=nothing) where {F, OP, T}
    Base.check_reducedims(R, A)
    length(A) == 0 && return R # isempty(::Broadcasted) iterates

    # add singleton dimensions to the output container, if needed
    if ndims(R) < ndims(A)
        dims = Base.fill_to_length(size(R), 1, Val(ndims(A)))
        R = reshape(R, dims)
    end

    # iteration domain, split in two: one part covers the dimensions that should
    # be reduced, and the other covers the rest. combining both covers all values.
    Rall = CartesianIndices(axes(A))
    Rother = CartesianIndices(axes(R))
    Rreduce = CartesianIndices(ifelse.(axes(A) .== axes(R), Ref(Base.OneTo(1)), axes(A)))
    # NOTE: we hard-code `OneTo` (`first.(axes(A))` would work too) or we get a
    #       CartesianIndices object with UnitRanges that behave badly on the GPU.
    @assert length(Rall) == length(Rother) * length(Rreduce)

    # allocate an additional, empty dimension to write the reduced value to.
    # this does not affect the actual location in memory of the final values,
    # but allows us to write a generalized kernel supporting partial reductions.
    R′ = reshape(R, (size(R)..., 1))

    # how many items do we want?
    #
    # items in a group work together to reduce values across the reduction dimensions;
    # we want as many as possible to improve algorithm efficiency and execution occupancy.
    wanted_items = nextpow(2, length(Rreduce))
    function compute_items(max_items)
        if wanted_items > max_items
            prevpow(2, max_items)
        else
            wanted_items
        end
    end

    # how many items can we launch?
    #
    # we might not be able to launch all those items to reduce each slice in one go.
    # that's why each items also loops across their inputs, processing multiple values
    # so that we can span the entire reduction dimension using a single item group.

    # group size is restricted by local memory
    agent = R.buf.agent
    pools = filter(pool->pool_segment(pool) == HSA.AMD_SEGMENT_GROUP, memory_pools(agent))
    max_items = if !isempty(pools)
        pool = first(pools)
        max_lmem_elements = pool_size(pool) ÷ sizeof(T)
        isa = first(isas(agent))
        Base.min(max_group_size(isa), compute_items(max_lmem_elements ÷ 2))
    else
        @warn "No group segment detected for agent $agent; assuming 64 elements\nThis message will not be shown again" maxlog=1
        64
    end
    # TODO: dynamic local memory to avoid two compilations

    #= TODO: let the runtime suggest a group size
    args = (f, op, init, Val(max_items), Rreduce, Rother, R′, A)
    kernel_args = rocconvert.(args)
    kernel_tt = Tuple{Core.Typeof.(kernel_args)...}
    kernel = rocfunction(partial_mapreduce_device, kernel_tt)
    reduce_items = compute_items(suggest_groupsize(kernel.fun, wanted_items).x)
    =#
    reduce_items = max_items

    # how many groups should we launch?
    #
    # even though we can always reduce each slice in a single item group, that may not be
    # optimal as it might not saturate the GPU. we already launch some groups to process
    # independent dimensions in parallel; pad that number to ensure full occupancy.
    other_groups = length(Rother)
    reduce_groups = cld(length(Rreduce), reduce_items)

    # determine the launch configuration
    items = reduce_items
    groups = reduce_groups*other_groups
    gridsize = items*groups

    # perform the actual reduction
    if reduce_groups == 1
        # we can cover the dimensions to reduce using a single group
        wait(@roc gridsize=gridsize groupsize=items partial_mapreduce_device(
            f, op, init, Val(items), Rreduce, Rother, R′, A))
    else
        # we need multiple steps to cover all values to reduce
        partial = similar(R, (size(R)..., reduce_groups))
        if init === nothing
            # without an explicit initializer we need to copy from the output container
            sz = prod(size(R))
            for i in 1:reduce_groups
                # TODO: async copies (or async fill!, but then we'd need to load first)
                #       or maybe just broadcast since that extends singleton dimensions
                copyto!(partial, (i-1)*sz+1, R, 1, sz)
            end
        end
        wait(@roc gridsize=gridsize groupsize=items partial_mapreduce_device(
            f, op, init, Val(items), Rreduce, Rother, partial, A))

        GPUArrays.mapreducedim!(identity, op, R′, partial; init=init)
    end

    return R
end


text:
using BenchmarkTools
using CatViews
import Base.Cartesian: @ntuple

x = CatView( @ntuple 1000 (n)->randn(10) )

@benchmark begin
    for xi in $x
        nothing
    end
end

ext: jl


text:
module argminUniProb
  export argminProb
  using Statistics

  """
      argminProb( r::Matrix{Float64} )

      compute minimum value per row, if several cols hold the minimum value, randomly choose among them.
      #Arguments
        - `r::Matrix{Float64}`: risk matrix each column represents one hypothesis.
      #Return
        - `Vector{Float64}`:    indices of maximal values.
  """
  function argminProb( r::Matrix{Float64} )
    index = Vector{Int64}( undef, size( r, 1 ) );
    for i in 1:1:size( index, 1 )
      index[i] = argminProb( r[i,:] );
    end
    return index;
  end

  """
      argminProb!( r::Matrix{Float64}, v::Vector{Int64} )

      compute minimum value per row, if several cols have maximum value randomly choose among them.
      #Arguments
        - `r::Matrix{Float64}`: risk matrix each column represents one hypothesis.
        - `v::Vector{Int64} `:  resulting vector.
      #Return
        - `nothing`:            nothing.
  """
  function argminProb!( r::Matrix{Float64}, v::AbstractVector )
    @assert size( v,1 ) == size( r, 2 );
    for i in 1:1:size( v, 1 )
      v[i] = argminProb( r[i,:] );
    end
  end

  """
      argminProb( values::Vector{Float64} )

      compute minimum, if several indices have minimum value randomly choose among them.
      #Arguments
        - `values::Vector{Float64}`: values to compute argmin on.
      #Return
        - `Vector{Float64}`:         indices of minimal values.
  """
  function argminProb( values::AbstractVector )
    minPositions = collect( 1:1:size( values, 1 ) )[ values .== minimum( values ) ];
    if size( minPositions, 1 ) > 1
      res = minPositions[ rand( 1:size( minPositions, 1 ) ) ];
    else
      res = minPositions[1];
    end
    return res;
  end
end


text:
#=

ДАНО: Робот - в произвольной клетке поля (без внутренних перегородок и маркеров)
РЕЗУЛЬТАТ: Робот - в исходном положении, и все клетки по периметру внешней рамки промакированы

=#

# глабальная массив coords, содержащий координаты стартовой клетки Робота
coords = []

using HorizonSideRobots
originRobot = Robot(animate = true, "2example.sit")
include("../functions.jl")

# get_line! создаёт линию из маркеров
function get_line!(robot::Robot,side::HorizonSide)::Nothing
    while (isborder(robot,side) == false)
        move!(robot,side)
        putmarker!(robot)
    end
end

# Point2_master! главная функция программы
function master2!(robot::Robot)::Nothing
    back_path = moving_in_angle!(robot)

    for side in [Nord,Ost,Sud,West]
        get_line!(robot,side)
    end

    moving_back_to_start!(robot, back_path)
end

master2!(originRobot)

max_stars_repo_name: Igetsuga/AlgorithmicsHomework

max_issues_repo_name: Igetsuga/AlgorithmicsHomework

max_forks_repo_name: Igetsuga/AlgorithmicsHomework


text:
module Luna
import FFTW
import Hankel
import Logging
import LinearAlgebra: mul!, ldiv!
Logging.disable_logging(Logging.BelowMinLevel)

"""
    HDF5LOCK

Lock on the HDF5 library for multi-threaded execution.
"""
const HDF5LOCK = ReentrantLock()

"""
    @hlock

Wait for HDF5LOCK, execute the expression, and release H5DFLOCK.

!!! warning
    For thread safety, any call to functions from HDF5.jl needs to be preceeded by @hlock.
"""
macro hlock(expr)
    quote
        try
            lock(HDF5LOCK)
            $(esc(expr))
        finally
            unlock(HDF5LOCK)
        end
    end
end

"""
    Luna.settings

Dictionary of global settings for `Luna`.
"""
settings = Dict{String, Any}("fftw_flag" => FFTW.PATIENT,
                             "fftw_threads" => 0)

"""
    set_fftw_mode(mode)

Set FFTW planning mode for all FFTW transform planning in `Luna`.

Possible values for `mode` are `:estimate`, `:measure`, `:patient`, and `:exhaustive`.
The initial value upon loading `Luna` is `:patient`

# Examples
```jldoctest
julia> Luna.set_fftw_mode(:patient)
0x00000020
```
"""
function set_fftw_mode(mode)
    s = uppercase(string(mode))
    flag = getfield(FFTW, Symbol(s))
    settings["fftw_flag"] = flag
end

"""
    set_fftw_threads(nthr)

Set number of threads to be used by FFTW. If set to `0`, the number of threads used by
FFTW is determined automatically (see [`Utils.FFTWthreads()`](@ref))
"""
function set_fftw_threads(nthr=0)
    settings["fftw_threads"] = nthr
    FFTW.set_num_threads(Utils.FFTWthreads())
end

function __init__()
    set_fftw_threads()
end

include("Utils.jl")
include("Scans.jl")
include("Output.jl")
include("Maths.jl")
include("PhysData.jl")
include("Grid.jl")
include("Modes.jl")
include("Fields.jl")
include("RK45.jl")
include("LinearOps.jl")
include("Capillary.jl")
include("Antiresonant.jl")
include("RectModes.jl")
include("Nonlinear.jl")
include("Ionisation.jl")
include("NonlinearRHS.jl")
include("Processing.jl")
include("Stats.jl")
include("Polarisation.jl")
include("Tools.jl")
include("Plotting.jl")
include("Raman.jl")
include("SFA.jl")
include("Interface.jl")

prop_capillary = Interface.prop_capillary
Pulses = Interface.Pulses

Scan = Scans.Scan
runscan = Scans.runscan
makefilename = Scans.makefilename
addvariable! = Scans.addvariable!

export Utils, Scans, Output, Maths, PhysData, Grid, RK45, Modes, Capillary, RectModes,
       Nonlinear, Ionisation, NonlinearRHS, LinearOps, Stats, Polarisation,
       Tools, Plotting, Raman, Antiresonant, Fields, Processing, Interface, SFA,
       prop_capillary, Pulses, Scan, runscan, makefilename, addvariable!

# for a tuple of TimeFields we assume all inputs are for mode 1
function doinput_sm(grid, inputs::Tuple{Vararg{T} where T <: Fields.TimeField}, FT)
    out = fill(0.0 + 0.0im, length(grid.ω))
    for field in inputs
        out .+= field(grid, FT)
    end
    return out
end

# for a single Fields.TimeField we assume a single input for mode 1
function doinput_sm(grid, inputs::Fields.TimeField, FT)
    doinput_sm(grid, (inputs,), FT)
end

function doinput_sm(grid, inputs::Tuple{Vararg{T} where T <: NamedTuple{<:Any, <:Tuple{Vararg{Any}}}}, FT)
    if any([i.mode ≠ 1 for i in inputs])
        error("For mode-averaged propagation, all inputs must be in 1st mode.")
    end
    inputs_flat = Tuple(Iterators.flatten([i.fields for i in inputs]))
    doinput_sm(grid, inputs_flat, FT)
end

function setup(grid::Grid.RealGrid, densityfun, responses, inputs, βfun!, aeff;
               norm! = NonlinearRHS.norm_mode_average(grid, βfun!, aeff))
    Utils.loadFFTwisdom()
    xo = Array{Float64}(undef, length(grid.to))
    FTo = FFTW.plan_rfft(xo, 1, flags=settings["fftw_flag"])
    transform = NonlinearRHS.TransModeAvg(grid, FTo, responses, densityfun, norm!, aeff)
    x = Array{Float64}(undef, length(grid.t))
    FT = FFTW.plan_rfft(x, 1, flags=settings["fftw_flag"])
    Eω = doinput_sm(grid, inputs, FT)
    inv(FT) # create inverse FT plans now, so wisdom is saved
    inv(FTo)
    Utils.saveFFTwisdom()
    Eω, transform, FT
end

function setup(grid::Grid.EnvGrid, densityfun, responses, inputs, βfun!, aeff;
               norm! = NonlinearRHS.norm_mode_average(grid, βfun!, aeff))
    Utils.loadFFTwisdom()
    x = Array{ComplexF64}(undef, length(grid.t))
    FT = FFTW.plan_fft(x, 1, flags=settings["fftw_flag"])
    xo = Array{ComplexF64}(undef, length(grid.to))
    FTo = FFTW.plan_fft(xo, 1, flags=settings["fftw_flag"])
    transform = NonlinearRHS.TransModeAvg(grid, FTo, responses, densityfun, norm!, aeff)
    Eω = doinput_sm(grid, inputs, FT)
    inv(FT) # create inverse FT plans now, so wisdom is saved
    inv(FTo)
    Utils.saveFFTwisdom()
    Eω, transform, FT
end

# for a tuple of NamedTuple's with tuple fields we assume all is well
function doinput_mm!(Eω, grid, inputs::Tuple{Vararg{T} where T <: NamedTuple{<:Any, <:Tuple{Vararg{Any}}}}, FT)
    for input in inputs
        out = @view Eω[:, input.mode]
        for field in input.fields
            out .+= field(grid, FT)
        end
    end
end

# for a tuple of TimeFields we assume all inputs are for mode 1
function doinput_mm!(Eω, grid, inputs::Tuple{Vararg{T} where T <: Fields.TimeField}, FT)
    doinput_mm!(Eω, grid, ((mode=1, fields=inputs),), FT)
end

# for a single Fields.TimeField we assume a single input for mode 1
function doinput_mm!(Eω, grid, inputs::Fields.TimeField, FT)
    doinput_mm!(Eω, grid, ((mode=1, fields=(inputs,)),), FT)
end

function setup(grid::Grid.RealGrid, densityfun, responses, inputs,
               modes::Modes.ModeCollection, components;
               full=false, norm! = NonlinearRHS.norm_modal(grid))
    ts = Modes.ToSpace(modes, components=components)
    Utils.loadFFTwisdom()
    xt = Array{Float64}(undef, length(grid.t))
    FTt = FFTW.plan_rfft(xt, 1, flags=settings["fftw_flag"])
    Eω = zeros(ComplexF64, length(grid.ω), length(modes))
    doinput_mm!(Eω, grid, inputs, FTt)
    x = Array{Float64}(undef, length(grid.t), length(modes))
    FT = FFTW.plan_rfft(x, 1, flags=settings["fftw_flag"])
    xo = Array{Float64}(undef, length(grid.to), ts.npol)
    FTo = FFTW.plan_rfft(xo, 1, flags=settings["fftw_flag"])
    transform = NonlinearRHS.TransModal(grid, ts, FTo,
                                 responses, densityfun, norm!,
                                 rtol=1e-3, atol=0.0, mfcn=300, full=full)
    inv(FT) # create inverse FT plans now, so wisdom is saved
    inv(FTo)
    Utils.saveFFTwisdom()
    Eω, transform, FT
end

function setup(grid::Grid.EnvGrid, densityfun, responses, inputs,
               modes::Modes.ModeCollection, components;
               full=false, norm! = NonlinearRHS.norm_modal(grid))
    ts = Modes.ToSpace(modes, components=components)
    Utils.loadFFTwisdom()
    xt = Array{ComplexF64}(undef, length(grid.t))
    FTt = FFTW.plan_fft(xt, 1, flags=settings["fftw_flag"])
    Eω = zeros(ComplexF64, length(grid.ω), length(modes))
    doinput_mm!(Eω, grid, inputs, FTt)
    x = Array{ComplexF64}(undef, length(grid.t), length(modes))
    FT = FFTW.plan_fft(x, 1, flags=settings["fftw_flag"])
    xo = Array{ComplexF64}(undef, length(grid.to), ts.npol)
    FTo = FFTW.plan_fft(xo, 1, flags=settings["fftw_flag"])
    transform = NonlinearRHS.TransModal(grid, ts, FTo,
                                 responses, densityfun, norm!,
                                 rtol=1e-3, atol=0.0, mfcn=300, full=full)
    inv(FT) # create inverse FT plans now, so wisdom is saved
    inv(FTo)
    Utils.saveFFTwisdom()
    Eω, transform, FT
end

function doinputs_fs!(Eωk, grid, spacegrid::Union{Hankel.QDHT,Grid.FreeGrid}, FT,
                   inputs::Tuple{Vararg{T} where T <: Fields.SpatioTemporalField})
    for field in inputs
        Eωk .+= field(grid, spacegrid, FT)
    end
end

function doinputs_fs!(Eωk, grid, spacegrid::Union{Hankel.QDHT,Grid.FreeGrid}, FT,
                   inputs::Fields.SpatioTemporalField)
    doinputs_fs!(Eωk, grid, spacegrid, FT, (inputs,))
end

function setup(grid::Grid.RealGrid, q::Hankel.QDHT,
               densityfun, normfun, responses, inputs)
    Utils.loadFFTwisdom()
    xt = zeros(Float64, length(grid.t), length(q.r))
    FT = FFTW.plan_rfft(xt, 1, flags=settings["fftw_flag"])
    Eω = zeros(ComplexF64, length(grid.ω), length(q.k))
    Eωk = q * Eω
    doinputs_fs!(Eωk, grid, q, FT, inputs)
    xo = Array{Float64}(undef, length(grid.to), length(q.r))
    FTo = FFTW.plan_rfft(xo, 1, flags=settings["fftw_flag"])
    transform = NonlinearRHS.TransRadial(grid, q, FTo, responses, densityfun, normfun)
    inv(FT) # create inverse FT plans now, so wisdom is saved
    inv(FTo)
    Utils.saveFFTwisdom()
    Eωk, transform, FT
end

function setup(grid::Grid.EnvGrid, q::Hankel.QDHT,
               densityfun, normfun, responses, inputs)
    Utils.loadFFTwisdom()
    xt = zeros(ComplexF64, length(grid.t), length(q.r))
    FT = FFTW.plan_fft(xt, 1, flags=settings["fftw_flag"])
    Eω = zeros(ComplexF64, length(grid.ω), length(q.k))
    Eωk = q * Eω
    doinputs_fs!(Eωk, grid, q, FT, inputs)
    xo = Array{ComplexF64}(undef, length(grid.to), length(q.r))
    FTo = FFTW.plan_fft(xo, 1, flags=settings["fftw_flag"])
    transform = NonlinearRHS.TransRadial(grid, q, FTo, responses, densityfun, normfun)
    inv(FT) # create inverse FT plans now, so wisdom is saved
    inv(FTo)
    Utils.saveFFTwisdom()
    Eωk, transform, FT
end

function setup(grid::Grid.RealGrid, xygrid::Grid.FreeGrid,
               densityfun, normfun, responses, inputs)
    Utils.loadFFTwisdom()
    x = xygrid.x
    y = xygrid.y
    xr = Array{Float64}(undef, length(grid.t), length(y), length(x))
    FT = FFTW.plan_rfft(xr, (1, 2, 3), flags=settings["fftw_flag"])
    Eωk = zeros(ComplexF64, length(grid.ω), length(y), length(x))
    doinputs_fs!(Eωk, grid, xygrid, FT, inputs)
    xo = Array{Float64}(undef, length(grid.to), length(y), length(x))
    FTo = FFTW.plan_rfft(xo, (1, 2, 3), flags=settings["fftw_flag"])
    transform = NonlinearRHS.TransFree(grid, xygrid, FTo,
                                       responses, densityfun, normfun)
    inv(FT) # create inverse FT plans now, so wisdom is saved
    inv(FTo)
    Utils.saveFFTwisdom()
    Eωk, transform, FT
end

function setup(grid::Grid.EnvGrid, xygrid::Grid.FreeGrid,
               densityfun, normfun, responses, inputs)
    Utils.loadFFTwisdom()
    x = xygrid.x
    y = xygrid.y
    xr = Array{ComplexF64}(undef, length(grid.t), length(y), length(x))
    FT = FFTW.plan_fft(xr, (1, 2, 3), flags=settings["fftw_flag"])
    Eωk = zeros(ComplexF64, length(grid.ω), length(y), length(x))
    doinputs_fs!(Eωk, grid, xygrid, FT, inputs)
    xo = Array{ComplexF64}(undef, length(grid.to), length(y), length(x))
    FTo = FFTW.plan_fft(xo, (1, 2, 3), flags=settings["fftw_flag"])
    transform = NonlinearRHS.TransFree(grid, xygrid, FTo,
                                       responses, densityfun, normfun)
    inv(FT) # create inverse FT plans now, so wisdom is saved
    inv(FTo)
    Utils.saveFFTwisdom()
    Eωk, transform, FT
end

linoptype(l::AbstractArray) = "constant"
linoptype(l) = "variable"

gridtype(g::Grid.RealGrid) = "field-resolved"
gridtype(g::Grid.EnvGrid) = "envelope"
gridtype(g) = "unknown"

simtype(g, t, l) = Dict("field" => gridtype(g),
                        "transform" => string(t),
                        "linop" => linoptype(l))

function dumps(t, l)
    io = IOBuffer()
    dump(io, t)
    tr = String(take!(io))
    io = IOBuffer()
    dump(io, l)
    lo = String(take!(io))
    Dict("transform" => tr, "linop" => lo)
end

function run(Eω, grid,
             linop, transform, FT, output;
             min_dz=0, max_dz=grid.zmax/2, init_dz=1e-4, z0=0.0,
             rtol=1e-6, atol=1e-10, safety=0.9, norm=RK45.weaknorm,
             status_period=1)

    Et = FT \ Eω

    function stepfun(Eω, z, dz, interpolant)
        Eω .*= grid.ωwin
        ldiv!(Et, FT, Eω)
        Et .*= grid.twin
        mul!(Eω, FT, Et)
        output(Eω, z, dz, interpolant)
    end

    # check_cache does nothing except for HDF5Outputs
    Eωc, zc, dzc = Output.check_cache(output, Eω, z0, init_dz)
    if zc > z0
        Logging.@info("Found cached propagation. Resuming...")
        Eω, z0, init_dz = Eωc, zc, dzc
    end

    output(Grid.to_dict(grid), group="grid")
    output(simtype(grid, transform, linop), group="simulation_type")
    output(dumps(transform, linop), group="dumps")

    RK45.solve_precon(
        transform, linop, Eω, z0, init_dz, grid.zmax, stepfun=stepfun,
        max_dt=max_dz, min_dt=min_dz,
        rtol=rtol, atol=atol, safety=safety, norm=norm,
        status_period=status_period)
end

end # module


text:
using Test
using StaticArrays
using Polyhedra

function test_planar_square(VT, d)
    ei = ntuple(i -> Polyhedra.basis(VT, d, i), 2)
    so = -ei[1] + ei[2]
    no = -ei[1] - ei[2]
    ne = ei[1] - ei[2]
    se = ei[1] + ei[2]
    expected = [no, so, se, ne]
    for ps in [
        [no, so, ne, se],
        [no, se, so, ne],
        [se, so, no, ne],
        expected
    ]
        hull = Polyhedra.planar_hull(vrep(ps))
        @test collect(points(hull)) == expected
        @test !hasallrays(hull)
        v = vrep(ps) + conichull(se)
        hull = Polyhedra.planar_hull(v)
        @test collect(points(hull)) == [ne, no, so]
        @test collect(rays(hull)) == [Ray(se)]
        @test !haslines(hull)
        v = vrep(ps) + conichull(se, se)
        hull = Polyhedra.planar_hull(v)
        @test collect(points(hull)) == [ne, no, so]
        @test collect(rays(hull)) == [Ray(se)]
        @test !haslines(hull)
        ne = ei[1] - ei[2]
        for v in [vrep(ps) + conichull(se, ne),
                  vrep(ps) + conichull(se, ne, ei[1])]
            hull = Polyhedra.planar_hull(v)
            @test collect(points(hull)) == [no, so]
            @test collect(rays(hull)) == [Ray(ne), Ray(se)]
            @test !haslines(hull)
        end
        v = vrep(ps) + conichull(se, ne, -se)
        hull = Polyhedra.planar_hull(v)
        @test collect(points(hull)) == [ne]
        @test collect(rays(hull)) == [Ray(ne)]
        @test collect(lines(hull)) == [Line(se)]
        for v in [vrep(ps) + conichull(se, ne, -ne),
                  vrep(ps) + conichull(se, Line(ne)),
                  vrep(ps) + conichull(ne, se, -ne)]
            hull = Polyhedra.planar_hull(v)
            @test collect(points(hull)) == [no]
            @test collect(rays(hull)) == [Ray(se)]
            @test collect(lines(hull)) == [Line(ne)]
        end
        for v in [vrep(ps) + conichull(se, ne, -ne, -se),
                  vrep(ps) + conichull(se, ne, -ei[1])]
            hull = Polyhedra.planar_hull(v)
            @test collect(points(hull)) == [no]
            @test !hasrays(hull)
            @test collect(lines(hull)) == [Line(ne), Line(se)]
        end
    end
end

@testset "Planar" begin
    for T in [Float64, Int]
        for (VT, d) in [(Vector{T}, 2), (SVector{2, T}, StaticArrays.Size(2))]
            test_planar_square(VT, d)
        end
    end
end


text:
# Script executed by Biowulf to run code

@everywhere include("/home/carsonc/StochasticGene/src/StochasticGene.jl")

include("/home/carsonc/StochasticGene/runfiles/scriptfunctions.jl")

@time fit_rna(parse(Int,ARGS[1]),ARGS[2],ARGS[3],parse(Int,ARGS[4]),parse(Float64,ARGS[5]),ARGS[6],ARGS[7],ARGS[8],ARGS[9],ARGS[10],parse(Int,ARGS[11]),parse(Bool,ARGS[12]),parse(Bool,ARGS[13]),parse(Bool,ARGS[14]))

# Arguments
# 1: nchains
# 2: gene
# 3: cell
# 4: cond
# 5: G
# 6: maxtime
# 7: infolder
# 8: resultfolder
# 9: datafolder
# 10: inlabel
# 11: label
# 12: nsets (number of rate sets)
# 13: runcycle (bool)
# 14: transient (bool)

max_stars_repo_name: nih-niddk-mbs/Transcription

max_issues_repo_name: nih-niddk-mbs/Transcription

max_forks_repo_name: nih-niddk-mbs/Transcription


text:
using DashBootstrapComponents, DashHtmlComponents

inputs = html_div([
    dbc_input(placeholder="Valid input...", valid=true, className="mb-3"),
    dbc_input(placeholder="Invalid input...", invalid=true),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
module DirectSweep
using BoundaryIntegralMethod
using PyCall
using LinearAlgebra

const linalg = PyNULL()
function __init__()
    copy!(linalg, pyimport_conda("scipy.linalg", "scipy"))
end

export directSweep

function directSweep(b::Boundary, n::Int, k0::Complex)
  @info "Doing initial precalculation for k=$k0"
  bpc0 = BoundaryIntegralMethod.doPrecalculation(b, n, k0)
  @info "Calculating A0"
  A0 = A(b, n, k0, bpc0)
  @info "Calculating A'"
  Apr = A′(b, n, k0, bpc0)
  try
    @info "Solving generalized eigenvalue problem"
    (vals, vecs) = linalg.eig(-A0, Apr)
    #TODO: Investigate performance of built in eigensolver
    # aEigen = eigen(-A0, Apr)
    # vals = aEigen.values
    # vecs = aEigen.vectors
    p = sortperm(abs.(vals))
    ks = k0.+vals[p]
    xs = vecs[:,p]'
    filterIdx = findall(k->imag(k) < 0 && imag(k) > -50,ks)
    ks = ks[filterIdx]
    xs = xs[filterIdx,:]
    (ks, xs)
  catch e
    @error "Got error attempting directSweep" exception=e
    rethrow(e)
  end
end
end


text:
export f

function f(x)
    function f_sub(y)
        norm(y)
    end
    2 * f_sub(x)
end


text:
module Kombinator
# TODO: introduce a graph problem in the abstract hierarchy? Graph-related functions could be defined on it.
# TODO: document the generic functions.
# TODO: refactor the package further to remove the old st_prim()-like names? Only solve(::CombinatorialInstance, ::CombinatorialAlgorithm)::CombinatorialSolution should be used. At least, don't export all these names.
# TODO: rename matching to "perfect matching" and introduce "imperfect matching".
# TODO: add a note to rather use LightGraph's algorithms instead of these in this package for anything serious.
# TODO: add a link to other packages like LightGraph for the corresponding algorithms, so that they seamlessly integrate with the others?
# TODO: how to add a link to functions like budgeted_msets_lp_select, budgeted_msets_lp_all? They can be quite useful (but do not bring much in terms of performance)

using DataStructures
using Hungarian
using JuMP
using LinearAlgebra
using LightGraphs
using Munkres
using Reexport

import Base: values, copy
import LightGraphs: src, dst
import JuMP: value, solve

include("helpers.jl")

include("interface/instance.jl")
include("interface/objective.jl")
include("interface/algorithm.jl")
include("interface/variants.jl")
include("interface/nonlinear.jl")

# include("BipartiteMatching/matching.jl")
# include("BipartiteMatching/matching_hungarian.jl")
# include("BipartiteMatching/matching_dp.jl")
# include("BipartiteMatching/matching_budgeted.jl")
# include("BipartiteMatching/matching_budgeted_dp.jl")
# include("BipartiteMatching/matching_budgeted_lagrangian.jl")

# Export all symbols. Code copied from JuMP.
symbols_to_exlude = [Symbol(@__MODULE__), :eval, :include, :symbols_to_exlude]

for sym in names(@__MODULE__, all=true)
    sym_string = string(sym)
    if sym in symbols_to_exlude || startswith(sym_string, "_")
        continue
    end
    if !(
        Base.isidentifier(sym) ||
        (startswith(sym_string, "@") && Base.isidentifier(sym_string[2:end]))
    )
        continue
    end
    @eval export $sym
end

# Reexport some symbols that are not included in the previous loop (imported
# from other modules).
export solve

# Include internal extensions.
include("UniformMatroid/UniformMatroid.jl")
@reexport using .UniformMatroid
include("SpanningTree/SpanningTree.jl")
@reexport using .SpanningTree
include("ElementaryPath/ElementaryPath.jl")
@reexport using .ElementaryPath
include("NLCOP/NLCOP.jl")
@reexport using .NLCOP

end


text:
# Kernels for unary array operations

using Knet: unary_ops

function cuda1src(f, j=f, ex="$f(xi)"; BLK=128, THR=128)
    sprint() do s
        for (T,F) in [("float","$(f)_32"),("double","$(f)_64")]
            print(s,
"""
__global__ void _$F(int n, $T *x, $T *y) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    $T xi = x[i];
    y[i] = $ex;
    i += blockDim.x * gridDim.x;
  }
}
#ifdef __cplusplus
extern "C" {
#endif
  void $F(int n, $T *x, $T *y) {
    _$F<<<$BLK,$THR>>>(n,x,y);
  }
#ifdef __cplusplus
}
#endif
""")
        end
    end
end

for a in unary_ops
    if !isa(a,Tuple); a=(a,); end
    print(cuda1src(a...))
end

# Kernels used by setindex! and getindex: fill, xfill, xcopy:

function cuda1fill(; BLK=256, THR=256)
    sprint() do s
        for (T,F) in [("float","32"),("double","64")]
            print(s,
"""
__global__ void _fill_$F(int n, $T x, $T *y) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    y[i] = x;
    i += blockDim.x * gridDim.x;
  }
}
#ifdef __cplusplus
extern "C" {
#endif
  void fill_$F(int n, $T x, $T *y) {
    _fill_$F<<<$BLK,$THR>>>(n,x,y);
  }
#ifdef __cplusplus
}
#endif
""")
        end
    end
end

print(cuda1fill())

function cuda1xfill(; BLK=256, THR=256)
    sprint() do s
        for (T,F) in [("float","32"),("double","64")]
            print(s,
"""
__global__ void _xfill_$F(int nrows, int ncols, $T x, $T *y, int incy) {
  int row, col, yidx;
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (1) {
    row = i % nrows;
    col = i / nrows;
    if (col >= ncols) break;
    yidx = row + col * incy;
    y[yidx] = x;
    i += blockDim.x * gridDim.x;
  }
}
#ifdef __cplusplus
extern "C" {
#endif
  void xfill_$F(int nrows, int ncols, $T x, $T *y, int incy) {
    _xfill_$F<<<$BLK,$THR>>>(nrows, ncols, x, y, incy);
  }
#ifdef __cplusplus
}
#endif
""")
        end
    end
end

print(cuda1xfill())

function cuda1xcopy(; BLK=256, THR=256)
"""
__global__ void _xcopy(int nrows, int ncols, const char *x, int incx, char *y, int incy) {
  int row, col, xidx, yidx;
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (1) {
    row = i % nrows;
    col = i / nrows;
    if (col >= ncols) break;
    xidx = row + col * incx;
    yidx = row + col * incy;
    y[yidx] = x[xidx];
    i += blockDim.x * gridDim.x;
  }
}
#ifdef __cplusplus
extern "C" {
#endif
  void xcopy(int nrows, int ncols, const void *x, int incx, void *y, int incy) {
    _xcopy<<<$BLK,$THR>>>(nrows,ncols,(char*)x,incx,(char*)y,incy);
  }
#ifdef __cplusplus
}
#endif
"""
end

print(cuda1xcopy())

permutedims3D_ops = [
#("permutedims3D_1_2_3","i","j","k"),#noop
("permutedims3D_1_3_2","i","k","j"),
("permutedims3D_2_1_3","j","i","k"),
("permutedims3D_2_3_1","k","i","j"),
("permutedims3D_3_1_2","j","k","i"),
("permutedims3D_3_2_1","k","j","i"),
]

function permutedims3Dsrc(f, i1, i2, i3; BLK=256, THR=256)
    sprint() do s
        for (T,F) in [("float","$(f)_32"),("double","$(f)_64")]
            print(s,
"""
__global__ void _$(F)_44($T* x, int dimx1, int dimx2, int dimx3, $T* y, int dimy1, int dimy2, int dimy3) {
  for (int v = threadIdx.x + blockIdx.x * blockDim.x; v < dimy1*dimy2*dimy3; v += blockDim.x * gridDim.x) {

		int i = v % dimy1;
		int j = ((v - i) / dimy1) % dimy2;
		int k = ((v - j * dimy1 - i) / (dimy1 * dimy2)) % dimy3;

		int srcIndex = $i1 + dimx1*$i2 + dimx1*dimx2*$i3;
		y[v] = x[srcIndex];
	}
}
#ifdef __cplusplus
extern "C" {
#endif
  void $(F)_44($T* x, int dimx1, int dimx2, int dimx3, $T* y, int dimy1, int dimy2, int dimy3) {
    _$(F)_44<<<$BLK,$THR>>>(x,dimx1,dimx2,dimx3,y,dimy1,dimy2,dimy3);
  }
#ifdef __cplusplus
}
#endif
""")
        end
    end
end

for a in permutedims3D_ops
    if !isa(a,Tuple); a=(a,); end
    print(permutedims3Dsrc(a...))
end


text:
# Single pendulums shall be solved numerically.
#
using OrdinaryDiffEq, Elliptic, Printf, DiffEqPhysics, Statistics

sol2q(sol) = [sol.u[i][j] for i in 1:length(sol.u), j in 1:length(sol.u[1])÷2]
sol2p(sol) = [sol.u[i][j] for i in 1:length(sol.u), j in length(sol.u[1])÷2+1:length(sol.u[1])]
sol2tqp(sol) = (sol.t, sol2q(sol), sol2p(sol))

# The exact solutions of single pendulums can be expressed by the Jacobian elliptic functions.
#
sn(u, k) = Jacobi.sn(u, k^2) # the Jacobian sn function

# Use PyPlot.
#
using PyPlot

colorlist = [
    "#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd",
    "#8c564b", "#e377c2", "#7f7f7f", "#bcbd22", "#17becf",
]
cc(k) = colorlist[mod1(k, length(colorlist))]

# plot the sulution of a Hamiltonian problem
#
function plotsol(sol::ODESolution)
    local t, q, p
    t, q, p = sol2tqp(sol)
    local d = size(q)[2]
    for j in 1:d
        j_str = d > 1 ? "[$j]" : ""
        plot(t, q[:,j], color=cc(2j-1), label="q$(j_str)", lw=1)
        plot(t, p[:,j], color=cc(2j),   label="p$(j_str)", lw=1, ls="--")
    end
    grid(ls=":")
    xlabel("t")
    legend()
end

# plot the solution of a Hamiltonian problem on the 2D phase space
#
function plotsol2(sol::ODESolution)
    local t, q, p
    t, q, p = sol2tqp(sol)
    local d = size(q)[2]
    for j in 1:d
        j_str = d > 1 ? "[$j]" : ""
        plot(q[:,j], p[:,j], color=cc(j), label="(q$(j_str),p$(j_str))", lw=1)
    end
    grid(ls=":")
    xlabel("q")
    ylabel("p")
    legend()
end

# plot the energy of a Hamiltonian problem
#
function plotenergy(H, sol::ODESolution)
    local t, q, p
    t, q, p = sol2tqp(sol)
    local energy = [H(q[i,:], p[i,:], nothing) for i in 1:size(q)[1]]
    plot(t, energy, label="energy", color="red", lw=1)
    grid(ls=":")
    xlabel("t")
    legend()
    local stdenergy_str = @sprintf("%.3e", std(energy))
    title("                    std(energy) = $stdenergy_str", fontsize=10)
end

# plot the numerical and exact solutions of a single pendulum
#
# Warning: Assume q(0) = 0, p(0) = 2k.   (for the sake of laziness)
#
function plotcomparison(k, sol::ODESolution)
    local t, q, p
    t, q, p = sol2tqp(sol)
    local y = sin.(q/2)
    local y_exact = k*sn.(t, k) # the exact solution

    plot(t, y,       label="numerical", lw=1)
    plot(t, y_exact, label="exact",     lw=1, ls="--")
    grid(ls=":")
    xlabel("t")
    ylabel("y = sin(q(t)/2)")
    legend()
    local error_str = @sprintf("%.3e", maximum(abs.(y - y_exact)))
    title("maximum(abs(numerical - exact)) = $error_str", fontsize=10)
end

# plot solution and energy
#
function plotsolenergy(H, integrator, Δt, sol::ODESolution)
    local integrator_str = replace("$integrator", r"^[^.]*\." => "")

    figure(figsize=(10,8))

    subplot2grid((21,20), ( 1, 0), rowspan=10, colspan=10)
    plotsol(sol)

    subplot2grid((21,20), ( 1,10), rowspan=10, colspan=10)
    plotsol2(sol)

    subplot2grid((21,20), (11, 0), rowspan=10, colspan=10)
    plotenergy(H, sol)

    suptitle("=====    $integrator_str,   Δt = $Δt    =====")
end

# Solve a single pendulum
#
function singlependulum(k, integrator, Δt; t0 = 0.0, t1 = 100.0)
    local H(p,q,params) = p[1]^2/2 - cos(q[1]) + 1
    local q0 = [0.0]
    local p0 = [2k]
    local prob = HamiltonianProblem(H, p0, q0, (t0, t1))

    local integrator_str = replace("$integrator", r"^[^.]*\." => "")
    @printf("%-25s", "$integrator_str:")
    sol = solve(prob, integrator, dt=Δt)
    @time local sol = solve(prob, integrator, dt=Δt)

    sleep(0.1)
    figure(figsize=(10,8))

    subplot2grid((21,20), ( 1, 0), rowspan=10, colspan=10)
    plotsol(sol)

    subplot2grid((21,20), ( 1,10), rowspan=10, colspan=10)
    plotsol2(sol)

    subplot2grid((21,20), (11, 0), rowspan=10, colspan=10)
    plotenergy(H, sol)

    subplot2grid((21,20), (11,10), rowspan=10, colspan=10)
    plotcomparison(k, sol)

    suptitle("=====    $integrator_str,   Δt = $Δt    =====")
end

# Single pendulum

k = rand()
integrator = VelocityVerlet()
Δt = 0.1
singlependulum(k, integrator, Δt, t0=-20.0, t1=20.0)

# Two single pendulums

H(q,p,param) = sum(p.^2/2 .- cos.(q) .+ 1)
q0 = pi*rand(2)
p0 = zeros(2)
t0, t1 = -20.0, 20.0
prob = HamiltonianProblem(H, q0, p0, (t0, t1))

integrator = VelocityVerlet()
Δt = 0.1
sol = solve(prob, integrator, dt=Δt)
@time sol = solve(prob, integrator, dt=Δt)

sleep(0.1)
plotsolenergy(H, integrator, Δt, sol)

SymplecticIntegrators = [
    SymplecticEuler(),
    VelocityVerlet(),
    VerletLeapfrog(),
    PseudoVerletLeapfrog(),
    McAte2(),
    Ruth3(),
    McAte3(),
    CandyRoz4(),
    McAte4(),
    CalvoSanz4(),
    McAte42(),
    McAte5(),
    Yoshida6(),
    KahanLi6(),
    McAte8(),
    KahanLi8(),
    SofSpa10(),
]

k = 0.999
Δt = 0.1
for integrator in SymplecticIntegrators
    singlependulum(k, integrator, Δt)
end

k = 0.999
Δt = 0.01
for integrator in SymplecticIntegrators[1:4]
    singlependulum(k, integrator, Δt)
end

k = 0.999
Δt = 0.001
singlependulum(k, SymplecticEuler(), Δt)

k = 0.999
Δt = 0.0001
singlependulum(k, SymplecticEuler(), Δt)

using DiffEqBenchmarks
DiffEqBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
module BlockDiagonals

using Base: @propagate_inbounds
using FillArrays: Zeros
using LinearAlgebra

export BlockDiagonal, blocks
export blocksize, blocksizes, nblocks

include("blockdiagonal.jl")
include("base_maths.jl")
include("linalg.jl")

end  # end module


text:
#=--------------------------------------------------------------------

ransac - Functions for robust estimation using RANSAC

Copyright (c) 2016 Peter Kovesi
pk@peterkovesi.com

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

The Software is provided "as is", without warranty of any kind.

PK February 2016

---------------------------------------------------------------------=#

export ransac
export ransacfithomography, ransacfitfundmatrix, ransacfitaffinefundmatrix
export ransacfitplane, ransacfitline
export iscolinear, fitline2d, fitline3d, fitplane

#---------------------------------------------------------------------
"""
ransac - Robustly fits a model to data with the RANSAC algorithm
```
Usage:

(M, inliers) = ransac(x, fittingfn, distfn, degenfn s, t, feedback,
                      maxDataTrials, maxTrials)

Arguments:
    x         - Data sets to which we are seeking to fit a model M
                It is assumed that x is of size [d x Npts]
                where d is the dimensionality of the data and Npts is
                the number of data points.

    fittingfn - Reference to a function that fits a model to s
                data from x.  It is assumed that the function is of the
                form:
                   M = fittingfn(x)
                Note it is possible that the fitting function can return
                multiple models (for example up to 3 fundamental matrices
                can be fitted to 7 pairs of matched points).  In this
                case it is assumed that the fitting function returns an
                array of models.
                If this function cannot fit a model it should return M as
                an empty matrix.

    distfn    - Reference to a function that evaluates the
                distances from the model to data x.
                It is assumed that the function is of the form:
                   (inliers, M) = distfn(M, x, t)
                This function must evaluate the distances between points
                and the model returning the indices of elements in x that
                are inliers, that is, the points that are within distance
                't' of the model.  Additionally, if M is an array of
                possible models distfn() will return the model that has the
                most inliers.  If there is only one model this function
                must still copy the model to the output.  After this call M
                will be a single object representing only one model.

    degenfn   - Reference to a function that determines whether a
                set of datapoints will produce a degenerate model.
                This is used to discard random samples that do not
                result in useful models.
                It is assumed that degenfn is a boolean function of
                the form:
                   r = degenfn(x)
                It may be that you cannot devise a test for degeneracy in
                which case you should write a dummy function that always
                returns a value of true and rely on fittingfn() to return
                an empty model should the data set be degenerate.

    s         - The minimum number of samples from x required by
                fittingfn to fit a model.

    t         - The distance threshold between a data point and the model
                used to decide whether the point is an inlier or not.

    feedback  - An optional boolean flag. If set to true the trial count and the
                estimated total number of trials required is printed out at
                each step.  Defaults to false.

maxDataTrials - Maximum number of attempts to select a non-degenerate
                data set. This parameter is optional and defaults to 100.

    maxTrials - Maximum number of iterations. This parameter is optional and
                defaults to 1000.

    p         - Desired probability of choosing at least one sample
                free from outliers, defaults to 0.99

Returns:
    M         - The model having the greatest number of inliers.
    inliers   - An array of indices of the elements of x that were
                the inliers for the best model.
```
For an example of the use of this function see ransacfithomography() or
ransacfitplane()
"""
#=
References:
   M.A. Fishler and  R.C. Boles. "Random sample concensus: A paradigm
   for model fitting with applications to image analysis and automated
   cartography". Comm. Assoc. Comp, Mach., Vol 24, No 6, pp 381-395, 1981

   Richard Hartley and Andrew Zisserman. "Multiple View Geometry in
   Computer Vision". pp 101-113. Cambridge University Press, 2001

May      2003 - Original version in MATLAB
February 2004 - Tidied up.
August   2005 - Specification of distfn changed to allow model fitter to
                return multiple models from which the best must be selected
Sept     2006 - Random selection of data points changed to ensure duplicate
                points are not selected.
December 2008 - Octave compatibility mods
June     2009 - Argument 'MaxTrials' corrected to 'maxTrials'!
February 2016 - Ported to Julia
=#

function ransac(x, fittingfn, distfn, degenfn, s, t, feedback = false,
                maxDataTrials = 1000, maxTrials = 1000, p = 0.99)

    (rows, npts) = size(x)

    if npts < s
        error("Not enough data points to form a model")
    end

    bestM = NaN      # Sentinel value allowing detection of solution failure.
    trialcount = 0
    bestscore =  0
    bestinliers = []  # Ensure global to function
    bestM = []
    N = 1            # Dummy initialisation for number of trials.

    while N > trialcount

        # Select at random s datapoints to form a trial model, M.
        # In selecting these points we have to check that they are not in
        # a degenerate configuration.
        degenerate = true
        count = 1
        while degenerate
            # Generate s random indicies in the range 1..npts
            ind = randperm(npts)[1:s]

            # Test that these points are not a degenerate configuration.
            degenerate = degenfn(x[:,ind])

            if !degenerate
                # Fit model to this random selection of data points.
                # Note that M may represent a set of models that fit the data in
                # this case M will be a cell array of models
                M = fittingfn(x[:,ind])

                # Depending on your problem it might be that the only way you
                # can determine whether a data set is degenerate or not is to
                # try to fit a model and see if it succeeds.  If it fails we
                # reset degenerate to true.
                if isempty(M)
                    degenerate = true
                end
            end

            # Safeguard against being stuck in this loop forever
            count += 1
            if count > maxDataTrials
                warn("Unable to select a nondegenerate data set")
                break
            end
        end

        # Once we are out here we should have some kind of model...
        # Evaluate distances between points and model returning the indices
        # of elements in x that are inliers.  Additionally, if M is a cell
        # array of possible models 'distfn' will return the model that has
        # the most inliers.  After this call M will be a non-cell object
        # representing only one model.
        (inliers, M) = distfn(M, x, t)

        # Find the number of inliers to this model.
        ninliers = length(inliers)

        if ninliers > bestscore    # Largest set of inliers so far...
            bestscore = ninliers   # Record data for this model
            bestinliers = copy(inliers)
            bestM = copy(M)

            # Update estimate of N, the number of trials to ensure we pick,
            # with probability p, a data set with no outliers.
            fracinliers =  ninliers/npts
            pNoOutliers = 1 -  fracinliers^s
            pNoOutliers = max(eps(), pNoOutliers)  # Avoid division by -Inf
            pNoOutliers = min(1-eps(), pNoOutliers)# Avoid division by 0.
            N = log(1-p)/log(pNoOutliers)
        end

        trialcount += 1
        if feedback
            @printf("trial %d out of %d         \r",trialcount, ceil(N))
        end

        # Safeguard against being stuck in this loop forever
        if trialcount > maxTrials
            @printf("RANSAC reached the maximum number of %d trials", maxTrials)
            break
        end
    end

    if feedback; @printf("\n"); end

    if !any(isnan(bestM))   # We got a solution
        M = bestM
        inliers = bestinliers
    else
        M = []
        inliers = []
        error("RANSAC was unable to find a useful solution")
    end

    return M, inliers
end

#-----------------------------------------------------------------------
"""
ransacfithomography - Fits 2D homography using RANSAC
```
Usage:   (H, inliers) = ransacfithomography(x1, x2, t)

Arguments:
         x1  - 2xN or 3xN set of homogeneous points.  If the data is
               2xN it is assumed the homogeneous scale factor is 1.
         x2  - 2xN or 3xN set of homogeneous points such that x1<->x2.
         t   - The distance threshold between data point and the model
               used to decide whether a point is an inlier or not.
               Note that point coordinates are normalised to that their
               mean distance from the origin is sqrt(2).  The value of
               t should be set relative to this, say in the range
               0.001 - 0.01
```
Note that it is assumed that the matching of x1 and x2 are putative and it
is expected that a percentage of matches will be wrong.
```
Returns:
         H       - The 3x3 homography such that x2 = H*x1.
         inliers - An array of indices of the elements of x1, x2 that were
                   the inliers for the best model.
```
See Also: ransac(), homography2d(), homography1d()
"""

#=
February 2004 - original version
July     2004 - error in denormalising corrected (thanks to Andrew Stein)
August   2005 - homogdist2d modified to fit new ransac specification.
=#

function ransacfithomography(x1i, x2i, t::Real)

    # Define some subfunctions

    #----------------------------------------------------------------------
    # Function to evaluate the symmetric transfer error of a homography with
    # respect to a set of matched points as needed by RANSAC.

    function homogdist2d(H, x, t)

        x1 = x[1:3,:]   # Extract x1 and x2 from x
        x2 = x[4:6,:]

        # Calculate, in both directions, the transfered points
        Hx1    = H*x1
        invHx2 = H\x2

        # Normalise so that the homogeneous scale parameter for all coordinates
        # is 1.
        x1     = hnormalise(x1)
        x2     = hnormalise(x2)
        Hx1    = hnormalise(Hx1)
        invHx2 = hnormalise(invHx2)

        d2 = sum((x1-invHx2).^2, 1)  + sum((x2-Hx1).^2, 1)
        inliers = find(abs(d2) .< t)
        return inliers, H
    end

    #----------------------------------------------------------------------
    # Function to determine if a set of 4 pairs of matched  points give rise
    # to a degeneracy in the calculation of a homography as needed by RANSAC.
    # This involves testing whether any 3 of the 4 points in each set is
    # colinear.

    function isdegenerate(x)

        (rows, npts) = size(x)
        assert(npts == 4)

        x1 = x[1:3,:]    # Extract x1 and x2 from x
        x2 = x[4:6,:]

        r =
        iscolinear(x1[:,1],x1[:,2],x1[:,3]) ||
        iscolinear(x1[:,1],x1[:,2],x1[:,4]) ||
        iscolinear(x1[:,1],x1[:,3],x1[:,4]) ||
        iscolinear(x1[:,2],x1[:,3],x1[:,4]) ||
        iscolinear(x2[:,1],x2[:,2],x2[:,3]) ||
        iscolinear(x2[:,1],x2[:,2],x2[:,4]) ||
        iscolinear(x2[:,1],x2[:,3],x2[:,4]) ||
        iscolinear(x2[:,2],x2[:,3],x2[:,4])

        return r
    end
    #----------------------------------------------------------

    # The main function

    if size(x1i) != size(x2i)
        error("Data sets x1 and x2 must have the same dimensions")
    end

    (rows,npts) = size(x1i)
    if !(rows == 2 || rows == 3)
        error("x1 and x2 must have 2 or 3 rows")
    end

    if npts < 4
        error("Must have at least 4 points to fit homography")
    end

    if rows == 2    # Pad data with homogeneous scale factor of 1
        x1 = [x1i; ones(1,npts)]
        x2 = [x2i; ones(1,npts)]
    else
        x1 = copy(x1i)
        x2 = copy(x2i)
    end

    # Normalise each set of points so that the origin is at centroid and
    # mean distance from origin is sqrt(2).  normalise2dpts also ensures the
    # scale parameter is 1.  Note that 'homography2d' will also call
    # 'normalise2dpts' but the code in 'ransac' that calls the distance
    # function will not - so it is best that we normalise beforehand.
    (x1, T1) = normalise2dpts(x1)
    (x2, T2) = normalise2dpts(x2)

    s = 4  # Minimum No of points needed to fit a homography.

    fittingfn = homography2d
    distfn    = homogdist2d
    degenfn   = isdegenerate
    # x1 and x2 are 'stacked' to create a 6xN array for ransac
    (H, inliers) = ransac([x1; x2], fittingfn, distfn, degenfn, s, t)

    # Now do a final least squares fit on the data points considered to
    # be inliers.
    H = homography2d(x1[:,inliers], x2[:,inliers])

    # Denormalise
    H = T2\H*T1

    return H, inliers

end # of ransacfithomography

#-----------------------------------------------------------------------
"""
ransacfitfundmatrix - Fits fundamental matrix using RANSAC
```
Usage:   (F, inliers) = ransacfitfundmatrix(x1, x2, t)

Arguments:
         x1  - 2xN or 3xN set of homogeneous points.  If the data is
               2xN it is assumed the homogeneous scale factor is 1.
         x2  - 2xN or 3xN set of homogeneous points such that x1<->x2.
         t   - The distance threshold between data point and the model
               used to decide whether a point is an inlier or not.
               Note that point coordinates are normalised to that their
               mean distance from the origin is sqrt(2).  The value of
               t should be set relative to this, say in the range
               0.001 - 0.01
```
Note that it is assumed that the matching of x1 and x2 are putative and it
is expected that a percentage of matches will be wrong.
```
Returns:
         F       - The 3x3 fundamental matrix such that x2'Fx1 = 0.
         inliers - An array of indices of the elements of x1, x2 that were
                   the inliers for the best model.
```
See also: ransac(), fundmatrix()
"""

# February 2004  Original version
# August   2005  Distance error function changed to match changes in RANSAC

function ransacfitfundmatrix(x1i, x2i, t::Real, feedback::Bool = false)

    # Define some subfunctions

    #--------------------------------------------------------------------------
    # Function to evaluate the first order approximation of the geometric error
    # (Sampson distance) of the fit of a fundamental matrix with respect to a
    # set of matched points as needed by RANSAC.  See: Hartley and Zisserman,
    # 'Multiple View Geometry in Computer Vision', page 270.
    #
    # Note that this code allows for F being an array of fundamental matrices of
    # which we have to pick the best one. (A 7 point solution can return up to 3
    # solutions)

    function funddist(F, x, t)

        npts = size(x, 2)
        x1 = x[1:3,:]    # Extract x1 and x2 from x
        x2 = x[4:6,:]

        if isa(F, Array{Any})  # We have several solutions each of which must be tested
	    nF = length(F)   # Number of solutions to test
	    bestF = F[1]     # Initial allocation of best solution
	    ninliers = 0     # Number of inliers

	    for k = 1:nF
	        x2tFx1 = zeros(npts)
	        for n = 1:npts
		    x2tFx1[n] = (x2[:,n:n]'*F[k]*x1[:,n:n])[1]
	        end

	        Fx1 = F[k]*x1
	        Ftx2 = F[k]'*x2

	        # Evaluate distances
	        d =  x2tFx1.^2 ./
		vec(Fx1[1,:].^2 + Fx1[2,:].^2 + Ftx2[1,:].^2 + Ftx2[2,:].^2)

	        inliers = find(abs(d) .< t)     # Indices of inlying points

	        if length(inliers) > ninliers   # Record best solution
		    ninliers = length(inliers)
		    bestF = copy(F[k])
		    bestInliers = copy(inliers)
	        end
	    end

        else     # We just have one solution
	    x2tFx1 = zeros(npts)
	    for n = 1:npts
	        x2tFx1[n] = (x2[:,n:n]'*F*x1[:,n:n])[1] # Ugly thanks to v0.5
	    end

	    Fx1 = F*x1
	    Ftx2 = F'*x2

	    # Evaluate distances (The call to vec() needed for 0.4)
	    d =  x2tFx1.^2 ./
	    vec(Fx1[1,:].^2 + Fx1[2,:].^2 + Ftx2[1,:].^2 + Ftx2[2,:].^2)

	    bestInliers = find(abs(d) .< t)     # Indices of inlying points
	    bestF = F                          # Copy F directly to bestF
        end

        return bestInliers, bestF
    end

    #----------------------------------------------------------------------
    # (Degenerate!) function to determine if a set of matched points will result
    # in a degeneracy in the calculation of a fundamental matrix as needed by
    # RANSAC.  This function assumes this cannot happen...

    function isdegenerate(x)
        return false
    end

    #----------------------------------------------------------------------
    # The main function

    if size(x1i) != size(x2i)
        error("Data sets x1 and x2 must have the same dimension")
    end

    (rows, npts) = size(x1i)
    if !(rows==2 || rows==3)
        error("x1 and x2 must have 2 or 3 rows")
    end

    if rows == 2    # Pad data with homogeneous scale factor of 1
        x1 = makehomogeneous(x1i)
        x2 = makehomogeneous(x2i)
    else
        x1 = copy(x1i)
        x2 = copy(x2i)
    end

    # Normalise each set of points so that the origin is at centroid and
    # mean distance from origin is sqrt(2).  normalise2dpts also ensures the
    # scale parameter is 1.  Note that fundmatrix() will also call
    # normalise2dpts() but the code in ransac() that calls the distance
    # function will not - so it is best that we normalise beforehand.
    (x1, T1) = normalise2dpts(x1)
    (x2, T2) = normalise2dpts(x2)

    s = 8  # Number of points needed to fit a fundamental matrix. Note that
           # only 7 are needed but the function fundmatrix() only
           # implements the 8-point solution.

    fittingfn = fundmatrix
    distfn    = funddist
    degenfn   = isdegenerate
    # x1 and x2 are 'stacked' to create a 6xN array for ransac
    (F, inliers) = ransac([x1; x2], fittingfn, distfn, degenfn, s, t, feedback)

    # Now do a final least squares fit on the data points considered to
    # be inliers.
    F = fundmatrix(x1[:,inliers], x2[:,inliers])

    # Denormalise
    F = T2'*F*T1

    return F, inliers

end # of ransacfitfundmatrix()

#-----------------------------------------------------------------------
"""
ransacfitaffinefundmatrix - Fits affine fundamental matrix using RANSAC
```
Usage:   (F, inliers) = ransacfitaffinefundmatrix(x1, x2, t)

Arguments:
         x1  - 2xN or 3xN set of homogeneous points.  If the data is
               2xN it is assumed the homogeneous scale factor is 1.
         x2  - 2xN or 3xN set of homogeneous points such that x1<->x2.
         t   - The distance threshold between data point and the model
               used to decide whether a point is an inlier or not.
               Note that point coordinates are normalised to that their
               mean distance from the origin is sqrt(2).  The value of
               t should be set relative to this, say in the range
               0.001 - 0.01
```
Note that it is assumed that the matching of x1 and x2 are putative and it
is expected that a percentage of matches will be wrong.
```
Returns:
         F       - The 3x3 fundamental matrix such that x2'Fx1 = 0.
         inliers - An array of indices of the elements of x1, x2 that were
                   the inliers for the best model.
```
See also: ransac(), fundmatrix(), affinefundmatrix()
"""

# February 2004  Original version
# August   2005  Distance error function changed to match changes in RANSAC

function ransacfitaffinefundmatrix(x1i, x2i, t::Real, feedback::Bool = false)

    # Define some subfunctions
    #--------------------------------------------------------------------------
    # Function to evaluate the first order approximation of the geometric error
    # (Sampson distance) of the fit of a fundamental matrix with respect to a
    # set of matched points as needed by RANSAC.  See: Hartley and Zisserman,
    # 'Multiple View Geometry in Computer Vision', page 270.

    function funddist(F, x, t)

        ( _ , npts) = size(x)
        x1 = x[1:3,:]    # Extract x1 and x2 from x
        x2 = x[4:6,:]

        x2tFx1 = zeros(npts)
        for n = 1:npts
	    x2tFx1[n] = (x2[:,n]'*F*x1[:,n])[1]
        end

        Fx1 = F*x1
        Ftx2 = F'*x2

        # Evaluate distances
        d =  x2tFx1.^2 ./
	vec(Fx1[1,:].^2 + Fx1[2,:].^2 + Ftx2[1,:].^2 + Ftx2[2,:].^2)

        inliers = find(abs(d) .< t)     # Indices of inlying points
        return inliers, F
    end

    #----------------------------------------------------------------------
    # (Degenerate!) function to determine if a set of matched points will result
    # in a degeneracy in the calculation of a fundamental matrix as needed by
    # RANSAC.  This function assumes this cannot happen...

    function isdegenerate(x)
        return false
    end
    #----------------------------------------------------------------------
    # The main function

    if size(x1i) != size(x2i)
        error("Data sets x1 and x2 must have the same dimension")
    end

    (rows, npts) = size(x1i)
    if !(rows ==2 || rows ==3)
        error("x1 and x2 must have 2 or 3 rows")
    end

    if rows == 2    # Pad data with homogeneous scale factor of 1
        x1 = makehomogeneous(x1i)
        x2 = makehomogeneous(x2i)
    else
        x1 = copy(x1i)
        x2 = copy(x2i)
    end

    # Normalise each set of points so that the origin is at centroid and
    # mean distance from origin is sqrt(2).  normalise2dpts also ensures the
    # scale parameter is 1.  Note that 'fundmatrix' will also call
    # 'normalise2dpts' but the code in 'ransac' that calls the distance
    # function will not - so it is best that we normalise beforehand.
    (x1, T1) = normalise2dpts(x1)
    (x2, T2) = normalise2dpts(x2)

    s = 4  # Number of points needed to fit an affine fundamental matrix.

    fittingfn = affinefundmatrix
    distfn    = funddist
    degenfn   = isdegenerate
    # x1 and x2 are 'stacked' to create a 6xN array for ransac
    (F, inliers) = ransac([x1; x2], fittingfn, distfn, degenfn, s, t, feedback)

    # Now do a final least squares fit on the data points considered to
    # be inliers.
    F = affinefundmatrix(x1[:,inliers], x2[:,inliers])

    # Denormalise
    F = T2'*F*T1

    return F, inliers

end # of ransacfitaffinefundmatrix

#-----------------------------------------------------------------------
"""
ransacfitplane - Fits plane to 3D array of points using RANSAC
```
Usage  (B, P, inliers) = ransacfitplane(XYZ, t, feedback)

This function uses the RANSAC algorithm to robustly fit a plane
to a set of 3D data points.

Arguments:
         XYZ - 3xNpts array of xyz coordinates to fit plane to.
         t   - The distance threshold between data point and the plane
               used to decide whether a point is an inlier or not.
         feedback - Optional flag 0 or 1 to turn on RANSAC feedback
                    information.

Returns:
          B - 4x1 array of plane coefficients in the form
              b[1]*X + b[2]*Y +b[3]*Z + b[4] = 0
              The magnitude of B is 1.
              This plane is obtained by a least squares fit to all the
              points that were considered to be inliers, hence this
              plane will be slightly different to that defined by P below.
          P - The three points in the data set that were found to
              define a plane having the most number of inliers.
              The three columns of P defining the three points
    inliers - The indices of the points that were considered
              inliers to the fitted plane.
```
See also:  ransac(), fitplane()
"""
#=
June 2003 - Original version.
Feb  2004 - Modified to use separate ransac function
Aug  2005 - planeptdist modified to fit new ransac specification
Dec  2008 - Much faster distance calculation in planeptdist (thanks to
            Alastair Harrison)
=#

function ransacfitplane(XYZ, t::Real, feedback::Bool = false)

    # Define some sub functions
    #------------------------------------------------------------------------
    # Function to define a plane given 3 data points as required by
    # RANSAC. In our case we use the 3 points directly to define the plane.
    function defineplane(X)
        return X
    end

    #------------------------------------------------------------------------
    # Function to calculate distances between a plane and a an array of points.
    # The plane is defined by a 3x3 matrix, P.  The three columns of P defining
    # three points that are within the plane.
    function planeptdist(P, X, t)

        n = cross(P[:,2]-P[:,1], P[:,3]-P[:,1]) # Plane normal.
        n = n/norm(n)                           # Make it a unit vector.

        npts = size(X,2)
        d = zeros(npts)   # d will be an array of distance values.

        # The following loop builds up the dot product between a vector from P[:,1]
        # to every X[:,i] with the unit plane normal.  This will be the
        # perpendicular distance from the plane for each point
        for i=1:3
	    d = d + vec(X[i,:]-P[i,1])*n[i]
        end

        inliers = find(abs(d) .< t)

        return inliers, P
    end
    #------------------------------------------------------------------------
    # Function to determine whether a set of 3 points are in a degenerate
    # configuration for fitting a plane as required by RANSAC.  In this case
    # they are degenerate if they are colinear.

    function isdegenerate(X)

        # The three columns of X are the coords of the 3 points.
        r = iscolinear(X[:,1],X[:,2],X[:,3])
    end

    #------------------------------------------------------------------------
    # The main function

    (rows, npts) = size(XYZ)

    if rows !=3
        error("Data is not 3D")
    end

    if npts < 3
        error("Too few points to fit plane")
    end

    s = 3  # Minimum No of points needed to fit a plane.

    fittingfn = defineplane
    distfn    = planeptdist
    degenfn   = isdegenerate

    (P, inliers) = ransac(XYZ, fittingfn, distfn, degenfn, s, t, feedback)

    # Perform least squares fit to the inlying points
    B = fitplane(XYZ[:,inliers])

    return B, P, inliers
end # of ransacfitplane

#-----------------------------------------------------------------------
"""
ransacfitline - Fits line to 3D array of points using RANSAC
```
Usage  (V, L, inliers) = ransacfitline(XYZ, t, feedback)

This function uses the RANSAC algorithm to robustly fit a line to a
set of 3D data points or to a set of 2D homogeneous points with scale
value of 1

Arguments:
         XYZ - 3xNpts array of xyz coordinates to fit line to.
         t   - The distance threshold between data point and the line
               used to decide whether a point is an inlier or not.
    feedback - Optional boolean flag to turn on RANSAC feedback
               information.

Returns:.
          V - Line obtained by a simple fitting on the points that
              are considered inliers.  The line goes through the
              calculated mean of the inlier points, and is parallel to
              the principal eigenvector.  The line is scaled by the
              square root of the largest eigenvalue.
              This line is returned as a nx2 matrix.  The first column
              is the beginning point, the second column is the end point
              of the line.
          L - The two points in the data set that were found to
              define a line having the most number of inliers.
              The two columns of L defining the two points.
    inliers - The indices of the points that were considered
              inliers to the fitted line.
```
See also:  ransac(), fitplane(), ransacfitplane()
"""
# Aug  2006 - created ransacfitline from ransacfitplane
#             author: Felix Duvallet

function ransacfitline(XYZ, t::Real, feedback::Bool = false)

    # First define some internal functions
    # ------------------------------------------------------------------------
    # Function to define a line given 2 data points as required by
    # RANSAC.  X is a 2 column matrix where each column defines a
    # point. Here we use the 2 points directly to define the line

    function defineline(X)
        return X
    end

    #------------------------------------------------------------------------
    # Function to calculate distances between a line and an array of points.
    # The line is defined by a 3x2 matrix, L.  The two columns of L defining
    # two points that are the endpoints of the line.
    #
    # A line can be defined with two points as:
    #        lambda*p1 + (1-lambda)*p2
    # Then, the distance between the line and another point (p3) is:
    #        norm( lambda*p1 + (1-lambda)*p2 - p3 )
    # where
    #                  (p2-p1).(p2-p3)
    #        lambda =  ---------------
    #                  (p1-p2).(p1-p2)
    #
    # lambda can be found by taking the derivative of:
    #      (lambda*p1 + (1-lambda)*p2 - p3)*(lambda*p1 + (1-lambda)*p2 - p3)
    # with respect to lambda and setting it equal to zero

    function lineptdist(L, X, t)

        p1 = L[:,1]
        p2 = L[:,2]

        npts = size(X,2)
        d = zeros(npts)

        for i = 1:npts
            p3 = X[:,i]

            lambda = dot((p2 - p1), (p2-p3)) / dot( (p1-p2), (p1-p2) )

            d[i] = norm(lambda*p1 + (1-lambda)*p2 - p3)
        end

        inliers = find(abs(d) .< t)
        return inliers, L
    end

    #------------------------------------------------------------------------
    # Function to determine whether a set of 2 points are in a degenerate
    # configuration for fitting a line as required by RANSAC.
    # In this case two points are degenerate if they are the same point
    # or if they are exceedingly close together.

    function isdegenerate(X)
        # Find the norm of the difference of the two points.
        # If this is zero we are degenerate
        return norm(X[:,1] - X[:,2]) < eps()
    end

    #------------------------------------------------------------------------
    # Now the main function ...

    (rows, npts) = size(XYZ)

    if rows !=3
        error("Data is not 3D")
    end

    if npts < 2
        error("Too few points to fit line")
    end

    s = 2  # Minimum No of points needed to fit a line.

    fittingfn = defineline
    distfn    = lineptdist
    degenfn   = isdegenerate

    (L, inliers) = ransac(XYZ, fittingfn, distfn, degenfn, s, t, feedback)

    # Find the line going through the mean, parallel to the major
    # eigenvector
    V = fitline3d(XYZ[:, inliers])

    return V, L, inliers

end # of ransacfitline

#-----------------------------------------------------------------------
"""
iscolinear - Are 3 points colinear
```
Usage:  r = iscolinear(p1, p2, p3, flag)

Arguments:
       p1, p2, p3 - Points in 2D or 3D.
       homog      - An optional boolean flag indicating that p1, p2, p3
                    are homogneeous coordinates with arbitrary scale.
                    The default is false, that is it is assumed that the
                    points are inhomogeneous, or that they are homogeneous
                    ewith qual scale.

Returns:
       r  -  true if points are co-linear, false otherwise
```
"""
# February 2004
# January  2005 - modified to allow for homogeneous points of arbitrary
#                 scale (thanks to Michael Kirchhof)

function iscolinear(p1i::Vector, p2i::Vector, p3i::Vector, homog::Bool = false)

    if (length(p1i) != length(p2i))  || (length(p1i) != length(p3i)) || !(length(p1i) == 2 || length(p1i) == 3)
        error("Points must have the same dimension of 2 or 3")
    end

    # If data is 2D, assume they are 2D inhomogeneous coords. Make them
    # homogeneous with scale 1.
    if length(p1i) == 2
        p1 = [p1i; 1]; p2 =[p2i; 1]; p3 = [p3i; 1]
    else
        p1 = p1i; p2 = p2i; p3 = p3i
    end

    if homog
	# Apply test that allows for homogeneous coords with arbitrary
        # scale.  p1 X p2 generates a normal vector to plane defined by
        # origin, p1 and p2.  If the dot product of this normal with p3
        # is zero then p3 also lies in the plane, hence co-linear.
	r = abs(dot(cross(p1, p2),p3)) < eps()
    else
	# Assume inhomogeneous coords, or homogeneous coords with equal
        # scale.
	r = norm(cross(p2-p1, p3-p1)) < eps()
    end

    return r
end

#-----------------------------------------------------------------------
"""
fitline2d - Least squares fit of a line to a set of 2D points
```
Usage:   (C, dist) = fitline2d(XY)

Where:   XY  - 2xNpts array of xy coordinates to fit line to data of
               the form
               [x1 x2 x3 ... xN
                y1 y2 y3 ... yN]

               XY can also be a 3xNpts array of homogeneous coordinates.

Returns: C    - 3x1 array of line coefficients in the form
                c[1]*X + c[2]*Y + c[3] = 0
         dist - Array of distances from the fitted line to the supplied
                data points.  Note that dist is only calculated if the
                function is called with two output arguments.
```
The magnitude of C is scaled so that line equation corresponds to
```
   sin(theta)*X + (-cos(theta))*Y + rho = 0
```
where theta is the angle between the line and the x axis and rho is the
perpendicular distance from the origin to the line.  Rescaling the
coefficients in this manner allows the perpendicular distance from any
point (x,y) to the line to be simply calculated as
```
   r = abs(c[1]*X + c[2]*Y + c[3])
```
If you want to convert this line representation to the classical form
```
     Y = a*X + b
use
  a = -c[1]/c[2]
  b = -c[3]/c[2]
```
Note, however, that this assumes c[2] is not zero
"""
# June      2003 - Original version
# September 2004 - Rescaling to allow simple distance calculation.
# November  2008 - Normalising of coordinates added to condition the solution.

function fitline2d(XYi)

    (rows,npts) = size(XYi)

    if npts < 2
        error("Too few points to fit line")
    end

    if rows == 2    # Add homogeneous scale coordinate of 1
        XY = makehomogneous(XYi)
    else
        XY = copy(XYi)
    end

    if npts == 2    # Pad XY with a third column of zeros
        XY = [XY zeros(3,1)]
    end

    # Normalise points so that centroid is at origin and mean distance from
    # origin is sqrt(2).  This conditions the equations
    (XYn, T) = normalise2dpts(XY)

    # Set up constraint equations of the form  XYn'*C = 0,
    # where C is a column vector of the line coefficients
    # in the form   c(1)*X + c(2)*Y + c(3) = 0.

    (u, d, v) = svd(XYn')   # Singular value decomposition.
    C = v[:,3]              # Solution is last column of v.

    # Denormalise the solution
    C = T'*C

    # Rescale coefficients so that line equation corresponds to
    #   sin(theta)*X + (-cos(theta))*Y + rho = 0
    # so that the perpendicular distance from any point [x,y] to the line
    # to be simply calculated as
    #   r = abs(c[1]*X + c[2]*Y + c[3])
    theta = atan2(C[1], -C[2])

    # Find the scaling (but avoid dividing by zero)
    if abs(sin(theta)) > abs(cos(theta))
        k = C[1]/sin(theta)
    else
        k = -C[2]/cos(theta)
    end

    C = C/k

    # Calculate the distances from the fitted line to the supplied
    # data points
    dist = abs(C[1]*XY[1,:] + C[2]*XY[2,:] + C[3])

    return C, dist
end

#-----------------------------------------------------------------------
"""
fitline3d - Fits a line to a set of 3D points
```
Usage:   L = fitline3d(XYZ)

Where: XYZ - 3xNpts array of XYZ coordinates
              [x1 x2 x3 ... xN;
               y1 y2 y3 ... yN;
               z1 z2 z3 ... zN]

Returns: L - 3x2 matrix consisting of the two endpoints of the line
             that fits the points.  The line is centered about the
             mean of the points, and extends in the directions of the
             principal eigenvectors, with scale determined by the
             eigenvalues.
```
"""
# Author: Felix Duvallet (CMU)
# August 2006

function fitline3d(XYZ)

    # Since the covariance matrix should be 3x3 (not NxN), need
    # to take the transpose of the points.

    # find mean of the points
    mu = mean(XYZ', 1)

    # covariance matrix
    C = cov(XYZ')

    # get the eigenvalues and eigenvectors
    (D, V) = eig(C)

    # largest eigenvector is in the last column
    col = size(V, 2)  # get the number of columns

    # get the last eigenvector column and the last eigenvalue
    eVec = V[:, col]
    eVal = D[col]

    L = zeros(3,2)
    # start point - center about mean and scale eVector by eValue
    L[:, 1] = mu' - sqrt(eVal)*eVec
    L[:, 2] = mu' + sqrt(eVal)*eVec      # end point

    return L
end
#-----------------------------------------------------------------------
"""
fitplane - Solves coefficients of plane fitted to 3 or more points
```
Usage:   B = fitplane(XYZ)

Where:   XYZ - 3xNpts array of xyz coordinates to fit plane to.
               If Npts is greater than 3 a least squares solution
               is generated.

Returns: B   - 4x1 array of plane coefficients in the form
               B[1]*X + B[2]*Y + B[3]*Z + B[4] = 0
               The magnitude of B is 1.
```
See also: ransacfitplane()
"""
# June 2003

function fitplane{T<:Real}(XYZ::Array{T,2})

    (rows, npts) = size(XYZ)

    if rows != 3
        error("Data is not 3D")
    end

    if npts < 3
        error("Too few points to fit plane")
    end

    # Set up constraint equations of the form  AB = 0,
    # where B is a column vector of the plane coefficients
    # in the form   B[1]*X + B[2]*Y + B[3]*Z + B[4] = 0
    A = [XYZ' ones(npts,1)]  # Build constraint matrix

    if npts == 3             # Pad A with zeros
        A = [A; zeros(1,4)]
    end

    (u, d, v) = svd(A)       # Singular value decomposition.

    return B = v[:,4]        # Solution is last column of v.
end


text:
__precompile__()

module PostSPH2

using Printf #To construct string message easily
using StaticArrays

include("..SaveVTK2.jl")

export readVtkArray,
    Cat,
    readVtkParticles,
    ForceVtk,
    FloatingVtkTranslation,
    MassVtk,
    readVtkNames,
    readVtkVariables,
    readBi4Array,
    readBi4Body,
    MkArray

##Hardcoded enum - Cat is "category"
@enum Cat begin
    Points
    Idp
    Vel
    Rhop
    Mass
    Press
    Vol
    Ace
    Vor
    Typ #Typ, since "Type" is illegal to assign
    Mk
end

##Read from binary files directly

## Make two dicts,   one for the enum -> string and one for the data type - Bi4 specific
const varNames = (:key, :offset)

const IdpSearch = [0x41; 0x52; 0x52; 0x41; 0x59; 0x03; 0x00; 0x00; 0x00; 0x49; 0x64; 0x70]
const IdpOffset = 11
const IdpKey = NamedTuple{varNames}([IdpSearch, IdpOffset])

const PosSearch = [0x41; 0x52; 0x52; 0x41; 0x59; 0x03; 0x00; 0x00; 0x00; 0x50; 0x6f; 0x73]
const PosOffset = 11
const PosKey = NamedTuple{varNames}([PosSearch, PosOffset])

const VelSearch = [0x41; 0x52; 0x52; 0x41; 0x59; 0x03; 0x00; 0x00; 0x00; 0x56; 0x65; 0x6c]
const VelOffset = 11
const VelKey = NamedTuple{varNames}([VelSearch, VelOffset])

const RhopSearch =
    [0x41; 0x52; 0x52; 0x41; 0x59; 0x04; 0x00; 0x00; 0x00; 0x52; 0x68; 0x6f; 0x70]
const RhopOffset = 12
const RhopKey = NamedTuple{varNames}([RhopSearch, RhopOffset])

const searchKeyBi4 =
    Dict{Cat,NamedTuple}(Idp => IdpKey, Points => PosKey, Vel => VelKey, Rhop => RhopKey)
const catTypeBi4 =
    Dict{Cat,DataType}(Idp => Int32, Points => Float32, Vel => Float32, Rhop => Float32)
const catArrayBi4 = Dict{Cat,Int64}(Idp => 1, Points => 2, Vel => 2, Rhop => 1)
const catColBi4 = Dict{Cat,Int64}(Idp => 1, Points => 3, Vel => 3, Rhop => 1)

##Lists files in directory and only returns applicable files, ie. "Part_XXXX.bi4"
function _dirFiles(first_file::Bool = false)
    if first_file == false
        files = readdir()
    else
        files = readdir()[1]
    end
    #Operation on dirFiles instantly
    filter!(x -> occursin(r"Part_\d{4}.bi4", x), files)
    return files
end

##transferData is smart since we change the value of an array in place and
# skip a lot of unnecessary allocation steps. It has been tuned to accept
# matrices and vectors, where the type is automatically deduced using "eltype".
_readel(ft::IOStream, typ) = read(ft, typ)
@inline _readel(ft::IOStream, typ::Type{<:SVector{3,T}}) where {T} =
    SVector(_readel(ft, T), _readel(ft, T), _readel(ft, T))

function _transferDataBi4(ft::IOStream, arrayVal::AbstractVector)
    typ = eltype(arrayVal)
    @inbounds for i in eachindex(arrayVal)
        arrayVal[i] = _readel(ft, typ)
    end
end

##If nCol is not 1, then type is static array
_typeMaker(typ, nCol) = nCol == 1 ? typ : SVector{nCol,typ}

#Command for plotting
#println(Char.(key))
#scatter(getindex.(a[1], 1), getindex.(a[1], 3))
#readBi4Array(typ::Cat,Bi4Files::Array{String,1}=_dirFiles()) = readBi4Array(typ, false, Bi4Files)
readBi4Array(typ::Cat, Bi4Files::Array{String,1}) = readBi4Array(typ, false, Bi4Files)
readBi4Array(typ::Cat, SeekNull::Bool, Bi4Files::String) =
    readBi4Array(typ, false, [Bi4Files])
readBi4Array(typ::Cat, Bi4Files::String) = readBi4Array(typ, false, [Bi4Files])
function readBi4Array(
    typ::Cat,
    SeekNull::Bool = false,
    Bi4Files::Array{String,1} = _dirFiles(),
)

    #if true start from top, else do not, use startPos from _Bi4Pos
    nBi4 = size(Bi4Files)[1]

    T = _typeMaker(catTypeBi4[typ], catColBi4[typ])
    j = fill(Array{T,1}(), nBi4) #Less allocs than Vector{Array{T}}(undef,nBi4)

    key = searchKeyBi4[typ].key
    offset = searchKeyBi4[typ].offset

    Threads.@threads for i = 1:nBi4
        ft = open(Bi4Files[i], read = true)
        rf = read(ft)
        startPos = Base._searchindex(rf, key, 1) + offset #1 byte offset
        seek(ft, startPos)

        read(ft, Int64)
        n = read(ft, Int32)
        read(ft, Int32)

        j[i] = zeros(T, n)
        _transferDataBi4(ft, j[i])

        close(ft)
    end
    return j
end

##StaticArrays is hard to use here since it is needed to offset with "Int32", between
# all searches unlike "readBi4Array"
function readBi4Particles(Bi4Files::Array{String,1} = _dirFiles())

    nBi4 = size(Bi4Files)[1]

    j = Vector{Array{Int32,1}}(undef, nBi4)

    ParticleString = ["CaseNp", "CaseNfixed", "CaseNmoving", "CaseNfloat", "CaseNfluid"]
    Threads.@threads for i = 1:nBi4
        ft = open(Bi4Files[i], read = true)
        j_inner = zeros(Int32, 5)
        for k in eachindex(j_inner)
            readuntil(ft, ParticleString[k])
            read(ft, Int32)
            j_inner[k] = read(ft, Int32)
        end
        j[i] = j_inner
        close(ft)
    end
    return j
end

## Function to read the time at current simulation step, as in "XXXX.out"
function readBi4Time(Bi4Files::Array{String,1} = _dirFiles())
    nBi4 = size(Bi4Files)[1]

    T = Float64
    j = Vector{T}(undef, nBi4)

    ParticleString = ["TimeStep"]
    Threads.@threads for i = 1:nBi4
        ft = open(Bi4Files[i], read = true)
        readuntil(ft, ParticleString[1])
        read(ft, Int32)
        j[i] = read(ft, T)
        close(ft)
    end
    return j
end
#Function to only find specific Idps
#for MovingSquare example "Bodies[2][1]""
#In the for loop the first index is the index of the relevant "typ" array,
#the second index is the sorting of this "typ" array corresponding to the "idp"
#array and "start:move" are the number of particles defined by the "Body", where
#0 and 1 indexing from C++ and Julia differences has been taken into account.
readBi4Body(Body, typ, Bi4Files::String) = readBi4Body(Body, typ, [Bi4Files])
function readBi4Body(Body, typ, Bi4Files::Array{String,1} = _dirFiles())
    start = getfield(Body, :beg) + 1    #First idp
    move = start + getfield(Body, :count) - 1  #Number of particles from first idp
    idp_vec = readBi4Array(Idp, false, Bi4Files)
    val_vec = readBi4Array(typ, false, Bi4Files)

    j = similar(val_vec)

    for i = 1:length(j)
        j[i] = val_vec[i][sortperm(idp_vec[i])][start:move]
    end

    return j
end

#Functions

end #PostSPH2


text:
## juliaserver.jl
#
## sets up server to listen on channel 8081
#
## in future it might be preferable to migrate to HTTP.jl
## which has more features than HttpServer.jl and is
## under active development.
## (Tom Price, Dec 2018)

## usage:
## cd("QpcrAnalysis")
## julia -e 'push!(LOAD_PATH,pwd()); include("../juliaserver.jl")'

import HTTP: serve, Request, Response, HandlerFunction, mkheaders, URI, URIs.splitpath
import JSON.json
import Memento: getlogger, gethandlers, debug, info, error
import QpcrAnalysis

## set up logging
logger = getlogger("QpcrAnalysis")
debug(logger, "******************** Julia Server is started ********************")

## headers
const RESPONSE_HEADERS = HTTP.mkheaders([
    "Server"           => "Julia/$VERSION",
    "Content-Type"     => "text/html; charset=utf-8",
    "Content-Language" => "en",
    "Date"             => Dates.format(now(Dates.UTC), Dates.RFC1123Format)])

global julia_running = false

info(logger, "Starting webserver listener on: http://127.0.0.1:8081")
global TCPREF = Ref{Base.TCPServer}()

## set up REST endpoints to dispatch to service functions
HTTP.listen( tcpref = TCPREF) do req::HTTP.Request
	global julia_running

    info(logger, "Julia webserver has received $(req.method) request to http://127.0.0.1:8081$(req.target)")

	## Server control
	if req.target == "/exit"
		close(TCPREF.x)
		# Exit server for debug puposes.
		exit(0)
		# should never come here
		return HTTP.Response(200, RESPONSE_HEADERS; body=JSON.json(Dict(:info => "Julia is exiting.. cya..")))
	elseif req.target == "/farewell"
		# Exit server with no spawn.
		exit(101)
		return HTTP.Response(200, RESPONSE_HEADERS; body=JSON.json(Dict(:info => "Julia is exiting.. cya..")))
	elseif req.target == "/end"
		# Polite way to exit listener. No agressive thread kill.
		info(logger, "Ending Julia webserver")
		close(TCPREF.x)
		return HTTP.Response(200, RESPONSE_HEADERS; body=JSON.json(Dict(:info => "Julia is exiting.")))
	elseif req.target == "/logerror"
		# setting logging level to errors only
		Memento.setlevel!(logger, "error")
		return HTTP.Response(200, RESPONSE_HEADERS; body=JSON.json(Dict(:info => "Julia logging level set to error.")))
	elseif req.target == "/logwarn"
		# setting logging level to warn
		Memento.setlevel!(logger, "warn")
		return HTTP.Response(200, RESPONSE_HEADERS; body=JSON.json(Dict(:info => "Julia logging level set to warn.")))
	elseif req.target == "/logdebug"
		# setting logging level to debug
		Memento.setlevel!(logger, "debug")
		return HTTP.Response(200, RESPONSE_HEADERS; body=JSON.json(Dict(:info => "Julia logging level set to debug.")))
	elseif req.target == "/isbusy"
		# check julia sever status
		return HTTP.Response(200, RESPONSE_HEADERS; body=JSON.json(Dict(:info => (julia_running?"Juliaserver is busy":"Juliaserver is ready"))))
	end

	if julia_running == true
		# Julia can process only one request at a time
		warn(logger, "Julia is still processing an existing request, reject new request with 503")
		return HTTP.Response(503, RESPONSE_HEADERS; body=JSON.json(Dict(:error => "Julia is busy")))
	end

	julia_running = true

    code =
        if req.method == "GET" ## per HTTP RFC, this is actually a POST request because it contains body data
            const nodes = HTTP.URIs.splitpath(req.target)
            if length(nodes) >= 3
                experiment_id = nodes[2]
                action        = nodes[3]
                request_body  = String(req.body)

                ## calls to http://localhost/experiments/0/
                ## will activate a slow test mode
                kwargs = Dict{Symbol,Bool}(
                    (experiment_id == "0") ? :verify => true : ())

                ## dispatch request to Julia engine
                debug(logger, "Calling QpcrAnalysis.dispatch()")
                success, response_body =
                    QpcrAnalysis.dispatch(action, request_body; kwargs...)
                debug(logger, "Returning from QpcrAnalysis.dispatch(): success=$success")
                code = (success) ? 200 : 500
            else ## length(nodes) < 3
                404
            end
        else ## not GET
            404
        end
    (code == 404) && (response_body = JSON.json(Dict(:error => "not found")))
    debug(logger, "Julia finish processing: status: $code, response body: $response_body")

    julia_running = false

    return HTTP.Response(code, RESPONSE_HEADERS; body=response_body)

end ## HTTP.serve

info(logger, "Server module quitting...")

#

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
using CSV
using StatsBase

include("common.jl")

_shifted_geomean(vals, shift) = StatsBase.geomean(vals .+ shift) - shift

df = CSV.read(RESULT_PATH)

# Remove configurations that might be in data set, but that we
# do not wish to include in our analysis.
df = filter(df) do row
    for (solver, method) in CONFIGS_TO_ANALYZE
        if row.solver == solver && row.formulation == method
            return true
        end
    end
    return false
end

df.fam_instance = df.family .* "-" .* df.instance

# Figure out the best primal objective value available in dataset
# for each instance.
df[!, :best_primal] .= 0.0
for i in 1:size(df, 1)
    row = df[i, :]
    best_primal_value = minimum(df[(df.family .== row.family) .& (df.instance .== row.instance), :true_primal])
    df.best_primal[i] = best_primal_value
end

# Add a column to record MIP gap, with respect to best primal solution in data set.
function _mip_gap(row)
    return 100 * abs(row.dual - row.best_primal) / abs(row.best_primal)
end
df[!, :gap] .= 0.0
for i in 1:size(df, 1)
    df.gap[i] = _mip_gap(df[i,:])
end

# Split off problem classes based on how easy they were. We will do separate analysis
# on each problem class.
PROBLEM_CLASS_DFS = Dict()
PROBLEM_CLASS_DFS["solved"] = filter(df) do row
    family, instance = row.family, row.instance
    slice = df[(df.family .== row.family) .& (df.instance .== row.instance), :]
    return all(slice.optimal .== "true")
end

PROBLEM_CLASS_DFS["contested"] = filter(df) do row
    family, instance = row.family, row.instance
    slice = df[(df.family .== row.family) .& (df.instance .== row.instance), :]
    return any(slice.optimal .== "true") && !all(slice.optimal .== "true")
end

PROBLEM_CLASS_DFS["unsolved"] = filter(df) do row
    family, instance = row.family, row.instance
    slice = df[(df.family .== row.family) .& (df.instance .== row.instance), :]
    return !any(slice.optimal .== "true")
end

for (class, subset_df) in PROBLEM_CLASS_DFS
    total_instance_count = length(unique(subset_df.fam_instance))
    total_instance_count == 0 && continue
    println("Class: $class ($total_instance_count instances)")
    println("="^(30))

    fails = Dict(pair => 0 for pair in CONFIGS_TO_ANALYZE)
    best_bound = Dict(pair => 0 for pair in CONFIGS_TO_ANALYZE)
    for fam_instance in unique(subset_df.fam_instance)
        slice = subset_df[subset_df.fam_instance .== fam_instance, :]
        if size(slice, 1) == 0 error("Unexpected case, aborting.") end
        best_solve_time = minimum(slice.solve_time)
        best_gap = minimum(slice.gap)
        for (solver, formulation) in CONFIGS_TO_ANALYZE
            sub_slice = slice[(slice.solver .== solver) .& (slice.formulation .== formulation), :]
            solve_times = sub_slice.solve_time
            @assert size(solve_times, 1) == 1
            gaps = sub_slice.gap
            solve_time = minimum(solve_times)
            gap = minimum(gaps)
            best_primal = minimum(sub_slice.best_primal)
            if solve_time >= 0.999 * TIME_LIMIT
                fails[(solver, formulation)] += 1
            end
            # Either attain best bound, or reach Gurobi's optimality cutoff,
            # which is 1e-4 * |primal_bound|
            if gap == best_gap || gap < 1e-4 * abs(best_primal)
                best_bound[(solver, formulation)] += 1
            end
        end
    end

    time_shift = minimum(subset_df.solve_time)
    gap_shift = max(1e-4, minimum(subset_df.gap))

    for (solver, formulation) in CONFIGS_TO_ANALYZE
        slice = subset_df[(subset_df.solver .== solver) .& (subset_df.formulation .== formulation), :]
        println("$solver -- $formulation")
        println("-"^(length(solver) + length(formulation) + 4))
        println("  * solve time: ", _shifted_geomean(slice.solve_time, time_shift), " sec")
        println("  * MIP gap:    ", _shifted_geomean(slice.gap, gap_shift))
        num_fails = fails[(solver, formulation)]
        num_best_bound = best_bound[(solver, formulation)]
        println("  * Fails:      ", num_fails, " / ", total_instance_count)
        println("  * Best bound: ", num_best_bound, " / ", total_instance_count)
    end
end

max_stars_repo_name: joehuchette/compact-mips-for-qp

max_issues_repo_name: joehuchette/compact-mips-for-qp

max_forks_repo_name: joehuchette/compact-mips-for-qp


text:
using Documenter

push!(LOAD_PATH, "../src/")
using InterProcessCommunication

DEPLOYDOCS = (get(ENV, "CI", nothing) == "true")

makedocs(
    sitename = "Inter-Process Communication",
    format = Documenter.HTML(
        prettyurls = DEPLOYDOCS,
    ),
    authors = "Éric Thiébaut and contributors",
    pages = ["index.md", "semaphores.md", "sharedmemory.md", "reference.md"]
)

if DEPLOYDOCS
    deploydocs(
        repo = "github.com/emmt/InterProcessCommunication.jl.git",
    )
end


text:
"""
Pre-fab acset schemas+semagrams schemas
"""
module Examples

export TheoryPetri, ReactionNet, ReactionNetSema,
  TheoryDirectedPortGraph, DirectedPortGraph, DirectedPortGraphSema,
  TheoryColoredDPG, ColoredDPG, ColoredDPGSema

using Catlab.Present, Catlab.CSetDataStructures
using ..Schema
using ..Boxes

@present TheoryReactionNet(FreeSchema) begin
  (T,S,I,O)::Ob
  it::Hom(I,T)
  is::Hom(I,S)
  ot::Hom(O,T)
  os::Hom(O,S)
  N::Data
  L::Data
  rate::Attr(T,N)
  concentration::Attr(S,N)
  species_label::Attr(S,L)
  transitions_label::Attr(T,L)
end

const ReactionNet = ACSetType(TheoryReactionNet)

@semagramschema ReactionNetSema(TheoryReactionNet) begin
  @box S Circle :species_label
  @box T Square :transitions_label
  @wire I(is,it)
  @wire O(ot,os)
  @data N Numeric
end

@present TheoryDirectedPortGraph(FreeSchema) begin
  Box::Ob
  IPort::Ob
  OPort::Ob
  Wire::Ob
  ibox::Hom(IPort,Box)
  obox::Hom(OPort,Box)
  src::Hom(Wire,OPort)
  tgt::Hom(Wire,IPort)

  String::Data
  label::Attr(Box,String)
end

const DirectedPortGraph = ACSetType(TheoryDirectedPortGraph)

@semagramschema DirectedPortGraphSema(TheoryDirectedPortGraph) begin
  @box Box Square :label
  @port IPort(ibox) "Input"
  @port OPort(obox) "Output"
  @wire Wire(src,tgt)
  @data String Stringlike
end

@present TheoryColoredDPG <: TheoryDirectedPortGraph begin
  Color::Data
  iportcolor::Attr(IPort,Color)
  oportcolor::Attr(OPort,Color)
  wirecolor::Attr(Wire,Color)
end

const ColoredDPG = ACSetType(TheoryColoredDPG)

@semagramschema ColoredDPGSema(TheoryColoredDPG) begin
  @box Box Square :label
  @port IPort(ibox) "Input"
  @port OPort(obox) "Output"
  @wire Wire(src,tgt)
  @data String Stringlike
  @data Color Stringlike
end

end


text:
function init_model(idealmodel::Type{<:IdealModel},components,userlocations,verbose)
    verbose && @info("""Now creating ideal model:
    $idealmodel""")
    return idealmodel(components;userlocations,verbose)
end

function eos(model::IdealModel, V, T, z=SA[1.0])
    negative_vt(V,T) && return nan_num(V,T,z)
    return N_A*k_B*sum(z)*T * a_ideal(model,V,T,z)
end

function volume(model::IdealModel,p,T,z=SA[1.0];phase=:unknown,threaded=false)
    return sum(z)*R̄*T/p
end


text:
using Quoridor
using Graphs
using Base.Test

@test out_edges(Q.Loc(4,1), g) == [Q.PlayerMove((4, 1),(4, 2))
                                   Q.PlayerMove((4, 1),(3, 1))
                                   Q.PlayerMove((4, 1),(6, 1))]