text:
#========== the macro! ==========#

"""
    @tullio C[i,k] := A[i,j] * B[j,k]
    @tullio F[i,k] := \$α * D[i].field[j] * E[col=k, row=j] + \$β

This is a replacement for `@einsum` which understands a bit more syntax.
The expression on the right is summed over all possible valued of the free index `k`,
and `:=` makes a new array `C`, while `=` and `+=` would write into an existing one.
Scalar arguments should have a dollar sign, like `\$α` or `A[i,\$γ]`.

    @tullio G[i,j] := A[i+x+1, j+y+1] * K[x,y]
    @tullio H[i,j] := A[2i+x, 2j+y]  (x in -1:1, y in -1:1)

Shifts and scaling of indices are allowed, including shifts by other indices.
Ranges can be provided as shown, for under-constrained indices.
If they are over-constrained, shifted indices run over the intersection allowed by all constraints,
while un-shifted indices demand agreement between them (e.g. `axes(A,2) == axes(B,1)` above).

    @tullio (*) L[i] := A[J[k]+2, i] / B[k]^2

This is a product instead of a sum, which could also enabled by writing `L[i] *= ...` (in-place).
You can use any reduction function such as `@tullio (max) M[i,j] := ...`.
Indexing by `J[k]+2` here demands `issubset(J, axes(A,1) .- 2)`.

    @tullio N[j] := sqrt <| M[i,j]^2

Pipe operators `|>` and `<|` apply a function after the sum, here `N ≈ map(norm, eachcol(M))`.
Underscores create functions, e.g. `|> sqrt(_ / V[i])` where clearly `i` must not have been summed.

See the readme for further options.
"""
macro tullio(exs...)
    _tullio(exs...; mod=__module__)
end

function _tullio(exs...; mod=Main)

    opts, ranges, ex = parse_options(exs...)
    if isnothing(ex) # then we simply updated global settings
        return (verbose=_VERBOSE[], fastmath=_FASTMATH[], threads=_THREADS[], grad=_GRAD[], avx=_AVX[], cuda=_CUDA[], tensor=_TENSOR[])
    end

    if opts.tensor && opts.redfun == :+ && isdefined(mod, :TensorOperations) && opts.grad != :Dual
        res = try_tensor(ex, ranges, DotDict(;mod = mod, opts...,
            arrays = Symbol[], indices = [], scalars = Symbol[]))
        if res != nothing # then forward & backward both handled by try_tensor
            return Expr(:block, res...) |> esc
        end
    end

    store = DotDict(; mod = mod, opts...,
        flags = Set{Symbol}(), # set while parsing input
    # Reduction
        redind = Symbol[],
    # Everything writes into leftarray[leftraw...], sometimes with a generated name
        leftraw = [],
        leftind = Symbol[],    # vcat(leftind, redind) is the complete list of loop indices
        leftarray = nothing,
        leftscalar = nothing, # only defined for scalar reduction
        leftnames = Symbol[],  # for NamedDims
    # Whole RHS, without finaliser, plus things extracted:
        right = nothing,
        finaliser = nothing,
        rightind = Symbol[],
        sharedind = Symbol[], # indices appearing on every RHS array, safe for ∇thread
        unsafeleft = Symbol[], # k in A[J[k]] never written to by different threads
        unsaferight = Symbol[], # same for gradient
        arrays = Symbol[],
        scalars = Symbol[],
        cost = 1,
    # Index ranges: first save all known constraints
        constraints = Dict{Symbol,Vector}(), # :k => [:(axes(A,2)), :(axes(B,1))] etc.
        notfree = Symbol[], # indices assigned values i = clamp(j, 1,3) within RHS
        shiftedind = Symbol[],
        pairconstraints = Tuple[], # (:i, :j, entangled range_i, range_j) from A[i+j] etc.
        axisdefs = Expr[],
    # Expressions:
        outpre = Expr[],  # preliminary steps
        outex = Expr[],   # the rest!
    )

    parse_input(ex, store)

    parse_ranges(ranges, store)

    index_ranges(store)

    output_array(store)

    ex = action_functions(store)

    opts.verbose == 2 && verboseprint(store)

    ex |> esc
end

#========== options, etc ==========#

OPTS = Dict(
    :verbose => Any[true, false, 2],
    :fastmath => [true, false],
    :threads => Integer,
    :grad => [false, :Base, :Dual],
    :avx => Integer,
    :cuda => Integer,
    :tensor => [true, false],
    )

_VERBOSE = Ref{Any}(false)
_FASTMATH = Ref(true)
_THREADS = Ref{Any}(true)
_GRAD = Ref{Any}(:Base)
_AVX = Ref{Any}(true)
_CUDA = Ref{Any}(256)
_TENSOR = Ref(true)

function parse_options(exs...)
    opts = Dict{Symbol,Any}(
        :redfun => :+,
        :verbose => _VERBOSE[],
        :fastmath => _FASTMATH[],
        :threads => _THREADS[],
        :grad => _GRAD[],
        :avx => _AVX[],
        :cuda => _CUDA[],
        :tensor => _TENSOR[],
        )
    expr = nothing
    nograd = Symbol[]
    ranges = Tuple[]
    for ex in exs
        # Actual options:
        if ex isa Expr && ex.head == :(=) && haskey(OPTS, ex.args[1])
            checklegal(ex.args[1], ex.args[2])
            opts[ex.args[1]] = ex.args[2]

        # Nograd keyword
        elseif ex isa Expr && ex.head == :(=) && ex.args[1] == :nograd
            if ex.args[2] isa Symbol
                push!(nograd, ex.args[2])
            elseif ex.args[2] isa Expr && ex.args[2].head == :tuple
                append!(nograd, ex.args[2].args)
            else
                error("this accepts nograd=A or nograd=(A,B,C)")
            end

        # Ranges specified outside:
        elseif ex isa Expr && ex.head == :call && ex.args[1] in [:in, :∈]
            push!(ranges, (ex.args[2], ex.args[3]))
        elseif ex isa Expr && ex.head == :tuple && ex.args[1] isa Expr && ex.args[1].args[1] in [:in, :∈]
            for el in ex.args
                el isa Expr && el.head == :call && el.args[1] in [:in, :∈] || error("expected (i ∈ 1:3) but got $el")
                push!(ranges, (el.args[2], el.args[3]))
            end

        # Reduction function
        elseif ex isa Symbol
            opts[:redfun] = ex

        # The main course!
        elseif ex isa Expr
            isnothing(expr) || error("too many expressions! recognised keywords are $(vcat(:nograd, keys(opts)...))")
            expr = ex
        else
            error("not sure what to do with input $ex")
        end
    end
    if isnothing(expr) # if run with no expression, it updates global options
        _VERBOSE[] = opts[:verbose]
        _FASTMATH[] = opts[:fastmath]
        _THREADS[] = opts[:threads]
        _GRAD[] = opts[:grad]
        _AVX[] = opts[:avx]
        _CUDA[] = opts[:cuda]
        _TENSOR[] = opts[:tensor]
    end
    (redfun=opts[:redfun],
        verbose=opts[:verbose],
        fastmath=opts[:fastmath],
        threads=opts[:threads],
        grad=opts[:grad],
        avx=opts[:avx],
        cuda=opts[:cuda],
        tensor=opts[:tensor],
        nograd=nograd,
    ), ranges, expr
end

checklegal(opt, val) =
    if OPTS[opt] isa Vector
        val in OPTS[opt] || error("keyword $opt accepts values [" * join(OPTS[opt], ", ") * "]")
    elseif val isa Expr || val isa Symbol
        # allows threads=64^3 to work
    elseif OPTS[opt] == Integer
        val isa Integer && val >= 0 || error("keyword $opt accepts false or a positive integer")
    end

verboseprint(store) =
    foreach(propertynames(store)) do k
        r = getproperty(store, k) # startswith(string(k), "out") fails?
        k ∉ [:outpre, :outex] && return printstyled("    $k = ", repr(r), "\n", color=:blue)
        printstyled("    $k =\n", color=:blue)
        foreach(ex -> printstyled(verbosetidy(ex) , "\n", color=:green), r)
    end

#========== symbols ==========#

# These only need not to clash with symbols in the input:
RHS, AXIS = :𝓇𝒽𝓈, :𝒶𝓍
ZED, TYP, ACC, KEEP, FINAL = :ℛ, :𝒯, :𝒜𝒸𝒸, :♻, :💀
EPS, DEL, EXPR = :𝜀, :𝛥, :ℰ𝓍
MAKE, ACT! = :ℳ𝒶𝓀ℯ, :𝒜𝒸𝓉!

# @gensym RHS MAKE ACT!
# @gensym AXIS ZED TYP ACC KEEP FINAL
# @gensym EPS DEL EXPR

SYMBOLS = [
    RHS, MAKE, ACT!, AXIS, ZED, TYP, ACC, KEEP, EPS, DEL, EXPR,
    Symbol(:∇, MAKE), Symbol(:∇, ACT!), Symbol(DEL, ZED), Symbol(AXIS, :i),
    ] # to test for leaks

#========== input parsing ==========#

function parse_input(expr, store)

    # Equals sign & friends:
    if @capture_(expr, left_ := right_ )
        push!(store.flags, :newarray)
    elseif @capture_(expr, left_ = right_ )
    elseif @capture_(expr, left_ += right_ )
        push!(store.flags, :plusequals)
        store.redfun == :+ || error("can't use += with reduction $(store.redfun)")
    elseif @capture_(expr, left_ *= right_ )
        push!(store.flags, :plusequals) # slightly abusing the name of the flag!
        if store.redfun == :+ # default, then we change it?
            store.verbose>0 && @info "inferring reduction by *, because of lhs *= rhs"
            store.redfun = :*
        elseif store.redfun == :*
        else
            error("can't use *= with reduction $(store.redfun)")
        end
    else error("can't understand input, expected A[] := B[] (or with =, +=, or *=) got $expr")
    end

    # Left hand side:
    if @capture_(left, Z_[leftraw__] ) || @capture_(left, [leftraw__] )
    elseif left isa Symbol # complete reduction, by writing into a new 1-el array
        push!(store.flags, :newarray, :scalar)
        store.leftscalar = left # because store.leftarray will be the array
        leftraw = [1,] # make a 1D array, not zero
        expr.head == :(+=) && push!(store.scalars, left)
    else
        error("can't understand LHS, expected A[i,j,k], got $left")
    end
    leftraw2 = tidyleftraw(leftraw, store)
    store.leftind = filter(i -> i isa Symbol, leftraw2) # this gives correct outer loop order

    isnothing(Z) && !(:newarray in store.flags) && error("can't write into an array whose name isn't given!")
    Zed = isnothing(Z) ? ZED : Z
    store.leftarray = Zed

    store.leftraw = finishleftraw(leftraw2, store)
    if :newarray in store.flags
        !allunique(store.leftind) && push!(store.flags, :zero) # making diagonals, etc.
        Zed in store.arrays && error("can't create a new array $Zed when this also appears on the right")
    else
        saveconstraints(Zed, leftraw, store, false) # this adds to leftind, e.g. A[2i+1] = ..., is that bad??
        detectunsafe(left, store.unsafeleft, store)
    end

    # Right hand side
    detectunsafe(right, store.unsaferight, store)
    right2 = MacroTools_postwalk(rightwalk(store), right)

    if right2 isa Expr && right2.head == :call && right2.args[1] in (:|>, :<|)
        if right2.args[1] == :|>
            store.finaliser = makefinaliser(right2.args[3], store)
            store.right = MacroTools_postwalk(dollarwalk(store), right2.args[2])
        elseif right.args[1] == :<|
            store.finaliser = makefinaliser(right2.args[2], store)
            store.right = MacroTools_postwalk(dollarwalk(store), right2.args[3])
        end
        if :scalar in store.flags
            # scalar threaded reduction won't work with nontrivial finalisers
            store.threads = false
        end
    else
        store.right = MacroTools_postwalk(dollarwalk(store), right2)
        store.finaliser = :identity
    end

    unique!(store.scalars)
    unique!(store.arrays)
    unique!(store.leftind)
    store.sharedind = unique!(setdiff(store.sharedind, store.notfree))
    store.rightind = unique!(setdiff(store.rightind, store.notfree))
    unique!(store.outpre) # kill mutiple assertions, and evaluate any f(A) only once

end

rightwalk(store) = ex -> begin
        @nospecialize ex
        # First, this will detect any assignment before it is used:
        if ex isa Expr && ex.head == :(=)
            if ex.args[1] isa Symbol
                push!(store.notfree, ex.args[1])
            elseif ex.args[1] isa Expr && ex.args[1].head == :tuple
                for i in ex.args[1].args
                    i isa Symbol && push!(store.notfree, i)
                end
            end
        end
        ex isa Expr && ex.head == :return && error("can't use return inside body")

        # Second, alter indexing expr. to pull out functions of arrays:
        @capture_(ex, A_[inds__]) || return ex

        if isnothing(arrayonly(A))
            Anew = Symbol(string("≪", A, "≫"))
            push!(store.outpre, :(local $Anew = $A))
            A = Anew
        end

        # Third, save letter A, and what axes(A) says about indices:
        push!(store.arrays, arrayonly(A))
        inds = primeindices(inds)
        saveconstraints(A, inds, store, true)

        # Re-assemble RHS with new A, and primes on indices taken care of.
        return :( $A[$(inds...)] )
    end # A1[i][k] should be seen later, with corrected A

arrayonly(A::Symbol) = A   # this is for RHS(i,j,k, A,B,C)
arrayonly(A::Expr) =
    if @capture_(A, B_[inds__]) || @capture_(A, B_.field_)
        return arrayonly(B)
    end # returns nothing from :(f(A)), signal to pull function out.

saveconstraints(A, inds, store, right=true) = begin
    A1 = arrayfirst(A)
    is = Symbol[]
    foreach(enumerate(inds)) do (d,ex)
        is_const(ex) && return
        containsany(ex, store.notfree) && return
        axis_i = length(inds)==1 ? :(eachindex($A1)) : :(axes($A1,$d))
        range_i, i = range_expr_walk(axis_i, ex)
        if i isa Symbol
            push!(is, i)
            ex isa Symbol || push!(store.shiftedind, i)
            v = get!(store.constraints, i, [])
            push!(v, dollarstrip(range_i))
        elseif i isa Tuple # from things like A[i+j]
            push!(is, filter(!isnothing, collect(i))...) # collect for Julia ⩽ 1.3
            push!(store.shiftedind, filter(!isnothing, collect(i))...)
            push!(store.pairconstraints, (i..., dollarstrip.(range_i)...))
        elseif isnothing(i) # from A[J[k]], but A[J[k]+i] goes via store.pairconstraints
            str = "extrema of index $ex must fit within $A1"
            push!(store.outpre, :(issubset($range_i, $axis_i) || error($str)))
        end
    end
    if right
        append!(store.rightind, is)
        if A1 in store.nograd # then don't care whether it sharesindices
        elseif isassigned(store.sharedind)
            shared = intersect(is, store.sharedind)
            empty!(store.sharedind)
            append!(store.sharedind, shared)
        else
            append!(store.sharedind, is)
        end
    else
        append!(store.leftind, is) # why can's this be the only path for store.leftind??
    end
    n = length(inds)
    if n==1
        str = "expected a 1-array $A1, or a tuple"
        push!(store.outpre, :( $A1 isa Tuple || ndims($A1) == 1 || error($str) ))
    else
        str = "expected a $n-array $A1" # already arrayfirst(A)
        push!(store.outpre, :( ndims($A1) == $n || error($str) ))
    end
end

arrayfirst(A::Symbol) = A  # this is for axes(A,d), axes(first(B),d), etc.
arrayfirst(A::Expr) =
    if (@capture_(A, Binds_.field_) && @capture_(Binds, B_[inds__]))
        return :( first($B).$field )
    elseif @capture_(A, B_[inds__])
        return :( first($B) )
    elseif @capture_(A, B_.field_)
        return A
    end

containsany(ex, list) = begin
    out = false
    MacroTools_postwalk(ex) do x
        if x in list
            out = true
        end
        x
    end
    out
end

primeindices(inds) = map(inds) do ex
    ex isa Expr && ex.head == Symbol("'") &&
        return Symbol(ex.args[1], "′") # normalise i''
    ex
end

dollarwalk(store) = ex -> begin
        @nospecialize ex
        ex isa Expr || return ex
        if ex.head == :call
            callcost(ex.args[1], store) # cost model for threading
        elseif ex.head == :$ # interpolation of $c things:
            ex.args[1] isa Symbol || error("you can only interpolate single symbols, not $ex")
            push!(store.scalars, ex.args[1])
            return ex.args[1]
        end
        ex
    end

dollarstrip(expr) = MacroTools_postwalk(expr) do @nospecialize ex
        ex isa Expr && ex.head == :$ && return ex.args[1]
        ex
    end

tidyleftraw(leftraw, store) = begin
    step1 = map(leftraw) do i
        if i isa Expr && i.head == :kw && :newarray in store.flags # then NamedDims wrapper is put on later
                push!(store.leftnames, i.args[1])
                return i.args[2]
        elseif i === :_ # underscores on left denote trivial dimensions
            return 1
        end
        i
    end
    primeindices(step1) # normalise i' to i′
end

finishleftraw(leftraw, store) = map(enumerate(leftraw)) do (d,i)
    if i isa Expr && i.head == :$
        :newarray in store.flags && error("can't fix indices on LHS when making a new array")
        i.args[1] isa Symbol || error("you can only interpolate single symbols, not $ex")
        push!(store.scalars, i.args[1])
        return i.args[1]

    elseif @capture_(i, J_[inds__]) # scatter operation, A[i,J[j,k]] := ...
        push!(store.nograd, J)
        rightwalk(store)(i) # array J viewed as part of RHS, and provides a range for j,k
        inds2 = filter(j->j isa Symbol, tidyleftraw(inds, store))
        append!(store.leftind, inds2) # but j,k aren't to be summed

        ex = :($J[$(tidyleftraw(inds, store)...)])
        if :newarray in store.flags
            ax_i = Symbol(AXIS, string("≪", ex, "≫")) # fake index name, to which to attach a size?
            push!(store.axisdefs, :(local $ax_i = $extremerange($J)))
            push!(store.flags, :zero)
        elseif !(:plusequals in store.flags) # A[i,J[j,k]] += ... doesn't zero
            push!(store.flags, :zero)
        end

        return ex # has primes dealt with
    end
    i
end

detectunsafe(expr, list, store) = MacroTools_postwalk(expr) do ex
        @capture_(ex, A_[inds__]) || return ex
        for i in inds
            MacroTools_postwalk(i) do x
                @capture_(x, B_[inner__]) || return x
                # Now we have found an array which indexes another one, mark its indices unsafe
                append!(list, filter(j -> j isa Symbol, inner))
                unique!(list)
                # and don't compute a gradient for the inner array
                B isa Symbol && push!(store.nograd, B)
                x
            end
        end
        ex
    end

makefinaliser(s::Symbol, store) = s
makefinaliser(expr::Expr, store) = begin
    underscore = false
    out = MacroTools_postwalk(expr) do ex
        if ex == :_
            underscore = true
            return RHS
        elseif @capture_(ex, A_[inds__])
            for i in inds
                i isa Symbol || continue
                i in store.leftind || error("index $i can't be used in finaliser")
            end
        end
        ex
    end
    if underscore
        return dollarstrip(:($RHS -> $out))
    else
        return dollarstrip(ex)
    end
end

function parse_ranges(ranges, store) # now runs after parse_input
    for (i,r) in ranges
        if i isa Expr && i.head == Symbol("'") # catch primes!
            i = Symbol(i.args[1], "′")
        end
        push!(store.rightind, i)
        v = get!(store.constraints, i, [])
        if r isa Expr && r.head == :call && r.args[1] == :(:) && length(r.args) == 3
            # for a literal range, write OneTo(10) or 0:9 directly into constraints
            if r.args[2] == 1 && r.args[3] isa Integer
                push!(v, :(Base.OneTo($(r.args[3]))))
                continue
            elseif r.args[2] isa Integer && r.args[3] isa Integer
                push!(v, r)
                continue
            end
        end
        # for axes(A,2) where A is already available, just save it
        if r isa Expr && r.head == :call && r.args[1] in (:axes, :eachindex) && r.args[2] in store.arrays
            push!(v, r)
            continue
        end
        # for anything else, treat it as a scalar argument
        if r isa Symbol
            push!(store.scalars, r)
            push!(v, r)
        else
            s = Symbol(string("≪", r, "≫"))
            push!(store.outpre, :($s = $r))
            push!(store.scalars, s)
            push!(v, s)
        end
    end
    unique!(store.rightind)
    unique!(store.scalars)
    store.redind = setdiff(store.rightind, store.leftind)
end

#========== index ranges ==========#

function index_ranges(store)

    todo = Set(vcat(store.leftind, store.redind))
    done = Dict{Symbol,Any}()

    for (i,j,r_i,r_j) in store.pairconstraints

        if isnothing(i) # case of A[j + I[k]]
            v = get!(store.constraints, j, [])
            push!(v, r_j)
        elseif isnothing(j)
            v = get!(store.constraints, i, [])
            push!(v, r_i)

        elseif haskey(store.constraints, i) && i in todo
            resolveintersect(i, store, done) # use existing knowledge to fix i's range
            pop!(todo, i)
            v = get!(store.constraints, j, []) # and then allow j's range to depend on that
            push!(v, r_j)
        elseif haskey(store.constraints, j) && j in todo
            resolveintersect(j, store, done)
            pop!(todo, j)
            v = get!(store.constraints, i, [])
            push!(v, r_i)
        end
    end

    for i in todo
        haskey(store.constraints, i) || error("unable to infer range of index $i")
        if i in store.shiftedind
            resolveintersect(i, store, done)
        else
            resolvestrict(i, store, done)
        end
    end

    append!(store.outex, store.axisdefs)

    if store.verbose > 0
        lex = map(i -> Expr(:(=), i, done[i]), store.leftind)
        push!(store.outex, :(@info "left index ranges" $(lex...)))
        if !isempty(store.redind)
            rex = map(i -> Expr(:(=), i, done[i]), store.redind)
            push!(store.outex, :(@info "reduction index ranges" $(rex...)))
        end
    end
end

resolvestrict(i, store, done) = begin
    res = first(store.constraints[i])
    ax_i = Symbol(AXIS, i)
    push!(store.axisdefs, :( local $ax_i = $res ))
    done[i] = res
    for alt in store.constraints[i][2:end] # in which case it shouldn't be a Set
        str = "range of index $i must agree"
        push!(store.axisdefs, :( $alt == $res || error($str) ))
    end
end

resolveintersect(i, store, done) = begin
    res = length(store.constraints[i])==1 ?
        first(store.constraints[i]) : # because intersect(1:3) isa Vector, wtf?
        :( intersect($(store.constraints[i]...)) )
    ax_i = Symbol(AXIS, i)
    push!(store.axisdefs, :( local $ax_i = $res ))
    done[i] = res
end

#========== output array + eltype ==========#

function output_array(store)
    if :newarray in store.flags

        push!(store.outex, :( local $RHS($(store.arrays...), $(store.rightind...)) = $(store.finaliser)($(store.right)) ))

        # Try inference first, usually fine, and avoids scalar evaluation on GPU
        allfirst = map(i -> :(first($(Symbol(AXIS, i)))), store.rightind)
        T0 = Symbol(TYP,0)
        warn = store.verbose>0 ? :(@warn "unable to infer eltype from RHS") : nothing
        push!(store.outex, quote
            local $T0 = Core.Compiler.return_type($RHS, typeof(($(store.arrays...), $(allfirst...))))
            local $TYP = if Base.isconcretetype($T0)
                $T0
            else
                $warn
                typeof($RHS($(store.arrays...), $(allfirst...)))
            end
        end)

        # This now checks for OffsetArrays, and allows A[i,1] := ...
        outaxes = map(store.leftraw) do i
            i isa Integer && i==1 && return :(Base.OneTo(1))
            i isa Symbol && return Symbol(AXIS, i)
            i isa Expr && @capture_(i, J_[inds__]) && return Symbol(AXIS, string("≪", i, "≫"))
            error("can't use index $i on LHS for a new array")
        end

        if !isdefined(store.mod, :OffsetArrays)
            outaxes = map(store.leftraw, outaxes) do i, ax
                ax == :(Base.OneTo(1)) && return ax
                i in store.shiftedind || @capture_(i, J_[inds__]) || return ax
                push!(store.outex, :( first($ax) == 1 || error("to allow indices not starting at 1, OffsetArrays must be visible in the caller's module")))
                return :(Base.OneTo($ax))
            end
        end

        simex = if isempty(store.arrays)
            # :( zeros($TYP, tuple($(outaxes...))) ) # Array{T} doesn't accept ranges... but zero() doesn't accept things like  @tullio [i,j] := (i,j)  i ∈ 2:3, j ∈ 4:5
            :( similar(1:0, $TYP, tuple($(outaxes...))) )
        else
            # parent() is a trick to avoid a NamedDims bug
            :( similar(parent($(store.arrays[1])), $TYP, tuple($(outaxes...),)) )
        end
        if isempty(store.leftnames)
            push!(store.outex, :( local $(store.leftarray) = $simex ))
        else
            nex = :(tuple($(QuoteNode.(store.leftnames)...)))
            push!(store.outex, :( local $(store.leftarray) = NamedDims.NamedDimsArray($simex, $nex) ))
        end

        # Deal with scalar += now: write into array, later read it out:
        if :scalar in store.flags && :plusequals in store.flags
            push!(store.outex, :($(store.leftarray)[$(store.leftraw...)] = $(store.leftscalar)))
        end
    end

    if :zero in store.flags
        push!(store.outex, :( $(store.leftarray) .= false )) # zero($TYP) won't work in-place
    end

end

#========== action functions ==========#

function action_functions(store)

    axisleft = map(i -> Symbol(AXIS, i), setdiff(store.leftind, store.unsafeleft))
    axisred = map(i -> Symbol(AXIS, i), union(store.redind, store.unsafeleft))
    axislist = vcat(axisleft, axisred)
    # Order of these is convenient for threader(), which divides axisleft up freely,
    # divides axisred up with re-starts.
    # This is independent of the grouping inner/outer for make_many_actors().

    #===== constructing loops =====#

    init = store.redfun == :* ? :(one($TYP)) :
        store.redfun == :max ? :(typemin($TYP)) :
        store.redfun == :min ? :(typemax($TYP)) :
        store.redfun == :& ? :(true) :
        store.redfun == :| ? :(false) :
        :(zero($TYP))

    # Right now this would allow *= only with reduction * too. Could separate them:
    # acc=0; acc = acc + rhs; Z[i] = ifelse(keep, acc, Z[i] * acc)
    # But then keep=true can't be used for blocking, which wants to continue the same as acc.

    ex_init = :( $ACC = ifelse(isnothing($KEEP), $init, $ZED[$(store.leftraw...)]) )
    # ex_init = :( $ACC = isnothing($KEEP) ? $init : $ZED[$(store.leftraw...)] ) # more allocations with @avx, not sure why

    ex_iter = :( $ACC = $(store.redfun)($ACC, $(store.right) ) )

    ex_write = if store.finaliser == :identity
        :( $ZED[$(store.leftraw...)] = $ACC )
    else
        # :( $ZED[$(store.leftraw...)] = isnothing($FINAL) ? $ACC : $(store.finaliser)($ACC) )
        :( $ZED[$(store.leftraw...)] = ifelse(isnothing($FINAL), $ACC, $(store.finaliser)($ACC)) )
    end

    ex_nored = if :plusequals in store.flags # meaning always keep = true and final = true, since there's no branch, no need to reduce :identity
        :( $ZED[$(store.leftraw...)] =  $(store.finaliser)($(store.redfun)($ZED[$(store.leftraw...)] ,$(store.right))) )
    else
        :( $ZED[$(store.leftraw...)] =  $(store.finaliser)($(store.right)) )
    end

    if isempty(store.redind)
        make_many_actors(ACT!,
            vcat(:($ZED::AbstractArray{$TYP}), store.arrays, store.scalars, axislist),
            nothing, store.leftind, nothing, Symbol[], ex_nored, nothing, store)
    else
        make_many_actors(ACT!,
            vcat(:($ZED::AbstractArray{$TYP}), store.arrays, store.scalars, axislist),
            nothing, store.leftind, ex_init, store.redind, ex_iter, ex_write, store)
    end

    ∇make = if :newarray in store.flags
        # make_many_actors and backward_definitions both push into store.outpre
        backward_definitions(store)
    else
        nothing
    end

    #===== action! =====#

    ST = :($storage_type($(store.leftarray), $(store.arrays...)))
    keep = (:plusequals in store.flags) ? :true : :nothing
    block = store.threads==false ? nothing :
        store.threads==true ? (BLOCK[] ÷ store.cost) :
        store.threads
    push!(store.outex, quote
        $threader($ACT!, $ST, $(store.leftarray),
            tuple($(store.arrays...), $(store.scalars...),),
            tuple($(axisleft...),), tuple($(axisred...),), $(store.redfun), $block, $keep)
        $(store.leftarray)
    end)
    store.verbose>0 && block != nothing && @info "threading threshold (from cost = $(store.cost))" block

    if :newarray in store.flags
        # then slurp up outex to make a function:
        ex = quote
            let $ACT! = $ACT!
                local @inline function $MAKE($(store.arrays...), $(store.scalars...), )
                    $(store.outex...)
                end
                $Eval($MAKE, $∇make)($(store.arrays...), $(store.scalars...), )
            end
        end

        # wrap pre and out in one let block so that ACT! doesn't escape:
        ex = :(let
            $(store.outpre...)
            $ex
        end)

        # and assign the result if necc:
        if store.leftarray != ZED
            push!(store.outex, :($(store.leftarray) = $ex ))
            return :($(store.leftarray) = $ex )
        elseif :scalar in store.flags
             push!(store.outex, :($(store.leftscalar) = $getonly($ex)))
             return :($(store.leftscalar) = $getonly($ex))
        else # case of [i,j] := ... with no name given
            # push!(store.outex, ex)
            return ex
        end

    else
        # in-place, no MAKE function, but still keep ACT! from escaping
        return :(let
            $(store.outpre...)
            $(store.outex...)
        end)
    end
end

"""
    make_many_actors(f!, args, ex1, [:i,], ex3, [:k,], ex5, ex6, store)

This makes several functions of this form,
decorated as necessary with `@inbouds` or `@avx` etc,
and with appropriate `storage_type` as the first argument.
```
f!(::Type, args..., keep=nothing) where {T}
    ex1
    ex2 = (for i in axis_i
        ex3
        ex4 = (for k in axis_k
            ex5
        end)
        ex6
    end)
end
```
"""
function make_many_actors(act!, args, ex1, outer::Vector, ex3, inner::Vector, ex5, ex6, store, note="")

    ex4 = recurseloops(ex5, inner)
    ex2 = recurseloops(:($ex3; $ex4; $ex6), outer)

    if store.fastmath && isempty(store.notfree)
        push!(store.outpre, quote
            local @inline function $act!(::Type, $(args...), $KEEP=nothing, $FINAL=true) where {$TYP}
                @inbounds @fastmath ($ex1; $ex2)
            end
        end)
    elseif isempty(store.notfree)
        push!(store.outpre, quote
            local @inline function $act!(::Type, $(args...), $KEEP=nothing, $FINAL=true) where {$TYP}
                @inbounds ($ex1; $ex2)
            end
        end)
    else
        push!(store.outpre, quote
            local @inline function $act!(::Type, $(args...), $KEEP=nothing, $FINAL=true) where {$TYP}
                ($ex1; $ex2)
            end
        end)
    end

    expre, exloop, expost = if isempty(outer)
        :($ex1; $ex3), ex4, ex6
    else
        ex1, ex2, nothing
    end

    # Disable @avx for scatter, https://github.com/chriselrod/LoopVectorization.jl/issues/145
    safe = if act! == ACT!
        isempty(store.unsafeleft)
    else # working on ∇act!
        isempty(store.unsaferight)
    end

    if safe && store.avx != false && isdefined(store.mod, :LoopVectorization)
        unroll = store.avx == true ? 0 : store.avx # unroll=0 is the default setting
        info1 = store.verbose>0 ? :(@info "running LoopVectorization actor $($note)") : nothing
        try lex = macroexpand(store.mod, quote

                local @inline function $act!(::Type{<:Array{<:Union{Base.HWReal, Bool}}}, $(args...), $KEEP=nothing, $FINAL=true) where {$TYP}
                    $expre
                    $info1
                    LoopVectorization.@avx unroll=$unroll $exloop
                    $expost
                end

            end) # macroexpand quote
            push!(store.outpre, lex)
            store.verbose == 2 && @info "success wtih LoopVectorization, unroll=$unroll $note"
        catch err
            store.verbose > 0 && @warn "LoopVectorization failed $note" err
        end
    end

    axouter = map(i -> Symbol(AXIS, i), outer)

    if store.cuda > 0 && isdefined(store.mod, :KernelAbstractions)

        kernel = gensym(:🇨🇺)
        asserts = map(ax -> :( first($ax)==1 || error("KernelAbstractions can't handle OffsetArrays here")), axouter)
        sizes = map(ax -> :(length($ax)), axouter)
        try
            if isempty(outer)
                store.verbose > 0 && @warn "using KernelAbstractions with no outer indices, this will be slow"
                outer = [Symbol(EPS, 1)] # fake index name, only appears in @index
                sizes = [:(one(Int))]    # iterate over 1:1
            end
            kex1 = macroexpand(store.mod, quote

                KernelAbstractions.@kernel function $kernel($(args...), $KEEP, $FINAL) where {$TYP}
                    ($(outer...),) = @index(Global, NTuple)
                    ($ex1; $ex3; $ex4; $ex6)
                end

            end)
            push!(store.outpre, kex1)
            if isdefined(store.mod, :CuArrays) && isdefined(store.mod, :CuArray) # old-style, CuArrays.jl
                info2 = store.verbose>0 ? :(@info "running KernelAbstractions + CuArrays actor $($note)") : nothing
                kex2 = quote

                    local @inline function $act!(::Type{<:CuArray}, $(args...), $KEEP=nothing, $FINAL=true) where {$TYP}
                        $info2
                        cu_kern! = $kernel(CUDA(), $(store.cuda))
                        $(asserts...)
                        $ACC = cu_kern!($(args...), $KEEP, $FINAL; ndrange=tuple($(sizes...)))
                        KernelAbstractions.wait($ACC)
                    end

                end
                push!(store.outpre, kex2)
            end
            if isdefined(store.mod, :CUDA) && isdefined(store.mod, :CuArray) # new-style, CUDA.jl, with CUDADevice()
                info2bis = store.verbose>0 ? :(@info "running KernelAbstractions + CUDA actor $($note)") : nothing
                kex2bis = quote

                    local @inline function $act!(::Type{<:CuArray}, $(args...), $KEEP=nothing, $FINAL=true) where {$TYP}
                        $info2bis
                        cu_kern! = $kernel(CUDADevice(), $(store.cuda))
                        $(asserts...)
                        $ACC = cu_kern!($(args...), $KEEP, $FINAL; ndrange=tuple($(sizes...)))
                        KernelAbstractions.wait($ACC)
                    end

                end
                push!(store.outpre, kex2bis)
            end
            info3 = store.verbose>0 ? :(@info "running KernelAbstractions CPU actor $($note)") : nothing
            kex3 = quote

                local @inline function $act!(::Type{<:Array}, $(args...), $KEEP=nothing, $FINAL=true) where {$TYP}
                    $info3
                    cpu_kern! = $kernel(CPU(), 4)
                    $(asserts...)
                    $ACC = cpu_kern!($(args...), $KEEP, $FINAL; ndrange=tuple($(sizes...)))
                    KernelAbstractions.wait($ACC)
                end

            end
            if store.threads==false
                # This CPU kernel can't be called by threader, and so threads=false
                # offers a way to control whether it gets used or not. By default, not.
                push!(store.outpre, kex3)
            end
            store.verbose == 2 && @info "success wtih KernelAbstractions $note"
        catch err
            store.verbose > 0 && @warn "KernelAbstractions failed $note" err
        end
    end

    if act! != ACT! && isempty(store.sharedind) && Threads.nthreads()>1
        store.verbose > 0 && @warn "can't parallelise this gradient, no shared indices $note"
    end
end

recurseloops(ex, list::Vector) =
    if isempty(list)
        return ex
    else
        i = first(list)
        r = Symbol(AXIS, i)
        ex = :(for $i in $r; $ex; end)
        return recurseloops(ex, list[2:end])
    end

#===== define gradient hooks =====#

function backward_definitions(store)
    store.grad == false && return nothing # no gradient wanted

    axisshared = map(i -> Symbol(AXIS, i), setdiff(store.sharedind, store.unsaferight)) # safe to multi-thread
    loopind = vcat(store.leftind, store.redind)
    axisnonshared = map(i -> Symbol(AXIS, i), union(setdiff(loopind, store.sharedind), store.unsaferight))

    axislist = vcat(axisshared, axisnonshared) # this defines the order of arguments of ∇act!

    ok = false
    if store.grad == :Dual && store.redfun == :+
        try
            insert_forward_gradient(axislist, store)
            ok = true
            store.verbose == 2 && @info "using ForwardDiff gradient"
        catch err
            store.verbose > 0 && @warn "ForwardDiff gradient failed" err
        end
    elseif store.grad == :Base
        try
            insert_symbolic_gradient(axislist, store)
            ok = true
            store.verbose == 2 && @info "success wtih Symbolic gradient"
        catch err
            store.verbose > 0 && @warn "Symbolic gradient failed" err
        end
    end

    ok == false && return nothing # failed to make a gradient

    dZ = Symbol(DEL, ZED)
    ∇make = Symbol(:∇, MAKE)
    ∇act! = Symbol(:∇, ACT!)

    gradarrays = map(A -> Symbol(DEL, A), store.arrays)
    # gradscalars = map(A -> Symbol(DEL, A), store.scalars)
    defineempties = map(store.arrays, gradarrays) do A, dA
        if A in store.nograd
            :(local $dA = nothing)
        else
            :( local $dA = fill!(similar($A, Base.promote_type(eltype($A), $TYP)), 0) )
        end
    end
    # append!(defineempties, map((x,dx) -> :($dx = zero(Base.promote_type(typeof($x), $TYP))), store.scalars, gradscalars))
    returns = vcat(gradarrays, map(_->:nothing, store.scalars)) # ?? needs a test!
    # returns = vcat(gradarrays, gradscalars)

    ST = :($storage_type($(gradarrays...), $(store.arrays...)))
    block = store.threads==false ? nothing :
        store.threads==true ? (BLOCK[] ÷ store.cost) :
        store.threads
    push!(store.outpre, quote

        local $∇make = let $∇act! = $∇act!
            local function $∇make($dZ::AbstractArray{$TYP}, $ZED, $(store.arrays...), $(store.scalars...), ) where {$TYP}
                $(defineempties...)
                $(store.axisdefs...)
                $∇threader($∇act!, $ST,
                    tuple($(gradarrays...), $dZ, $ZED, $(store.arrays...), $(store.scalars...),),
                    tuple($(axisshared...),), tuple($(axisnonshared...), ), $block)
                return ($(returns...),)
            end
        end

    end)

    return ∇make
end

fillarrayreplace(rhs, dZ) = MacroTools_postwalk(rhs) do @nospecialize ex
        @capture_(ex, A_[inds__]) && A==dZ || return ex
        return Symbol(dZ, :_value)
    end

#========== the end ==========#


text:
"""
    SEBox{N}(axes; [r])

The N-dimensional structuring element with all elements connected. This is a special case of
[`SEMask`](@ref) that ImageMorphology algorithms might provide
optimized implementation.

It is recommended to use [`strel_box`](@ref) and [`strel_type`](@ref):

```jldoctest; setup=:(using ImageMorphology; using ImageMorphology.StructuringElements)
julia> se = strel_box((3, 3)) # C8 connectivity
3×3 SEBoxArray{2, UnitRange{$Int}} with indices -1:1×-1:1:
 1  1  1
 1  1  1
 1  1  1

julia> strel_type(se)
SEBox{2, UnitRange{$Int}}((-1:1, -1:1), (1, 1))

julia> se = centered(collect(se)) # converted to normal centered array
3×3 OffsetArray(::Matrix{Bool}, -1:1, -1:1) with eltype Bool with indices -1:1×-1:1:
 1  1  1
 1  1  1
 1  1  1

julia> strel_type(se)
SEMask{2}()
```
"""
struct SEBox{N,R} <: MorphologySE{N}
    axes::NTuple{N,R}
    r::Dims{N}
    function SEBox{N,R}(axes::NTuple{N,R}, r::Dims{N}) where {N,R<:AbstractUnitRange{Int}}
        if !all(r -> first(r) == -last(r), axes)
            throw(ArgumentError("axes must be symmetric along each dimension"))
        end
        return new{N,R}(axes, r)
    end
end
function SEBox{N}(ax::NTuple{N,R}; r=map(R -> length(R) ÷ 2, ax)) where {N,R}
    r = r isa Integer ? ntuple(_ -> r, N) : r
    return SEBox{N,R}(ax, r)
end

"""
    SEBoxArray(se::SEBox)

The instantiated array object of [`SEBox`](@ref SEBox).
"""
struct SEBoxArray{N,R<:AbstractUnitRange{Int}} <: MorphologySEArray{N}
    axes::NTuple{N,R}
    r::Dims{N}
end
SEBoxArray(se::SEBox{N,R}) where {N,R} = SEBoxArray{N,R}(se.axes, se.r)
strel_type(A::SEBoxArray{N}) where {N} = SEBox{N}(A.axes; r=A.r)
strel_size(se::SEBoxArray) = @. 1 + 2 * se.r

@inline Base.axes(A::SEBoxArray) = A.axes
@inline Base.size(A::SEBoxArray) = map(length, axes(A))
@inline function Base.getindex(A::SEBoxArray, inds::Int...)
    @inbounds for i in 1:length(inds)
        if abs(inds[i]) > A.r[i]
            return false
        end
    end
    return true
end

"""
    strel_box(A; r=1)
    strel_box(size; r=size .÷ 2)

Construct the N-dimensional structuring element (SE) with all elements in the local window
connected.

If image `A` is provided, then the SE size will be `(2r+1, 2r+1, ...)` with default
half-size `r=1`. If `size` is provided, the default `r` will be `size .÷ 2`. The default
`dims` will be all dimensions, that is, `(1, 2, ..., length(size))`.

```jldoctest; setup=:(using ImageMorphology; using ImageMorphology.StructuringElements)
julia> img = rand(64, 64);

julia> strel_box(img)
3×3 SEBoxArray{2, UnitRange{$Int}} with indices -1:1×-1:1:
 1  1  1
 1  1  1
 1  1  1

julia> strel_box(img; r=2)
5×5 SEBoxArray{2, UnitRange{$Int}} with indices -2:2×-2:2:
 1  1  1  1  1
 1  1  1  1  1
 1  1  1  1  1
 1  1  1  1  1
 1  1  1  1  1

julia> strel_box((5,5); r=(1,2))
5×5 SEBoxArray{2, UnitRange{$Int}} with indices -2:2×-2:2:
 0  0  0  0  0
 1  1  1  1  1
 1  1  1  1  1
 1  1  1  1  1
 0  0  0  0  0
```

!!! note "specialization and performance"
    The box shape `SEBox` is a special type for which many morphology algorithms may provide
    efficient implementations. For this reason, if one tries to collect an `SEBoxArray` into
    other array types (e.g. `Array{Bool}` via `collect`), then a significant performance
    drop is very likely to occur.

See also [`strel`](@ref) and [`strel_box`](@ref).
"""
function strel_box(A::AbstractArray{T,N}, dims=coords_spatial(A); r::Union{Nothing,Dims{N},Int}=nothing) where {T,N}
    dims = _to_dims(Val(N), dims)
    sz, r = if isnothing(r)
        ntuple(i -> !isempty(dims) && in(i, dims) ? 3 : 1, N), 1
    elseif r isa Dims{N}
        ntuple(i -> !isempty(dims) && in(i, dims) ? 2r[i] + 1 : 1, N), r
    elseif r isa Integer
        ntuple(i -> !isempty(dims) && in(i, dims) ? 2r + 1 : 1, N), r
    end
    return strel_box(sz, dims)
end
function strel_box(sz::Dims{N}, dims=ntuple(identity, N); r::Union{Nothing,Dims{N},Int}=nothing) where {N}
    dims = _to_dims(Val(N), dims)
    all(isodd, sz) || throw(ArgumentError("size should be odd integers"))
    radius = if isnothing(r)
        ntuple(i -> !isempty(dims) && in(i, dims) ? sz[i] ÷ 2 : 0, N)
    elseif r isa Dims{N}
        r
    elseif r isa Integer
        ntuple(i -> !isempty(dims) && in(i, dims) ? r : 0, N)
    end
    ax = map(r -> (-r):r, sz .÷ 2)
    return SEBoxArray(SEBox{N}(ax; r=radius))
end

ext: jl


text:
export SolverProjGradDesc

struct SolverProjGradDesc{T<:Real} <:AbstractNextBasisSolver
    maxiter::Int
    batchsize::Int
    projiter::Int
    stepsize::T
    ftol::T
end

function SolverProjGradDesc(;
    maxiter::Int=100,
    batchsize::Int=20,
    projiter::Int=10,
    stepsize::T=.01,
    ftol::T=1e-4
    ) where T
    SolverProjGradDesc{T}(maxiter, batchsize, projiter, stepsize, ftol)
end

function (self::SolverProjGradDesc{T})(prob::Problem{L, B, CS},
                            r::Vector{T},
                            uinit::Vector{T}) where {L, B<:AbUniDirBasis, CS, T}
    u = copy(uinit)
    n = sizen(prob)
    obj = r'*calc(B, prob.X, u)
    iter = 0

    for i in 1:self.maxiter
        for j in 1:self.batchsize:n
            j2 = min(n, j+self.batchsize-1)
            X = view(prob.X, j:j2, :)
            J = jacobian(B, X, u)
            u -= self.stepsize * J' * r[j:j2]

            # Project
            if (iter .% self.projiter == 0) || (j2 == n)
                u = projectΩ(prob, u)
            end
            iter += 1
        end
        # Stopping
        Δobj = r'*calc(B, prob.X, u) - obj
        obj += Δobj
        abs(Δobj) < self.ftol && break
        printfmtln("\t\t\t obj={:8.4f}", obj)
    end
    u, obj
end


text:
using Infiltrator

# c function documentation at https://mccdaq.github.io/daqhats/c.html
# Global functions and data - https://mccdaq.github.io/daqhats/c.html

"""
	function printError(resultCode)

Print error code text from error code number
"""
function printError(resultCode)
	# map resultCode to descriptive string
	resultDict = Dict{Int32, String}(
		0 => "Success, no errors",
		-1 => "A parameter passed to the function was incorrect.",
		-2 => "The device is busy.",
		-3 => "There was a timeout accessing a resource.",
		-4 => "There was a timeout while obtaining a resource lock.",
		-5 => "The device at the specified address is not the correct type.",
		-6 => "A needed resource was not available.",
		-7 => "Could not communicate with the device.",
		-10 => "Some other error occurred.")
	if resultCode != 0
		# @show(resultCode)
		error(resultDict[resultCode])
	end
end

"""
	function source_convert(source::Symbol)

Convert the source symbol to an ccall integer after checking for validity
"""
function source_convert(source::Symbol)
	# ensure source is valid
	source_set = Set{Symbol}([:SOURCE_LOCAL, :SOURCE_MASTER, :SOURCE_SLAVE])
	if !issubset(Set{Symbol}([source]), source_set)
		error("source is $source, but must be one of $source_set")
	end
	source_dict = Dict{Symbol, UInt8}(
	:SOURCE_LOCAL  => 0,
	:SOURCE_MASTER => 1,
	:SOURCE_SLAVE  => 2)

	return source_dict[source]
end

"""
	function mode_convert(mode::Symbol)

Convert the mode symbol to a ccall integer after checking for validity
"""
function mode_convert(mode::Symbol)
	# ensure trigger is valid
	mode_set = Set{Symbol}([:TRIG_RISING_EDGE, :TRIG_FALLING_EDGE,
		:TRIG_ACTIVE_HIGH, :TRIG_ACTIVE_LOW])
	if !issubset(Set{Symbol}([mode]), mode_set)
		error("mode must be one of $mode_set")
	end

	mode_dict = Dict{Symbol, UInt8}(
	:TRIG_RISING_EDGE  => 0,
	:TRIG_FALLING_EDGE => 1,
	:TRIG_ACTIVE_HIGH  => 2,
	:TRIG_ACTIVE_LOW   => 3)
	return mode_dict[mode]
end

# mcc172 functions - https://mccdaq.github.io/daqhats/c.html

"""
	mcc172_open(address::Integer)

Open a connection to the MCC 172 device at the specified address.
"""
function mcc172_open(address::Integer)
	resultCode = ccall((:mcc172_open, "libdaqhats.so"),
	Cint, (UInt8,), address)
	printError(resultCode)
	return resultCode
end

"""
	mcc172_close(address::Integer)

Close a connection to an MCC 172 device and free allocated resources.
"""
function mcc172_close(address::Integer)
	resultCode = ccall((:mcc172_close, "libdaqhats.so"),
	Cint, (UInt8,), address)
	printError(resultCode)
	return resultCode
end

"""
	mcc172_is_open(address::Integer)

Check if an MCC 172 is open.

Return `true' if open, 'false' if not open.
"""
function mcc172_is_open(address::Integer)
# 1 if open, 0 if not open
	resultCode = ccall((:mcc172_is_open, "libdaqhats.so"),
	Cint, (UInt8,), address)
	return Bool(resultCode)
end

# create struct at top level for device information
mutable struct MCC172DeviceInfo
    NUM_AI_CHANNELS::UInt8    	# The number of analog input channels. UInt8(2)
    AI_MIN_CODE::Int32			# The minimum uncalibrated ADC code. (-8,388,608 or -2^23)
    AI_MAX_CODE::Int32			# The maximum uncalibrated ADC code. (8,388,607 or 2^23-1)
    AI_VOLTAGE_MIN::Float64		# The input voltage corresponding to the minimum code. (-5V)
    AI_VOLTAGE_MAX::Float64		# The input voltage corresponding to the maximum code. (+5V - 1*LSB)
    AI_RANGE_MIN::Float64		# The minimum voltage of the input range. (-5v)
    AI_RANGE_MAX::Float64		# The maximum voltage of the input range. (+5v)
end

"""
	mcc172_info()

Return device information for all MCC 172s.
"""
function mcc172_info()
#   This line recalculates the data from the ccall.
#	mcc172_device_info = MCC172DeviceInfo(2, -2^23, 2^23-1, -5.0, 5.0-10/2^24, -5.0, 5.0)

	output_ptr = ccall((:mcc172_info, "libdaqhats"),
	Ptr{MCC172DeviceInfo}, ())
#	if output_ptr == C_NULL  # Assumed, need to check
#		throw(OutOfMemoryError())
#	end
	mcc172_device_info = unsafe_load(Ptr{MCC172DeviceInfo}(output_ptr))
	return mcc172_device_info
end

"""
	mcc172_blink_led(address::Integer, count::Integer)

Blink the LED on the MCC 172.

Passing 0 for count will result in the LED blinking continuously until
the board is reset or `mcc172_blink_led()` is called again with a non-zero
value for count.

count: The number of times to blink (0 - 255)
"""
function mcc172_blink_led(address::Integer, count::Integer)
	resultCode = ccall((:mcc172_blink_led, "libdaqhats.so"),
	Cint, (UInt8, UInt8), address, count)
	printError(resultCode)
end

"""
	mcc172_firmware_version(address::Integer)

Return the board firmware version.

The version will be in BCD hexadecimal with the high byte as the major
version and low byte as minor, i.e. 0x0103 is version 1.03.
"""
function mcc172_firmware_version(address::Integer)
	version = Ref{UInt16}()
	resultCode = ccall((:mcc172_firmware_version, "libdaqhats.so"),
	Cint, (UInt8, Ref{Cushort}), address, version)
	printError(resultCode)
	return version[]
end

"""
	mcc172_serial(address::Integer)

Read the MCC 172 serial number.
"""
function mcc172_serial(address::Integer)
	serial = Vector{UInt8}(undef, 9)  # initialize serial
	resultCode = ccall((:mcc172_serial, "libdaqhats.so"),
	Cint, (UInt8, Ptr{Cchar}), address, serial)
	printError(resultCode)
	return unsafe_string(pointer(serial))
end

"""
	mcc172_calibration_date(address::Integer)

Read the MCC 172 calibration date.
Date is a string (format “YYYY-MM-DD”)
"""
function mcc172_calibration_date(address::Integer)
	# calDate = "YYYY-MM-DD"
	calDate = Vector{UInt8}(undef,11) # initialized string
	resultCode = ccall((:mcc172_calibration_date, "libdaqhats.so"),
	Cint, (UInt8, Ptr{Cchar}), address, calDate)
	printError(resultCode)
	# return calDate
	calDate[end] = 0
	return unsafe_string(pointer(calDate))
end

"""
	mcc172_calibration_coefficient_read(address::Integer, channel::Integer)

Read the MCC 172 calibration coefficients for a single channel.

The coefficients are applied in the library as:
`calibrated_ADC_code = (raw_ADC_code - offset) * slope`
"""
function mcc172_calibration_coefficient_read(address::Integer, channel::Integer)
	# calibrated_ADC_code = (raw_ADC_code - offset) * slope
	slope = offset = Ref{Cdouble}()
	resultCode = ccall((:mcc172_calibration_coefficient_read, "libdaqhats.so"),
	Cint, (UInt8, UInt8, Ref{Cdouble}, Ref{Cdouble}), address, channel, slope, offset)
	printError(resultCode)
	return (slope[], offset[])
end

"""
	mcc172_calibration_coefficient_write(address::Integer, channel::Integer, slope::AbstractFloat, offset::AbstractFloat)

Temporarily write the MCC 172 calibration coefficients for a single channel.

The user can apply their own calibration coefficients by writing to these
values. The values will reset to the factory values from the EEPROM
whenever `mcc172_open()` is called. This function will fail and
return RESULT_BUSY if a scan is active when it is called.

The coefficients are applied in the library as:
`calibrated_ADC_code = (raw_ADC_code - offset) * slope`
"""
function mcc172_calibration_coefficient_write(address::Integer, channel::Integer, slope::AbstractFloat, offset::AbstractFloat)
	# calibrated_ADC_code = (raw_ADC_code - offset) * slope
	resultCode = ccall((:mcc172_calibration_coefficient_write, "libdaqhats.so"),
	Cint, (UInt8, UInt8, Cdouble, Cdouble), address, channel, slope, offset)
	printError(resultCode)
	return nothing
end

"""
	mcc172_iepe_config_read(address::Integer, channel::Integer)

Read the MCC 172 IEPE configuration for a single channel.

	address: The board address (0 - 7). Board must already be opened.
	channel: The channel number (0 - 1).
Return:

	config: Receives the configuration for the specified channel:
		0: IEPE power off
		1: IEPE power on

"""
function mcc172_iepe_config_read(address::Integer, channel::Integer)
	# 0: IEPE power off; 1: IEPE power on
	config = Ref{UInt8}()
	resultCode = ccall((:mcc172_iepe_config_read, "libdaqhats.so"),
	Cint, (UInt8, UInt8, Ref{Cuchar}), address, channel, config)
	printError(resultCode)
	return config[]
end

"""
	mcc172_iepe_config_write(address::Integer, channel::Integer, config::Integer)

Write the MCC 172 IEPE configuration for a single channel.

Writes the new IEPE configuration for a channel. This function will fail
and return RESULT_BUSY if a scan is active when it is called.

	address: The board address (0 - 7). Board must already be opened.
	channel: The channel number (0 - 1)
Return:

	config: IEPE configuration for specified channel:
			false: IEPE power off, true: IEPE power on
"""
function mcc172_iepe_config_write(address::Integer, channel::Integer, config::Union{Bool, Integer})
	resultCode = ccall((:mcc172_iepe_config_write, "libdaqhats.so"),
	Cint, (UInt8, UInt8, UInt8), address, channel, config)
	printError(resultCode)
	return nothing
end

"""
	mcc172_a_in_sensitivity_read(address::Integer, channel::Integer)

Read the MCC 172 analog input sensitivity scaling factor for a single channel.

The sensitivity is specified in mV / engineering unit. The default value
when opening the library is 1000, resulting in no scaling of the input voltage.
"""
function mcc172_a_in_sensitivity_read(address::Integer, channel::Integer)
	# The sensitivity is specified in mV / Engineering Unit
	sensitivity = Ref{Float64}()
	resultCode = ccall((:mcc172_a_in_sensitivity_read, "libdaqhats.so"),
	Cint, (UInt8, UInt8, Ref{Cdouble}), address, channel, sensitivity)
	printError(resultCode)
	return sensitivity[]
end

"""
	mcc172_a_in_sensitivity_write(address::Integer, channel::Integer, sensitivity::AbstractFloat)

Write the MCC 172 analog input sensitivity scaling factor for a single channel.

This applies a scaling factor to the analog input data so it returns
values that are meaningful for the connected sensor.

The sensitivity is specified in mV / engineering unit. The default value
when opening the library is 1000, resulting in no scaling of the input
voltage. Changing this value will not change the values reported by
mcc172_info() since it is simply sensor scaling applied to the data
before returning it.

Examples:

A seismic sensor with a sensitivity of 10 V/g. Set the sensitivity to
10,000 and the returned data will be in units of g.

A vibration sensor with a sensitivity of 100 mV/g. Set the sensitivity
to 100 and the returned data will be in units of g.

This function will fail and return RESULT_BUSY if a scan is active when
it is called.
"""
function mcc172_a_in_sensitivity_write(address::Integer, channel::Integer, sensitivity::AbstractFloat)
	# The sensitivity is specified in mV / Engineering Unit
	resultCode = ccall((:mcc172_a_in_sensitivity_write, "libdaqhats.so"),
	Cint, (UInt8, UInt8, Cdouble), address, channel, sensitivity)
	printError(resultCode)
	return nothing
end

"""
	mcc172_a_in_clock_config_read(address::Integer)

Read the sampling clock configuration.

This function will return the sample clock configuration and rate. If the clock is configured for local or master source, then the rate will be the internally adjusted rate set by the user. If the clock is configured for slave source, then the rate will be measured from the master clock after the synchronization period has ended. The synchronization status is also returned.

The clock source will be one of the following values:

	SOURCE_LOCAL = 0: The clock is generated on this MCC 172 and not shared with other MCC 172s.
	SOURCE_MASTER = 1: The clock is generated on this MCC 172 and is shared as the master clock for other MCC 172s.
	SOURCE_SLAVE = 2: No clock is generated on this MCC 172, it receives its clock from the master MCC 172.
The data rate will not be valid in slave mode if synced is equal to 0.
The device will not detect a loss of the master clock when in slave mode;
it only monitors the clock when a sync is initiated.

Return Result code, RESULT_SUCCESS if successful
Parameters:

	address: The board address (0 - 7). Board must already be opened.
	clock_source::UInt8 Receives the ADC clock source, one of the source type values.
	sample_rate_per_channel::Float64 Receives the sample rate in samples per second per channel
	sync::Bool Receives the syncronization status (false: sync in progress, true: sync complete)
"""
function mcc172_a_in_clock_config_read(address::Integer)
	# The clock can be SOURCE_LOCAL, SOURCE_MASTER, OR SOURCE_SLAVE
	clock_source = synced = Ref{UInt8}()
	sample_rate_per_channel = Ref{Float64}(0.0)
	resultCode = ccall((:mcc172_a_in_clock_config_read, "libdaqhats.so"),
	Cint, (UInt8, Ref{Cuchar}, Ref{Cdouble}, Ref{Cuchar}),
	address, clock_source, sample_rate_per_channel, synced)
	printError(resultCode)

	# synced[]: Receives  syncronization status (0x00: sync in progress, 0x01: sync complete)
	if synced[] == 0x00
		sync = false
	elseif synced[] == 0x01
		sync = true
	else
		error("synced[] must be 0x00 or 0x01")
	end
	return (clock_source[], sample_rate_per_channel[], sync)
end

"""
	mcc172_a_in_clock_config_write(address::Integer, clock_source::Symbol, sample_rate_per_channel::Real)

Write the sampling clock configuration.

This function will configure the ADC sampling clock. The default
configuration after opening the device is local mode, 51.2 KHz data rate.

The clock_source must be one of:

	SOURCE_LOCAL (0): The clock is generated on this MCC 172 and not shared with other MCC 172s.
	SOURCE_MASTER (1): The clock is generated on this MCC 172 and is shared as the master clock for other MCC 172s. All other MCC 172s must be configured for local or slave clock.
	SOURCE_SLAVE (2): No clock is generated on this MCC 172, it receives its clock from the master MCC 172.

The ADCs will be synchronized so they sample the inputs at the same time.
This requires 128 clock cycles before the first sample is available.
When using a master - slave clock configuration for multiple MCC 172s
there are additional considerations:

There should be only one master device; otherwise, you will be connecting
multiple outputs together and could damage a device.
Configure the clock on the slave device(s) first, master last. The
synchronization will occur when the master clock is configured, causing
the ADCs on all the devices to be in sync.

If you change the clock configuration on one device after configuring the
 master, then the data will no longer be in sync. The devices cannot
 detect this and will still report that they are synchronized. Always
 write the clock configuration to all devices when modifying the configuration.
Slave devices must have a master clock source or scans will never complete.
A trigger must be used for the data streams from all devices to start
on the same sample.

The MCC 172 can generate an ADC sampling clock equal to 51.2 kHz divided
by an integer between 1 and 256. The `data_rate_per_channel` will be
internally converted to the nearest valid rate. The actual rate can be
read back using `mcc172_a_in_clock_config_read()`. When used in slave clock
configuration, the device will measure the frequency of the incoming
master clock after the synchronization period is complete. Calling
`mcc172_a_in_clock_config_read()` after this will return the measured data rate.

	Return: Result code, RESULT_SUCCESS if successful
Parameters:

	address: The board address (0 - 7). Board must already be opened.
	clock_source: The ADC clock source, one of the source type values.
	sample_rate_per_channel: The requested sample rate in samples per second per channel
"""
function mcc172_a_in_clock_config_write(address::Integer, clock_source::Symbol, sample_rate_per_channel::Real)
	# The clock can be: :SOURCE_LOCAL, :SOURCE_MASTER, OR :SOURCE_SLAVE
	# :SOURCE_LOCAL  = 0,
	# :SOURCE_MASTER = 1,
	# :SOURCE_SLAVE  = 2
	# Only one master device ot damage can occur
	# Configure clock on slave devices first and master last. This will cause the ADC's on all devices to be in sync.
	# If a clock is configured after the master the devices will no longer be in sync.
	# Slave devices must have a master or scans will not complete
	# Trigger must be used to start on same sample

	clock_source_num = source_convert(clock_source)
	mcc172_a_in_clock_config_write(address, clock_source_num, sample_rate_per_channel)
end

"""
	function mcc172_a_in_clock_config_write(address::Integer, clock_source::Integer, sample_rate_per_channel::Real)

Write the sampling clock configuration with `clock_source` being an integer 0, 1, or 2 as in help above.
"""
function mcc172_a_in_clock_config_write(address::Integer, clock_source::Integer, sample_rate_per_channel::Real)
	resultCode = ccall((:mcc172_a_in_clock_config_write, "libdaqhats.so"),
	Cint, (UInt8, UInt8, Cdouble), address, clock_source, sample_rate_per_channel)
	printError(resultCode)
	return resultCode
end

"""
	mcc172_trigger_config(address::Integer, source::Symbol, mode::Symbol)

Configure the digital trigger.

The analog input scan may be configured to start saving the acquired
data when the digital trigger is in the desired state. A single device
trigger may also be shared with multiple boards. This command sets the
trigger source and mode.

The trigger source must be one of:

	SOURCE_LOCAL (0): The trigger terminal on this MCC 172 is used and not shared with other MCC 172s.
	SOURCE_MASTER (1): The trigger terminal on this MCC 172 is used and is shared as the master trigger for other MCC 172s.
	SOURCE_SLAVE (2): The trigger terminal on this MCC 172 is not used, it receives its trigger from the master MCC 172.
The trigger mode must be one of:

	TRIG_RISING_EDGE: Start the scan on a rising edge of TRIG.
	TRIG_FALLING_EDGE: Start the scan on a falling edge of TRIG.
	TRIG_ACTIVE_HIGH: Start the scan any time TRIG is high.
	TRIG_ACTIVE_LOW: Start the scan any time TRIG is low.

Due to the nature of the filtering in the A/D converters there is an
input delay of 39 samples, so the data coming from the converters at any
time is delayed by 39 samples from the current time. This is most
noticeable when using a trigger - there will be approximately 39 samples
prior to the trigger event in the captured data.

Care must be taken when using master / slave triggering; the input
trigger signal on the master will be passed through to the slave(s),
but the mode is set independently on each device. For example, it is
possible for the master to trigger on the rising edge of the signal
and the slave to trigger on the falling edge.

	Return: Result code, RESULT_SUCCESS if successful.
Parameters:

	address: The board address (0 - 7). Board must already be opened.
	source: The trigger source, one of the source type values.
	mode: The trigger mode, one of the trigger mode values.
"""
function mcc172_trigger_config(address::Integer, source::Symbol, mode::Symbol)
	source_num = source_convert(source)
	mode_num = mode_convert(mode)
	# @show(source_num, mode_num)

	# ADC is pretriggered by 39 samples
	resultCode = mcc172_trigger_config(address, source_num, mode_num)
	return resultCode
end

"""
	function mcc172_trigger_config(address::Integer, source::Integer, mode::Integer)

Configure the digital trigger with source & mode as an Integers.
"""
function mcc172_trigger_config(address::Integer, source::Integer, mode::Integer)
	# ADC is pretriggered by 39 samples
	resultCode = ccall((:mcc172_trigger_config, "libdaqhats.so"),
	Cint, (UInt8, UInt8, UInt8), address, source, mode)
	printError(resultCode)
	return resultCode
end

"""
	mcc172_a_in_scan_start(address::Integer, channel_mask::UInt8, samples_per_channel::UInt32, options::UInt32)

Start capturing analog input data from the specified channels.

The scan runs as a separate thread from the user’s code. The function
will allocate a scan buffer and read data from the device into that buffer.
The user reads the data from this buffer and the scan status using the
`mcc172_a_in_scan_read()` function. `mcc172_a_in_scan_stop()` is used to stop
a continuous scan, or to stop a finite scan before it completes. The user
must call `mcc172_a_in_scan_cleanup()` after the scan has finished and all
desired data has been read; this frees all resources from the scan and
allows additional scans to be performed.

The scan cannot be started until the ADCs are synchronized, so this
function will not return until that has completed. It is best to wait
for sync using `mcc172_a_in_clock_config_read()` before starting the scan.

The scan state has defined terminology:

	Active: mcc172_a_in_scan_start() has been called and the device may be acquiring data or finished with the acquisition. The scan has not been cleaned up by calling mcc172_a_in_scan_cleanup(), so another scan may not be started.
	Running: The scan is active and the device is still acquiring data. Certain functions will return an error because the device is busy.
The valid options are:

	OPTS_NOSCALEDATA: Returns ADC code (a value between AI_MIN_CODE and AI_MAX_CODE) rather than voltage.
	OPTS_NOCALIBRATEDATA: Return data without the calibration factors applied.
	OPTS_EXTTRIGGER: Hold off the scan (after calling mcc172_a_in_scan_start()) until the trigger condition is met.
	OPTS_CONTINUOUS: Scans continuously until stopped by the user by calling mcc172_a_in_scan_stop() and writes data to a circular buffer. The data must be read before being overwritten to avoid a buffer overrun error. samples_per_channel is only used for buffer sizing.
	OPTS_EXTCLOCK: option is not supported for this device and will return an error.

The options parameter is set to 0 or OPTS_DEFAULT for default operation,
which is scaled and calibrated data, no trigger, and finite operation.

Multiple options may be specified by ORing the flags. For instance,
specifying `OPTS_NOSCALEDATA | OPTS_NOCALIBRATEDATA` will return the values
read from the ADC without calibration or converting to voltage.

The buffer size will be allocated as follows:

	Finite mode: Total number of samples in the scan
	Continuous mode (buffer size is per channel): Either samples_per_channel or the value in the following table, whichever is greater

Sample Rate	Buffer Size (per channel)

	200-1024 	S/s		1 kS
	1280-10.24 	kS/s	10 kS
	12.8, 25.6, 51.2 	kS/s	100 kS
Specifying a very large value for samples_per_channel could use too much
of the Raspberry Pi memory. If the memory allocation fails, the function
will return `RESULT_RESOURCE_UNAVAIL`. The allocation could succeed, but
the lack of free memory could cause other problems in the Raspberry Pi.
If you need to acquire a high number of samples then it is better to run
the scan in continuous mode and stop it when you have acquired the
desired amount of data. If a scan is already active this function will
return RESULT_BUSY.

Return:

	Result code,
		RESULT_SUCCESS if successful,
		RESULT_BUSY if a scan is already active.
Parameters:

	address: The board address (0 - 7). Board must already be opened.
	channel_mask: A bit mask of the channels to be scanned. Set each bit to enable the associated channel (0x01 - 0x03.)
	samples_per_channel: The number of samples to acquire for each channel in the scan (finite mode,) or can be used to set a larger scan buffer size than the default value (continuous mode.)
	options: The options bitmask.
"""
function mcc172_a_in_scan_start(address::Integer, channel_mask::UInt8, samples_per_channel::UInt32, options::Integer)

	# start capturing analog input from selected channels in separate thread
	# see documentation at https://mccdaq.github.io/daqhats/c.html#
	# Result code, RESULT_SUCCESS if successful, RESULT_BUSY if a scan is already active.
	resultCode = ccall((:mcc172_a_in_scan_start, "libdaqhats.so"),
	Cint, (UInt8, UInt8, UInt32, UInt32),
	address, channel_mask, samples_per_channel, options)
	printError(resultCode)
end

"""
	function mcc172_a_in_scan_start(address::Int32, channel_mask::UInt8, samples_per_channel::UInt32, options::Set{Symbol})

Put the options as symbols in a set and this program will perform the option masking and call the scan start program
Note options must be a `Set{Symbol}([:OPTS_DEFAULT, :OPTS_NOSCALEDATA, :OPTS_NOCALIBRATEDATA, :OPTS_EXTCLOCK, :OPTS_EXTTRIGGER, :OPTS_CONTINUOUS])`

	eg one value will be entered as Set([:OPTS_DEFAULT])
	two values as Set([:OPTS_NOCALIBRATEDATA; :OPTS_NOSCALEDATA]) etc.
"""
function mcc172_a_in_scan_start(address::Integer, channel_mask::UInt8, samples_per_channel::UInt32, options::Set{Symbol})
	# same name but keywords put in as variable args
	options_dict = Dict{Symbol, UInt32}(
	:OPTS_DEFAULT => 0x0000,         # Default behavior.
	:OPTS_NOSCALEDATA => 0x0001,     # Read / write unscaled data.
	:OPTS_NOCALIBRATEDATA => 0x0002, # Read / write uncalibrated data.
	:OPTS_EXTCLOCK => 0x0004,        # Use an external clock source.
	:OPTS_EXTTRIGGER => 0x0008,      # Use an external trigger source.
	:OPTS_CONTINUOUS => 0x0010)      # Run until explicitly stopped.

	# ensure option is valid
	options_set = Set{Symbol}([:OPTS_DEFAULT, :OPTS_NOSCALEDATA, :OPTS_NOCALIBRATEDATA,
	:OPTS_EXTCLOCK, :OPTS_EXTTRIGGER, :OPTS_CONTINUOUS])
	# @show(options, options_set)
	if !issubset(options, options_set)
		error("options must be a subset of $options_set")
	end

	optionmask = 0x0000
	for i in options
		optionmask = optionmask | options_dict[i]
	end
	mcc172_a_in_scan_start(address, channel_mask, samples_per_channel, optionmask)
end

"""
	mcc172_a_in_scan_buffer_size(address::Integer)

Returns the size of the internal scan data buffer.

An internal data buffer is allocated for the scan when
`mcc172_a_in_scan_start()` is called. This function returns the total
size of that buffer in samples.

Return:

	Result code,
		RESULT_SUCCESS if successful,
		RESULT_RESOURCE_UNAVAIL if a scan is not currently active,
		RESULT_BAD_PARAMETER if the address is invalid or buffer_size_samples is NULL.
Parameters:

	address: The board address (0 - 7). Board must already be opened.
	buffer_size_samples: Receives the size of the buffer in samples. Each sample is a double.
"""
function mcc172_a_in_scan_buffer_size(address::Integer)
	# Returns the size of the internal scan data buffer
	# Result code,
	# RESULT_SUCCESS if successful,
	# RESULT_RESOURCE_UNAVAIL if a scan is not currently active,
	# RESULT_BAD_PARAMETER if the address is invalid or buffer_size_samples is NULL.
	buffer_size_samples = Ref{UInt32}()
	resultCode = ccall((:mcc172_a_in_scan_buffer_size, "libdaqhats.so"),
	Cint, (UInt8, Ref{UInt32}), address, buffer_size_samples)
	printError(resultCode)
	return buffer_size_samples[]
end

"""
	mcc172_a_in_scan_status(address::Integer)
Reads status and number of available samples from an analog input scan.

The scan is started with `mcc172_a_in_scan_start()` and runs in a
background thread that reads the data from the board into an internal
scan buffer. This function reads the status of the scan and amount of
data in the scan buffer.

Input Parameter:

	address: The board address (0 - 7). Board must already be opened.

Return:

	Result code,
		RESULT_SUCCESS if successful,
		RESULT_RESOURCE_UNAVAIL if a scan has not been started under this instance of the device.
	status: Receives the scan status, an ORed combination of the flags (use mcc172_status_decode() for Array of codes):
		STATUS_HW_OVERRUN: The device scan buffer was not read fast enough and data was lost.
		STATUS_BUFFER_OVERRUN: The thread scan buffer was not read by the user fast enough and data was lost.
		STATUS_TRIGGERED: The trigger conditions have been met.
		STATUS_RUNNING: The scan is running.
	samples_per_channel: Receives the number of samples per channel available in the scan thread buffer.
"""
function mcc172_a_in_scan_status(address::Integer)
	# Reads status and number of available samples from an analog input scan.
	# Status is an || combination of flags:
	# STATUS_HW_OVERRUN (0x0001) A hardware overrun occurred.
	# STATUS_BUFFER_OVERRUN (0x0002) A scan buffer overrun occurred.
	# STATUS_TRIGGERED (0x0004) The trigger event occurred.
	# STATUS_RUNNING   (0x0008) The scan is running (actively acquiring data.)

	status = Ref{UInt16}()			# Initialize
	samples_per_channel_available = Ref{UInt32}()	# Initialize

	resultCode = ccall((:mcc172_a_in_scan_status, "libdaqhats.so"),
	Cint, (UInt8, Ref{UInt16}, Ref{UInt32}),
	address, status, samples_per_channel_available)

	printError(resultCode)
	return (resultCode, status[], samples_per_channel_available[])
end

"""
	mcc172_a_in_scan_read(address::Integer, samples_per_channel::UInt32, mcc172_num_channels::Integer, timeout::Float64)

Reads status and multiple samples from an analog input scan.

The scan is started with `mcc172_a_in_scan_start()` and runs in a
background thread that reads the data from the board into an internal
scan buffer. This function reads the data from the scan buffer,
and returns the current scan status.

Input Parameters:

	address: The board address (0 - 7). Board must already be opened.
	samples_per_channel: The number of samples per channel to read. Specify -1 to read all available samples in the scan thread buffer, ignoring timeout. If buffer does not contain enough space then the function will read as many samples per channel as will fit in buffer.
	mcc172_num_channels: number of channels required to calculate buffer size
	timeout: The amount of time in seconds to wait for the samples to be read. Specify a negative number to wait indefinitely or 0 to return immediately with whatever samples are available (up to the value of samples_per_channel or buffer_size_samples.)

Return Parameters:

	Result code:
		RESULT_SUCCESS if successful,
		RESULT_RESOURCE_UNAVAIL if a scan is not active.
	status: Receives the scan status, an ORed combination of the flags (use mcc172_status_decode() for Array of codes):
		STATUS_HW_OVERRUN: The device scan buffer was not read fast enough and data was lost.
		STATUS_BUFFER_OVERRUN: The thread scan buffer was not read by the user fast enough and data was lost.
		STATUS_TRIGGERED: The trigger conditions have been met.
		STATUS_RUNNING: The scan is running.
	buffer: The user data buffer that receives the samples.
	buffer_size_samples: The size of the buffer in samples. Each sample is a double.
	samples_read_per_channel: Returns the actual number of samples read from each channel.
	"""
function mcc172_a_in_scan_read(address::Integer, samples_per_channel::Int32, mcc172_num_channels::Integer, timeout::Real)
	# Reads status and number of available samples from an analog input scan.
	# Status is an || combination of flags:
	# STATUS_HW_OVERRUN (0x0001) A hardware overrun occurred.
	# STATUS_BUFFER_OVERRUN (0x0002) A scan buffer overrun occurred.
	# STATUS_TRIGGERED (0x0004) The trigger event occurred.
	# STATUS_RUNNING   (0x0008) The scan is running (actively acquiring data.)

	status = Ref{UInt16}()					# Initialize
	#@show(samples_per_channel)
	#@infiltrate
	if samples_per_channel == -1
		result_code, status, samples_per_channel = mcc172_a_in_scan_status(address)
	end
	#@show(samples_per_channel)
	buffer_size_samples::Int32 = samples_per_channel * mcc172_num_channels
	buffer = Vector{Float64}(undef, buffer_size_samples)
	samples_read_per_channel = Ref{UInt32}() # Initialize
	#@show(status, buffer_size_samples, buffer, samples_per_channel)
	resultCode = ccall((:mcc172_a_in_scan_read, "libdaqhats.so"),
	Cint, (UInt8, Ref{UInt16}, UInt32, Cdouble, Ptr{Cdouble}, UInt32, Ref{UInt32}),
	address, status, samples_per_channel, timeout, buffer, buffer_size_samples, samples_read_per_channel)
	printError(resultCode)
	return resultCode, status[], buffer, Int(samples_read_per_channel[])
end

"""
	mcc172_a_in_scan_channel_count(address::Integer)

Return the number of channels in the current analog input scan.

This function returns 0 if no scan is active.

Return: The number of channels, 0 - 2.

Parameters: address: The board address (0 - 7). Board must already be opened.
"""
function mcc172_a_in_scan_channel_count(address::Integer)
	numChannels = ccall((:mcc172_a_in_scan_channel_count, "libdaqhats.so"),
	Cint, (UInt8,), address)
	return numChannels
end

"""
	mcc172_a_in_scan_stop(address::Integer)

Stops an analog input scan.

The scan is stopped immediately. The scan data that has been read into
the scan buffer is available until `mcc172_a_in_scan_cleanup()` is called.
"""
function mcc172_a_in_scan_stop(address::Integer)
	resultCode = ccall((:mcc172_a_in_scan_stop, "libdaqhats.so"),
	Cint, (UInt8,), address)
	printError(resultCode)
	return resultCode
end

"""
	mcc172_a_in_scan_cleanup(address::Integer)

Stops an analog input scan.

The scan is stopped immediately. The scan data that has been read into
the scan buffer is available until `mcc172_a_in_scan_cleanup()` is called.
"""
function mcc172_a_in_scan_cleanup(address::Integer)
	resultCode = ccall((:mcc172_a_in_scan_cleanup, "libdaqhats.so"),
	Cint, (UInt8,), address)
	printError(resultCode)
	return resultCode
end

max_stars_repo_name: Spectrum-Tec/MccDaqHats

max_issues_repo_name: Spectrum-Tec/MccDaqHats

max_forks_repo_name: Spectrum-Tec/MccDaqHats


text:
# This file was generated, do not modify it. # hide
DecisionTreeClassifier = @load DecisionTreeClassifier pkg=DecisionTree
tree_model = DecisionTreeClassifier()

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
#=
# A random node reordering.
=#
export reorder_grid_random, random_order

function random_order(N, seed)
    rng = Random.MersenneTwister(seed); # Use a specific RNG seed to make results consistent
    return randperm(rng, N); # Randomly permute node order
end

function reorder_grid_random(grid, seed = 17)
    dim = size(grid.allnodes,1);
    nnodes = size(grid.allnodes,2);
    nel = size(grid.loc2glb,2);
    nfc = size(grid.face2glb,2);

    norder= random_order(nnodes, seed);

    newnodes = zeros(size(grid.allnodes));
    newbdry = copy(grid.bdry);
    newloc2glb = copy(grid.loc2glb);
    newglbvertex = copy(grid.glbvertex);
    newface2glb = copy(grid.face2glb);
    newfaceVertex2glb = copy(grid.faceVertex2glb);

    for mi=1:length(norder)
        for d=1:dim
            newnodes[d,norder[mi]] = grid.allnodes[d,mi];
        end
    end

    for bid=1:length(newbdry)
        for i=1:length(newbdry[bid])
            newbdry[bid][i] = norder[grid.bdry[bid][i]];
        end
    end

    for ei=1:nel
        for ni=1:size(newloc2glb,1)
            newloc2glb[ni,ei] = norder[grid.loc2glb[ni,ei]];
        end
        for ni=1:size(newglbvertex,1)
            newglbvertex[ni,ei] = norder[grid.glbvertex[ni,ei]];
        end
    end

    for fi=1:nfc
        for ni=1:size(newface2glb,1)
            newface2glb[ni,fi] = norder[grid.face2glb[ni,fi]];
        end
        for ni=1:size(newfaceVertex2glb,1)
            newfaceVertex2glb[ni,fi] = norder[grid.faceVertex2glb[ni,fi]];
        end
    end

    return Femshop.Grid(newnodes, newbdry, grid.bdryface, grid.bdrynorm, grid.bids, newloc2glb, newglbvertex, newface2glb, newfaceVertex2glb);
end


text:
"""
    sample_names(x::PopObj)
View individual/sample names in a `PopObj`

Equivalent to `PopObj.samples.name`
"""
sample_names(x::PopObj) = x.samples.name

"""
    summary(x::PopObj)
Prints a summary of the information contained in a PopObj
"""
function Base.summary(x::PopObj)
    y = PopOpt(x)
    println("Object of type $(typeof(x)):")
    println("\nLongitude:")
    println("$(y.samples.longitude[1:3])" , " \u2026 ", "$(y.samples.longitude[end-2:end])", "\n")
    println("Latitude:")
    println("$(y.samples.latitude[1:3])" , " \u2026 ", "$(y.samples.latitude[end-2:end])", "\n")
    println("Number of individuals: $(length(y.samples.name))")
    println(y.samples.name[1:3] , " \u2026 ", y.samples.name[end-2:end], "\n")
    println("Number of loci: $(size(y.loci,2))")
    println(string.(names(y.loci))[1:3], " \u2026 " , string.(names(y.loci))[end-2:end], "\n" )
    println("Ploidy:")
    println("$(y.samples.ploidy[1:3])", " \u2026 ", "$(y.samples.ploidy[end-2:end])", "\n")
    println("Population names and counts:")
    print(populations(x), "\n")
    println("\nAvailable .samples fields: .name, .population, .ploidy, .longitude, .latitude")
end

"""
    isolate_genotypes(x::PopObj; samples::Union{String, Array, Nothing}, loci::Union{String, Array, Nothing})
View the genotypes of specific samples for specific loci in a `PopObj`.
Default shows all genotypes for all individuals.

`isolate_genotypes(nancycats, loci = "fca8")`

`isolate_genotypes(nancycats, samples = "N226", loci = ["fca8", "fca23"])`
"""
function isolate_genotypes(x::PopObj; samples::Union{String, Array, Nothing}= nothing, loci::Union{String, Array, Nothing}= nothing)
    if loci == nothing && samples == nothing
        @warn "please specify either loci= or samples=, otherwise use PopObj.loci"
    end
    y = PopOpt(x)
    df = deepcopy(y.loci)
    insertcols!(df, 1, :name => y.samples.name) ;
    insertcols!(df, 2, :population => categorical(y.samples.population))
    if samples != nothing
        if typeof(samples) == String
            samples ∉ y.samples.name && error("individual $samples not found in PopObj")
            tmp = df[df.name .== samples, :]
        else
            tmp = df[df.name .== samples[1], :]
            for ind in samples[2:end]
                ind ∉ y.samples.name && println("NOTICE: individual \"$ind\" not found in PopObj!")
                tmp = vcat(tmp, df[df.name .== ind, :])
            end
            println()
        end
    else
        tmp = df
    end
    if loci != nothing
        if typeof(loci) == String
            loci ∉ string.(names(y.loci)) && error("locus $loci not found in PopObj")
            return tmp[!, [:name, :population, Symbol(loci)]]
        else
            for locus in loci[2:end]
                locus ∉ string.(names(y.loci)) && println("NOTICE: locus \"$locus\" not found in PopObj!")
            end
            println()
            return tmp[!, append!([:name, :population], Symbol.(loci))]
        end
    else
        return tmp
    end
end

"""
    locations(x::PopObj)
View location data (`.longitude` and `.latitude`) in a `PopObj`

Use `locations!` to add spatial data to a `PopObj`
"""
function locations(x::PopObj)
    y = PopOpt(x)
    DataFrame(name = y.samples.name,
              population = y.samples.population,
              longitude = y.samples.longitude,
              latitude = y.samples.latitude)
end

"""
    locations!(x::PopObj; lat::Array, long::Array)
Add location data (latitude `lat`, longitude `long`) to `PopObj`. Takes decimal
degrees or decimal minutes format. **Must** use `-` symbol instead of cardinal directions.
Location data must be in order of `ind`. Replaces existing `PopObj` location data.
- Decimal Degrees : `-11.431`
- Decimal Minutes : `"-11 43.11"` (must use space and double-quotes)

If conversion is not necessary, can directly assign `PopObj.samples.longitude` and `PopObj.samples.latitude`
"""
function locations!(x::PopObj; lat::Array, long::Array)
    # test for decimal degrees vs decimal minutes
    if occursin(" ", string(lat[1])) == false && occursin(" ", string(long[1])) == false
        x.samples.longitude = lat ;
        x.samples.latitude = long ;
    else
        # make sure decimal minutes are Strings
        if typeof(lat) != Array{String,1}
            lat = string.(lat)
        end
        if typeof(long) != Array{String,1}
            long = string.(long)
        end
        # convert lat to decimal degrees
        latConverted = []
        for value in lat
            tmp = split(value, " ")
            if parse(Float64,tmp[1]) < 0   # if negative, subtract
                decideg = parse(Float64, tmp[1]) - round((parse(Float64,tmp[2])/60), digits = 3)
                push!(latConverted, decideg)
            else                           # if positive, add
                decideg = parse(Float64, tmp[1]) + round((parse(Float64,tmp[2])/60), digits = 3)
                push!(latConverted, decideg)
            end
        end
        # convert long to decimal degrees
        longConverted = []
        for value in long
            tmp = split(value, " ")
            if parse(Float64,tmp[1]) < 0
                decideg = parse(Float64, tmp[1]) - round((parse(Float64,tmp[2])/60), digits = 3)
                push!(longConverted, decideg)
            else
                decideg = parse(Float64, tmp[1]) + round((parse(Float64,tmp[2])/60), digits = 3)
                push!(longConverted, decideg)
            end
        end
        x.samples.longitude = latConverted
        x.samples.latitude = longConverted
        return x.samples

    end
end

"""
    populations(x::PopObj; listall::Bool = false)
View unique population ID's in a `PopObj`.

`listall = true`, displays `ind` and their `population` instead (default = `false`).
"""
function populations(x::PopObj; listall::Bool = false)
    y = PopOpt(x)
    if listall == true
        DataFrame(name = y.samples.name, population = y.samples.population)
    else
        count = [sum(y.samples.population .== i) for i in unique(y.samples.population)]
        count_conv = Int32.(count)
        popcounts = DataFrame(population = unique(y.samples.population) |> categorical,
                              count = count_conv)
    end
end

"""
    populations!(x::PopObj; rename::Dict, replace::Union{Tuple, NamedTuple})
Assign population names to a `PopObj`. There are two modes of operation:

## Rename

Rename the population ID's of `PopObj.population` using a `Dict` of `[population] => replacement`

Example:

`potatopops = Dict(1 => "Idaho", 2 => "Russet")`

`populations!(potatoes, rename = potatopops)`

## Replace (overwrite)
Completely replace the population names of a `PopObj`. Will generate an array of
population names from a tuple of (counts, names) where `counts` is an array of the
number of samples per population and `names` is an array of the names of the
populations. Can also use a named tuple.

Example assigning names for three populations in a `PopObj` named "Starlings":
Assuming population sizes are 15, 32, 11 and we want to name them "North", "South", "East"

`populations!(Starlings, replace = ([15,32,11], ["North","South", "East"])`
`populations!(Starlings, replace = (counts = [15,32,11], names = ["North","South", "East"])`

"""
function populations!(x::PopObj; rename::Dict=Dict(), replace::Union{Tuple,NamedTuple}=(0,0))
    if length(keys(rename)) != 0
        for eachkey in keys(rename)
            eachkey ∉ x.samples.population && @warn "$eachkey not found in PopObj"
            replace!(x.samples.population, eachkey => rename[eachkey])
        end
        return populations(x,listall = true)
    elseif replace != (0,0)
        if typeof(replace) <: NamedTuple
            popid_array = [fill(i, j) for (i,j) in zip(replace.names, replace.counts)]
        else typeof(replace) <: Tuple
            popid_array = [fill(i, j) for (i,j) in zip(replace[2], replace[1])]
        end
        flat_popid = Iterators.flatten(popid_array) |> collect
        length(flat_popid) != size(x.samples, 1) && error("length of names ($(length(flat_popid))) does not match sample number ($(length(x.samples.name)))")
        x.samples.population = flat_popid
        return populations(x, listall = true)
    else
        error("specify rename = Dict() to rename populations, or replace=(counts,names) to replace all names")
    end
end

const population = populations
const population! = populations!

#### Find missing ####

"""
    missing(x::PopObj)
Identify and count missing loci in each sample of a `PopObj`. Returns a tuple
of `DataFrames`: loci per sample, number per loci.

Example:

`missing_ind,missing_loc = missing(gulfsharks)`
"""
function Base.missing(x::PopObj)
    y = PopOpt(x)
    df = deepcopy(y.loci)
    insertcols!(df, 1, :ind => y.samples.name)
    # missing per individual
    nmissing = []
    missing_array = []
    for each in 1:length(df[:,1])
        miss_idx = findall(i -> i === missing, df[each,:])
        push!(nmissing, miss_idx |> length)
        push!(missing_array, String.(miss_idx))
    end
    sample_df = DataFrame(name = x.samples.name,
                       population = x.samples.population,
                       missing = nmissing,
                       loci = missing_array
                       )
    # missing per locus
    f(x) = map(eachcol(x)) do col
        count(i->i===missing, col)
    end

    loci_df = DataFrame(locus = string.(names(y.loci)), missing = f(y.loci))
    return (by_sample = sample_df, by_loci = loci_df)
end

##### Removal #####

"""
    remove_loci!(x::PopObj, loci::Union{String, Array{String,1}})
Removes selected loci from a `PopObj`.

Examples:

`remove_loci!(nancycats, "fca8")`

`remove_loci!(tulips, ["fca8", "fca23"])`
"""
function remove_loci!(x::PopObj, loci::Union{String,Array{String,1}})
    y = PopOpt(x)
    # get individuals indices
    if typeof(loci) == String
        sym_loci = Symbol(loci)
        sym_loci ∉ names(y.loci) && error("Locus \"$loci\" not found")
        return select!(x.loci, Not(sym_loci))
    else
        sym_loci = Symbol.(loci)
        for each in sym_loci
            if each ∉ names(y.loci)
                println("NOTICE: locus \"$each\" not found")
                continue
            end
        end
        return select!(x.loci, Not(sym_loci))
    end
end

"""
    remove_samples!(x::PopObj, samp_id::Union{Array{String,1}})
Removes selected samples from a `PopObj`.

Examples:

`remove_samples!(nancycats, "N100")`

`remove_samples!(nancycats, ["N100", "N102", "N211"])`
"""
function remove_samples!(x::PopObj, samp_id::Union{String, Array{String,1}})
    y = PopOpt(x)
    # get samp_id indices
    if typeof(samp_id) == String
        samp_id ∉ y.samples.name && error("sample \"$samp_id\" not found")
        idx = findfirst(i -> i == samp_id, y.samples.name)
    else
        idx = []
        for ind in samp_id
            if ind ∉ y.samples.name
                println("NOTICE: sample \"$ind\" not found!")
                continue
            end
            push!(idx, findfirst(i -> i == ind, y.samples.name) )
        end
        println()
    end
    deleterows!(x.samples, idx)
    deleterows!(x.loci, idx)
    return x
end


text:
using ECAEstimator
using Test

@testset "ECAEstimator.jl" begin
    # Write your tests here.
end

max_stars_repo_name: giannamars/ECAEstimator

max_issues_repo_name: giannamars/ECAEstimator

max_forks_repo_name: giannamars/ECAEstimator


text:
## ExoplanetsSysSim/examples/syssim_dist.jl
## (c) 2015 Eric B. Ford

# Calculaties the distance between simulation computed from an input "parameter file" and the summary statistics already stored in comp_file
# usage:  syssim_dist.jl param_file comp_file
# output: scalar_dist d_1 d_2 ...
#         where scalar_dist is the scalar distance between the simulated model's summary statistics and the precomputed summary statistics,
#         and d_1 d_2 ... represent a list of the individual distance metrics between the two sets of summary statistics

verbose = true

if length(ARGS)<2
   warn("usage:  syssim_dist.jl param_file comp_file")
   exit(255)
end

using Compat
import Compat: UTF8String, ASCIIString

param_file = convert(Compat.ASCIIString,ARGS[1])
if !isfile(param_file)
   warn(string("Can't read ",param_file))
   exit(255)
end

comp_file = convert(Compat.ASCIIString,ARGS[2])
if !isfile(comp_file)
   warn(string("Can't read ",comp_file))
   exit(255)
end

if verbose println("# Loading ExoplanetsSysSim module."); end
using ExoplanetsSysSim

if verbose println("# Reading parameter file ", param_file,"."); end
include(param_file)
if verbose println("# Calling setup_sim_param( ", ARGS[3:end], " )."); end
sim_param = setup_sim_param( convert(Array{Compat.ASCIIString,1},ARGS[3:end]) )
if verbose println("# Active parameter values: ", make_vector_of_sim_param(sim_param) ); end

#add_param_fixed(sim_param,"rng_seed",1234)   # If you want to be able to reproduce simulations
if haskey(sim_param,"rng_seed")
   seed = get_int(sim_param,"rng_seed")
     if verbose println("# Seeding RNG with seed = ", seed); end
   srand(seed)
end

if verbose println("# Loading summary statistics from ", comp_file,"."); end
ss_comp = load_summary_stats(comp_file)

num_repeats = 1
for k in 1:num_repeats
if num_repeats > 1  print("\n\n"); end
if verbose println("# Generating physical catalog."); end
@time cat_phys = generate_kepler_physical_catalog(sim_param)
@time cat_phys = generate_kepler_physical_catalog(sim_param)

if verbose println("# Generating observational catalog."); end
@time cat_obs = observe_kepler_targets_single_obs(cat_phys,sim_param)
@time cat_obs = observe_kepler_targets_single_obs(cat_phys,sim_param)

ss = CatalogSummaryStatistics()
if haskey(sim_param,"use obs data only") &&  convert(Bool,get(sim_param,"use obs data only"))
  if verbose println("# Calculating summary statistics based on all data."); end
  @time ss = calc_summary_stats_obs_demo(cat_obs,sim_param)
else
  if verbose println("# Calculating summary statistics based on observational data only."); end
  @time ss = calc_summary_stats_sim_pass_one_demo(cat_obs,cat_phys,sim_param)
  @time ss = calc_summary_stats_sim_pass_two_demo(cat_obs,cat_phys,ss,sim_param)
end

if verbose println("# Calculating distance (pass 1)"); end
@time d1 = calc_distance_vector_demo(ss_comp, ss, 1, sim_param)
if verbose println("# Calculating distance (pass 2)"); end
@time d2 = calc_distance_vector_demo(ss_comp, ss, 2, sim_param)

d = [d1; d2]
@time dist = calc_scalar_distance(d)
print(dist, "   ")
print(STDOUT,join([d1;d2],' '))
print("\n")
end


text:
# T: coeff type
# Degree: degree
# Dim1: dimension of the support + 1
type LagrangeRefSpace{T,Degree,Dim1,NF} <: RefSpace{T,NF} end

numfunctions{T,D}(s::LagrangeRefSpace{T,D,2}) = D+1
numfunctions{T}(s::LagrangeRefSpace{T,0,3}) = 1
numfunctions{T}(s::LagrangeRefSpace{T,1,3}) = 3

valuetype{T}(ref::LagrangeRefSpace{T}, charttype) =
        SVector{numfunctions(ref), Tuple{T,T}}

# Evaluate constant lagrange elements on anything
(ϕ::LagrangeRefSpace{T,0}){T}(tp) = SVector(((one(T),zero(T)),))

# Evaluate linear Lagrange elements on a segment
function (f::LagrangeRefSpace{T,1,2}){T}(mp)
    u = mp.bary[1]
    j = jacobian(mp)
    SVector((u,-1/j), (1-u,1/j))
end

# Evaluete linear lagrange elements on a triangle
function (f::LagrangeRefSpace{T,1,3}){T}(t)
    u,v,w, = barycentric(t)
    SVector(u, v, w)
end

"""
    f(tangent_space, Val{:withcurl})

Compute the values of the shape functions together with their curl.
"""
function (f::LagrangeRefSpace{T,1,3}){T}(t, ::Type{Val{:withcurl}})
    # Evaluete linear Lagrange elements on a triange, together with their curl
    j = jacobian(t)
    u,v,w, = barycentric(t)
    p = t.patch
    SVector(
        (u, (p[3]-p[2])/j),
        (v, (p[1]-p[3])/j),
        (w, (p[2]-p[1])/j)
    )
end

# Evaluate constant Lagrange elements on a triangle, with their curls
function (f::LagrangeRefSpace{T,0,3}){T}(t, ::Type{Val{:withcurl}})
    i = one(T)
    z = zero(cartesian(t))
    (
        (i,z,),
    )
end

function curl(ref::LagrangeRefSpace, sh, el)
    sh1 = Shape(sh.cellid, mod1(sh.refid+1,3), -sh.coeff)
    sh2 = Shape(sh.cellid, mod1(sh.refid+2,3), +sh.coeff)
    return [sh1, sh2]
end

function strace(x::LagrangeRefSpace, cell, localid, face)

    Q = zeros(scalartype(x),2,3)

    p1 = neighborhood(face, 1)
    p2 = neighborhood(face, 0)

    u1 = carttobary(cell, cartesian(p1))
    u2 = carttobary(cell, cartesian(p2))

    P1 = neighborhood(cell, u1)
    P2 = neighborhood(cell, u2)

    vals1 = x(P1)
    vals2 = x(P2)

    for j in 1:numfunctions(x)
        Q[1,j] = vals1[j]
        Q[2,j] = vals2[j]
    end

    Q
end

function restrict{T}(refs::LagrangeRefSpace{T,0}, dom1, dom2)
    Q = eye(T, numfunctions(refs))
end

function restrict{T}(f::LagrangeRefSpace{T,1}, dom1, dom2)

    D = numfunctions(f)
    Q = zeros(T, D, D)

    # for each point of the new domain
    for i in 1:D
        v = dom2.vertices[i]

        # find the barycentric coordinates in dom1
        uvn = carttobary(dom1, v)

        # evaluate the shape functions in this point
        x = neighborhood(dom1, uvn)
        fx = f(x)

        for j in 1:D
            Q[j,i] = fx[j][1]
        end
    end

    return Q
end


text:
function star1(n_wait)
    scores = [3, 7]

    elf1 = 1
    elf2 = 2

    while length(scores) < 10 + n_wait
        new_score = scores[elf1] + scores[elf2]
        scores = append!(scores, reverse!(digits(new_score)))

        elf1 = mod(elf1 + scores[elf1], length(scores)) + 1
        elf2 = mod(elf2 + scores[elf2], length(scores)) + 1
    end

    print.(scores[ (1:10) .+ n_wait ])
    println()

end

function star2(input::String)
    scores = [3, 7]

    elf1 = 1
    elf2 = 2

    looking_for = parse.(Int, collect(input))
    n_check = length(looking_for)
    done = false

    while !done
        new_score = scores[elf1] + scores[elf2]

        for new_recipe in reverse!(digits(new_score))
            push!(scores, new_recipe)
            # check if last digits have come up
            if length(scores) >= n_check && scores[end-n_check+1:end] == looking_for
                # println("Found $input in $scores")
                return length(scores)-n_check
            end
        end

        elf1 = mod(elf1 + scores[elf1], length(scores)) + 1
        elf2 = mod(elf2 + scores[elf2], length(scores)) + 1

    end
end


text:
# Test applycolourmaps.jl

# Tests could be more comprehensive (of course).  In some cases all
# that is done is to execute the function in some way just to make
# sure no exceptions are thrown as a result of something breaking with
# a new version of Julia.

println("testing applycolourmaps")

# applycolourmap
img = sineramp()
(lmap,) = cmap("L1")
rnge = [10, 200]
rgbimg = applycolourmap(img, lmap)
rgbimg = applycolourmap(img, lmap, rnge)

# applycycliccolourmap
(ang, mask) = circlesineramp()
(cycmap,) = cmap("C1")
rgbimg = applycycliccolourmap(ang, cycmap)
rgbimg = applycycliccolourmap(ang, cycmap, amp = ang, cyclelength = pi)

# applydivergingcolourmap
(divmap,) = cmap("D1")
refval = 0
rgbim = applydivergingcolourmap(img, divmap, refval)

# ternaryimage
bands = [2, 3, 1]
histcut = 2
(rows,cols) = size(img)
cimg = zeros(rows,cols,3)
for n = 1:3
    cimg[:,:,n] = img
end
rgbimg = ternaryimage(cimg, bands=bands, histcut=histcut)
rgbimg = ternaryimage(cimg, bands=bands, histcut=histcut, RGB=false)