text:
module SimpleNeuralNetwork

import Base: repr, show

export NeuralLayer, NeuralNetwork, train!, predict, repr, show

function derivative(fun::Function; epsilon = 1e-3)
    function (x)
        (fun(x + epsilon) - fun(x - epsilon)) / (2*epsilon)
    end
end

function sqr_cost(y_predict::Array, y::Array)
    sqrt(sum((y_predict - y).^2))
end

type NeuralLayer
    weight::Union(Matrix, Array{None, 1}) # Let conversion do the work....
    _nodes_value::Array{Float64, 1}

    function NeuralLayer(w, nv)
        new(w, nv)
    end
end

type NeuralNetwork
    structure::Array{Int, 1}
    act_fun::Function
    act_diff::Function
    layers::Array{NeuralLayer, 1}
    out_fun::Function
    out_diff::Function
    _deltas::Array{Vector{FloatingPoint}, 1} # Array for storing dels

    # Constructor for NeuralNetwork type
    function NeuralNetwork(struct; act_fun = tanh, act_diff = None, out_fun = x-> x, out_diff = None)
        if act_diff == None
            act_diff = derivative(act_fun)
        end
        if out_diff == None
            out_diff = derivative(out_fun)
        end
        layers = NeuralLayer[]
        for ind = 1:(length(struct)-1)
            dim_in = struct[ind]
            dim_out = struct[ind+1]
            b = sqrt(6) / sqrt(dim_in + dim_out)
            w = 2b*rand(dim_out, dim_in + 1) - b
            nodes_value = push!([1.0], rand(dim_in)...)
            temp_layer = NeuralLayer(w, nodes_value)
            push!(layers, temp_layer)
        end
        # append the output layer.
        w = []
        nodes_value = rand(struct[end])
        temp_layer = NeuralLayer(w, nodes_value)
        push!(layers, temp_layer)
        _deltas = Vector{Float64}[[1.] for i = 1:length(struct)]
        nn = new(struct, act_fun, act_diff, layers, out_fun, out_diff, _deltas)
    end
end

function repr(nn::NeuralNetwork)
    struct = join([string(i) for i in nn.structure], "x")
    msg = join(["It is a ", struct, " NeuralNetwork.\n"] , "")
    msg = join([msg, "Activate Function: ", string(nn.act_fun), '\n'], "")
    msg = join([msg, "Output Function: ", string(nn.out_fun), '\n'], "")
    msg
end

function show(nn::NeuralNetwork)
    print(repr(nn))
    println()
end

function predict(nn::NeuralNetwork, data::Array{Float64, 2})
    predict_results = Array{Float64, 1}[]
    for data_id = 1:size(data)[1]
        v = data[data_id, :][:]
        #println(v)
        forward_prob!(nn, v)
        push!(predict_results, nn.out_fun(nn.layers[end]._nodes_value))
    end
    return predict_results
end

function forward_prob!(nn::NeuralNetwork, x)
    # forward_prob! will update nodes_value for all layers.
    nn.layers[1]._nodes_value = [1, x]
    n_layers = length(nn.structure)
    for layer_id = 1:(n_layers-2)
        current_layer = nn.layers[layer_id]
        next_layer = nn.layers[layer_id + 1]
        temp = current_layer.weight * current_layer._nodes_value
        next_layer._nodes_value = [1., nn.act_fun(temp)]
    end
    # Compute the node values of the last layer without pass through activation function.
    current_layer = nn.layers[end - 1]
    next_layer = nn.layers[end]
    temp = current_layer.weight * current_layer._nodes_value
    next_layer._nodes_value = temp[:]
    return
end

function back_prob!(nn::NeuralNetwork, x, y)
    # back_prob! will update nn._detas.
    forward_prob!(nn, x)
    nn._deltas[end] = -2.*(y - nn.out_fun(nn.layers[end]._nodes_value)).*nn.out_diff(nn.layers[end]._nodes_value)

    n_layers = length(nn.structure)
    for layer_id = n_layers-1:-1:1
        delta_next = nn._deltas[layer_id + 1]
        w_this = nn.layers[layer_id].weight[:, 2:end]
        nodes_value_this = nn.layers[layer_id]._nodes_value
        dd = nn.act_diff(nodes_value_this[2:end])
        temp = transpose(w_this) * delta_next
        nn._deltas[layer_id] = temp[:] .* dd
    end
    return
end

function train!(nn::NeuralNetwork, X::Matrix, Y::Array; epos = 10000, cost_fun = sqr_cost, tol = 0.0001, learning_rate = 0.1)
    y_predict = predict(nn, X)
    #return y_predict
    err_p = cost_fun(y_predict, Y)
    #return err_p
    n_obs = size(X)[1]
    #return n_obs
    errors = Float64[]
    for iter = 1:epos
        ind = rand(1:n_obs)
        x = X[ind, :][:]
        y = Y[ind]
        back_prob!(nn, x, y)
        #return nn
        for layer_id in 1:(length(nn.layers)-1)
            next_id = layer_id + 1
            nl = nn.layers[layer_id]
            gradient = nn._deltas[next_id]*nl._nodes_value'
            nl.weight -= learning_rate .* gradient
        end

        y_predict = predict(nn, X)
        err_t = cost_fun(y_predict, Y)
        push!(errors, err_t)
        if iter % 1000 == 0
            if err_p - err_t > 0 && err_p - err_t < tol
                println("Terminating training process due to no significant improvement.")
                println("At iteration No. ", iter)
                break
            end
        end
    end
    return errors
end

end # module


text:
Base.:~(lhs::Number, rhs::Sym) = ~(promote(lhs, rhs)...)
Base.:~(lhs::Sym, rhs::Number) = ~(promote(lhs, rhs)...)
Base.:~(lhs::Sym, rhs::Sym) = asSymbolic(sympy.py.Eq(lhs.py, rhs.py))

ext: jl

max_stars_repo_name: jverzani/SymPyCall

max_issues_repo_name: jverzani/SymPyCall

max_forks_repo_name: jverzani/SymPyCall


text:
module Strapping

using Tables, StructTypes
# using DebugPrinting

struct Error <: Exception
    msg::String
end

"""
    Strapping.construct(T, tbl)
    Strapping.construct(Vector{T}, tbl)

Given a [Tables.jl](https://github.com/JuliaData/Tables.jl)-compatible input table source `tbl`,
construct an instance of `T` (single object, first method), or `Vector{T}` (list of objects, 2nd method).

The 1st method will throw an error if the input table is empty, and warn if there are more rows
than necessary to construct a single `T`.

The 2nd method will return an empty list for an empty input source, and construct as many `T` as are
found until the input table is exhausted.

`Strapping.construct` utilizes the [StructTypes.jl](https://github.com/JuliaData/StructTypes.jl) package
for determining the `StructTypes.StructType` trait of `T` and constructing an instance appropriately:
  * `StructTypes.Struct`/`StructTypes.Mutable`: field reflection will be used to retrieve values from the input table row, with field customizations respected, like excluded fields, field-specific keyword args, etc.
  * `StructTypes.DictType`: each column name/value of the table row will be used as a key/value pair to be passed to the `DictType` constructor
  * `StructTypes.ArrayType`: column values will be "collected" as an array to be passed to the `ArrayType` constructor
  * `StructTypes.StringType`/`StructTypes.NumberType`/`StructTypes.BoolType`/`StructTypes.NullType`: only the first value of the row will be passed to the scalar type constructor

Note that for `StructTypes.DictType` and `StructTypes.ArrayType`, "aggregate" value/eltypes are not allowed, since
the entire row is treated as key/value pairs or array elements. That means, for example, I can't have a table with rows like `tbl = [(a=1, b=2)]` and try to do `Strapping.construct(Dict{Symbol, Dict{Int, Int}}, tbl)`. It first attempts to map column names to the outer `Dict` keys, (`a` and `b`), but then tries to map the values `1` and `2`
to `Dict{Int, Int}` and fails.

For structs with `ArrayType` fields, the first row values will be used for other scalar fields, and subsequent rows
will be iterated for the `ArrayType` field values. For example, I may wish to construct a type like:
```julia
struct TestResult
    id::Int
    values::Vector{Float64}
end
StructTypes.StructType(::Type{TestResult}) = StructTypes.Struct()
StructTypes.idproperty(::Type{TestResult}) = :id
```

and my input table would look something like, `tbl = (id=[1, 1, 1], values=[3.14, 3.15, 3.16])`. I can then construct my type like:

```julia
julia> Strapping.construct(TestResult, tbl)
TestResult(1, [3.14, 3.15, 3.16])
```

Note that along with defining the `StructTypes.StructType` trait for `TestResult`, I also needed to define `StructTypes.idproperty` to signal which field of my struct is a "unique key" identifier. This enables Strapping to distinguish which rows belong to a particular instance of `TestResult`. This allows the slightly more complicated example of returning multiple `TestResult`s from a single table:

```julia
julia> tbl = (id=[1, 1, 1, 2, 2, 2], values=[3.14, 3.15, 3.16, 40.1, 0.01, 2.34])
(id = [1, 1, 1, 2, 2, 2], values = [3.14, 3.15, 3.16, 40.1, 0.01, 2.34])

julia> Strapping.construct(Vector{TestResult}, tbl)
2-element Array{TestResult,1}:
 TestResult(1, [3.14, 3.15, 3.16])
 TestResult(2, [40.1, 0.01, 2.34])
```

Here, we actually have _two_ `TestResult` objects in our `tbl`, and Strapping uses the `id` field to identify object owners for a row. Note that currently the table rows need to be sorted on the `idproperty` field, i.e. rows belonging to the same object must appear consecutively in the input table rows.

Now let's discuss "aggregate" type fields. Let's say I have a struct like:

```julia
struct Experiment
    id::Int
    name::String
    testresults::TestResult
end
StructTypes.StructType(::Type{Experiment}) = StructTypes.Struct()
StructTypes.idproperty(::Type{Experiment}) = :id
```

So my `Experiment` type also as an `id` field, in addition to a `name` field, and an "aggregate" field of `testresults`. How should the input table source account for `testresults`, which is itself a struct made up of its own `id` and `values` fields? The key here is "flattening" nested structs into a single set of table column names, and utilizing the `StructTypes.fieldprefix` function, which allows specifying a `Symbol` prefix to identify an aggregate field's columns in the table row. So, in the case of our `Experiment`, we can do:

```julia
StructTypes.fieldprefix(::Type{Experiment}, nm::Symbol) = nm == :testresults ? :testresults_ : :_
```

Note that this is the default definition, so we don't really need to define this, but for illustration purposes, we'll walk through it. We're saying that for the `:testresults` field name, we should expect its column names in the table row to start with `:testresults_`. So the table data for an `Experiment` instance, would look something like:

```julia
tbl = (id=[1, 1, 1], name=["exp1", "exp1", "exp1"], testresults_id=[1, 1, 1], testresults_values=[3.14, 3.15, 3.16])
```

This pattern generalizes to structs with multiple aggregate fields, or aggregate fields that themselves have aggregate fields (nested aggregates); in the nested case, the prefixes are concatenated, like `testresults_specifictestresult_id`.
"""
function construct end

function construct(::Type{T}, source; silencewarnings::Bool=false, kw...) where {T}
    rows = Tables.rows(source)
    state = iterate(rows)
    state === nothing && throw(Error("can't construct `$T` from empty source"))
    row, st = state
    x, state = construct(T, rows, row, st; kw...)
    state === nothing || (silencewarnings && println("warning: additional source rows left after reading `$T`"))
    return x
end

function construct(::Type{Vector{T}}, source; kw...) where {T}
    rows = Tables.rows(source)
    state = iterate(rows)
    A = Vector{T}(undef, 0)
    while state !== nothing
        row, st = state
        x, state = construct(T, rows, row, st; kw...)
        push!(A, x)
    end
    return A
end

# collection handler: will iterate results for collection fields
function construct(::Type{T}, rows, row, st; kw...) where {T}
    x = construct(StructTypes.StructType(T), T, row; kw...)
    idprop = StructTypes.idproperty(T)
    # @print 2 idprop
    if idprop !== :_
        id = Tables.getcolumn(row, idprop)
        state = iterate(rows, st)
        if state !== nothing
            while state !== nothing
                row, st = state
                Tables.getcolumn(row, idprop) == id || break
                construct!(StructTypes.StructType(T), x, row; kw...)
                state = iterate(rows, st)
            end
        end
    else
        state = iterate(rows, st)
    end
    return x, state
end

# aggregate handlers
# construct versions construct initial object
# construct! versions take existing object and append additional elements to collection fields
function construct(::StructTypes.CustomStruct, ::Type{T}, row, prefix=Symbol(), offset=Ref(0); kw...) where {T}
    S = StructTypes.lowertype(T)
    x = construct(StructTypes.StructType(S), S, row, prefix, offset; kw...)
    return StructTypes.construct(T, x)
end

function construct(::StructTypes.Struct, ::Type{T}, row, prefix=Symbol(), offset=Ref(0); kw...) where {T}
    # @print 3 (T, prefix, offset)
    return StructTypes.construct(T) do i, nm, TT
        x = construct(StructTypes.StructType(TT), TT, row, Symbol(prefix, StructTypes.fieldprefix(T, nm)), offset, offset[] + 1, Symbol(prefix, nm); kw...)
        # @print 3 x
        return x
    end
end

function construct(::ST, ::Type{T}, row, prefix, offset, i, nm; kw...) where {ST <: Union{StructTypes.Struct, StructTypes.Mutable, StructTypes.CustomStruct}, T}
    # @print 3 (T, prefix, offset, i, nm)
    construct(ST(), T, row, prefix, offset; kw...)
end

function construct!(::StructTypes.CustomStruct, x::T, row, prefix=Symbol(), offset=Ref(0); kw...) where {T}
    y = StructTypes.lower(x)
    return construct!(StructTypes.StructType(y), y, row, prefix, offset; kw...)
end

function construct!(::Union{StructTypes.Struct, StructTypes.Mutable}, x::T, row, prefix=Symbol(), offset=Ref(0); kw...) where {T}
    # @print 3 (T, prefix, offset)
    return StructTypes.foreachfield(x) do i, nm, TT, v
        y = construct!(StructTypes.StructType(TT), v, row, Symbol(prefix, StructTypes.fieldprefix(T, nm)), offset, offset[] + 1, Symbol(prefix, nm); kw...)
        # @print 3 y
        return y
    end
end

function construct!(::ST, x::T, row, prefix, offset, i, nm; kw...) where {ST <: Union{StructTypes.Struct, StructTypes.Mutable, StructTypes.CustomStruct}, T}
    # @print 3 (T, prefix, offset, i, nm)
    construct!(ST(), x, row, prefix, offset; kw...)
end

function construct(::StructTypes.Mutable, ::Type{T}, row, prefix=Symbol(), offset=Ref(0); kw...) where {T}
    # @print 3 (T, prefix, offset)
    x = T()
    StructTypes.mapfields!(x) do i, nm, TT
        y = construct(StructTypes.StructType(TT), TT, row, Symbol(prefix, StructTypes.fieldprefix(T, nm)), offset, offset[] + 1, Symbol(prefix, nm); kw...)
        # @print 3 y
        return y
    end
    return x
end

# fallback for scalars when called from "aggregate" method
construct(ST, T, row, prefix=Symbol(), offset=Ref(0); kw...) = construct(ST, T, row, prefix, offset, 1, Tables.columnnames(row)[1]; kw...)
construct!(ST, x, row, prefix=Symbol(), offset=Ref(0); kw...) = nothing
function construct!(ST, x, row, prefix, offset, i, nm; kw...)
    offset[] += 1
    nothing
end

function construct(::StructTypes.NumberType, ::Type{T}, row, prefix, offset, i, nm; kw...) where {T}
    offset[] += 1
    StructTypes.construct(T, StructTypes.numbertype(T)(getvalue(row, T, i, nm)))
end

# fallback for scalars
function construct(ST, ::Type{T}, row, prefix, offset, i, nm; kw...) where {T}
    offset[] += 1
    StructTypes.construct(T, getvalue(row, T, i, nm))
end

function construct(::StructTypes.Struct, ::Type{Any}, row, prefix, offset, i, nm; kw...)
    offset[] += 1
    Tables.getcolumn(row, nm)
end

# scalar handlers (take a `col` argument)
function getvalue(row, ::Type{T}, col::Int, nm::Symbol) where {T}
    # @print 3 (T, col, nm)
    Tables.getcolumn(row, T, col, nm)
end
function getvalue(row, ::Type{T}, col::Int, nm::Int) where {T}
    # @print 3 (T, col, nm)
    Tables.getcolumn(row, nm)
end

# default aggregate
construct(::StructTypes.Struct, ::Type{Any}, row, prefix=Symbol(), offset=Ref(0); kw...) =
    construct(StructTypes.DictType(), Dict{String, Any}, row; kw...)

# note these don't take `prefix`/`offset` args because we don't expect them to ever be called with a prefix/offset (i.e. only called from line 117)
construct(::StructTypes.DictType, ::Type{T}, row; kw...) where {T} = construct(StructTypes.DictType(), T, row, Symbol, Any; kw...)
construct(::StructTypes.Struct, ::Type{T}, row; kw...) where {T <: NamedTuple} = construct(StructTypes.DictType(), T, row, Symbol, Any; kw...)
construct(::StructTypes.DictType, ::Type{T}, row; kw...) where {T <: NamedTuple} = construct(StructTypes.DictType(), T, row, Symbol, Any; kw...)
construct(::StructTypes.DictType, ::Type{Dict}, row; kw...) = construct(StructTypes.DictType(), Dict, row, String, Any; kw...)
construct(::StructTypes.DictType, ::Type{T}, row; kw...) where {T <: AbstractDict} = construct(StructTypes.DictType(), T, row, keytype(T), valtype(T); kw...)

const AggregateStructTypes = Union{StructTypes.Struct, StructTypes.Mutable, StructTypes.DictType, StructTypes.ArrayType}
@noinline aggregatefieldserror(T) = error("$(StructTypes.StructType(T)) $T not allowed as aggregate field when calling Strapping.construct")

function construct(::StructTypes.DictType, ::Type{T}, row, ::Type{K}, ::Type{V}; kw...) where {T, K, V}
    # check for aggregate types as V
    # the problem is we're already treating the DictType as the aggregate, so
    # it "slurps" all columns as Dict keys, making it
    # impossible to distinguish between multiple aggregate type V
    (V === Any || !(StructTypes.StructType(V) isa AggregateStructTypes)) || aggregatefieldserror(T)
    x = Dict{K, V}()
    for (i, nm) in enumerate(Tables.columnnames(row))
        val = construct(StructTypes.StructType(V), V, row, Symbol(), Ref(0), i, nm; kw...)
        if K === Symbol
            x[nm] = val
        else
            x[StructTypes.construct(K, String(nm))] = val
        end
    end
    return StructTypes.construct(T, x; kw...)
end

# construct(::StructTypes.Struct, ::Type{T}, row, prefix, offset, i, nm; kw...) where {T <: NamedTuple} = construct(StructTypes.DictType(), T, row, prefix, offset, i, nm; kw...)
construct(::StructTypes.DictType, ::Type{T}, row, prefix, offset, i, nm; kw...) where {T} = construct(StructTypes.DictType(), T, row, prefix, offset, i, nm, Symbol, Any; kw...)
construct(::StructTypes.Struct, ::Type{T}, row, prefix, offset, i, nm; kw...) where {T <: NamedTuple} = construct(StructTypes.DictType(), T, row, prefix, offset, i, nm, Symbol, Any; kw...)
construct(::StructTypes.DictType, ::Type{T}, row, prefix, offset, i, nm; kw...) where {T <: NamedTuple} = construct(StructTypes.DictType(), T, row, prefix, offset, i, nm, Symbol, Any; kw...)
construct(::StructTypes.DictType, ::Type{Dict}, row, prefix, offset, i, nm; kw...) = construct(StructTypes.DictType(), Dict, row, prefix, offset, i, nm, String, Any; kw...)
construct(::StructTypes.DictType, ::Type{T}, row, prefix, offset, i, nm; kw...) where {T <: AbstractDict} = construct(StructTypes.DictType(), T, row, prefix, offset, i, nm, keytype(T), valtype(T); kw...)

function construct(::StructTypes.DictType, ::Type{T}, row, prefix, offset, i, nm, ::Type{K}, ::Type{V}; kw...) where {T, K, V}
    (V === Any || !(StructTypes.StructType(V) isa AggregateStructTypes)) || aggregatefieldserror(T)
    x = Dict{K, V}()
    nms = Tables.columnnames(row)
    for j = i:length(nms)
        nm = nms[j]
        val = construct(StructTypes.StructType(V), V, row, prefix, offset, j, nm; kw...)
        if K === Symbol
            x[nm] = val
        else
            x[StructTypes.construct(K, String(nm))] = val
        end
    end
    return StructTypes.construct(T, x; kw...)
end

function construct!(::StructTypes.DictType, x::T, row; kw...) where {T}
    for (i, nm) in enumerate(Tables.columnnames(row))
        v = x[nm]
        V = typeof(v)
        construct!(StructTypes.StructType(V), v, row, Symbol(), Ref(0), i, nm; kw...)
    end
    return
end

construct(::StructTypes.ArrayType, ::Type{T}, row; kw...) where {T} =
    _construct(StructTypes.ArrayType(), T, row, Base.IteratorEltype(T) == Base.HasEltype() ? eltype(T) : Any; kw...)
construct(::StructTypes.ArrayType, ::Type{Tuple}, row; kw...) =
    _construct(StructTypes.ArrayType(), Tuple, row, Any; kw...)

function _construct(::StructTypes.ArrayType, ::Type{T}, row, ::Type{eT}; kw...) where {T, eT}
    # disallow aggregate eT, same problem as DictType; we're treating the ArrayType as our aggregate
    (eT === Any || !(StructTypes.StructType(eT) isa AggregateStructTypes)) || aggregatefieldserror(T)
    nms = Tables.columnnames(row)
    N = length(nms)
    x = Vector{eT}(undef, N)
    for (i, nm) in enumerate(nms)
        x[i] = construct(StructTypes.StructType(eT), eT, row, Symbol(), Ref(0), i, nm; kw...)
    end
    return StructTypes.construct(T, x; kw...)
end

function construct!(::StructTypes.ArrayType, x::T, row; kw...) where {T}
    for (i, nm) in enumerate(Tables.columnnames(row))
        v = x[i]
        V = typeof(v)
        construct!(StructTypes.StructType(V), v, row, Symbol(), Ref(0), i, nm; kw...)
    end
    return
end

function construct(::StructTypes.ArrayType, ::Type{T}, row, ::Type{eT}; kw...) where {T <: Tuple, eT}
    return StructTypes.construct(T) do i, nm, TT
        construct(StructTypes.StructType(TT), TT, row, Symbol(), Ref(0); kw...)
    end
end

function construct(::StructTypes.ArrayType, ::Type{T}, row, prefix, offset, i, nm; kw...) where {T}
    eT = Base.IteratorEltype(T) == Base.HasEltype() ? eltype(T) : Any
    return StructTypes.construct(T, [construct(StructTypes.StructType(eT), eT, row, prefix, offset, i, nm; kw...)]; kw...)
end

function construct!(::StructTypes.ArrayType, x::T, row, prefix, offset, i, nm; kw...) where {T}
    eT = Base.IteratorEltype(T) == Base.HasEltype() ? eltype(T) : Any
    push!(x, construct(StructTypes.StructType(eT), eT, row, prefix, offset, i, nm; kw...))
    return
end

# deconstruct

"""
    Strapping.deconstruct(x::T)
    Strapping.deconstruct(x::Vector{T})

The inverse of `Strapping.construct`, where an object instance `x::T` or `Vector` of objects `x::Vector{T}` is "deconstructed" into a Tables.jl-compatible row iterator. This works following the same patterns outlined in `Strapping.construct` with regards to `ArrayType` and aggregate fields. Specifically, `ArrayType` fields will cause multiple rows to be outputted, one row per collection element, with other scalar fields being repeated in each row. Note this requires no more than one collection field, since otherwise it's ambiguous on how to repeat rows. Similarly for aggregate fields, the field prefix will be used (`StructTypes.fieldprefix`) and nested aggregates will all be flattened into a single list of column names with aggregate prefixes.

In general, this allows outputting any "object" as a 2D table structure that could be stored in any Tables.jl-compatible sink format, e.g. csv file, sqlite table, mysql database table, arrow file, etc.

Note that currently the code requires collection fields to be indexable (i.e. supports `getindex` and `length`) as well as
to have `eltype` defined. Also note that if any collection fields are empty, it will result in the output columns having `Union{T, Missing}`
types since the non-collection fields will still be output as a single row, with collection field(s) being `missing`.
"""
function deconstruct end

mutable struct FieldNode
    index::Int
    name::Symbol
    subfield::Union{FieldNode, Nothing}
end

Base.copy(x::FieldNode) = FieldNode(x.index, x.name, x.subfield === nothing ? nothing : copy(x.subfield))

# recurse until find empty subfield and set `sub`
function setsubfield!(fn::FieldNode, sub::FieldNode)
    if fn.subfield === nothing
        fn.subfield = sub
    else
        setsubfield!(fn.subfield, sub)
    end
    return
end

getfieldvalue(x::T, ind, fn) where {T} = getfieldvalue(StructTypes.StructType(T), x, ind, fn)

function getfieldvalue(::StructTypes.DictType, x, ind, fn)
    val = x[fn.name]
    return getfieldvalue(val, ind, fn.subfield)
end

getfieldvalue(::StructTypes.ArrayType, x, ind, fn) = isempty(x) ? missing : getfieldvalue(x[ind], 0, fn)

function getfieldvalue(::Union{StructTypes.Struct, StructTypes.Mutable}, x, ind, fn)
    val = Core.getfield(x, fn.index)
    return getfieldvalue(val, ind, fn.subfield)
end

function getfieldvalue(::StructTypes.CustomStruct, x, ind, fn)
    y = StructTypes.lower(x)
    return getfieldvalue(StructTypes.StructType(y), y, ind, fn)
end

function getfieldvalue(ST, x, ind, fn)
    @assert fn === nothing
    return x
end

# single object or vector => Tables.rows iterator
struct DeconstructedRowsIterator{T, A}
    values::A
    lens::Vector{Int}
    len::Int
    names::Vector{Symbol}
    types::Vector{Type}
    fields::Vector{FieldNode}
    lookup::Dict{Symbol, Int}
end

Base.eltype(x::DeconstructedRowsIterator{T}) where {T} = DeconstructedRow{T}
Base.length(x::DeconstructedRowsIterator) = x.len

# single object w/ plain scalar fields: 1 obj => 1 row
# single object w/ scalar field and collection field: 1 obj => N rows, where N = length(collection)
# single object w/ aggregate object field: 1 obj => 1 row
# single object w/ aggregate object collection field: 1 obj => N rows, where N = length(collection)
@inline function Base.iterate(x::DeconstructedRowsIterator, (i, j, k)=(1, 1, 1))
    k > x.len && return nothing
    if j == x.lens[i]
        return DeconstructedRow(x.values[i], j, x.names, x.fields, x.lookup), (i + 1, 1, k + 1)
    else
        return DeconstructedRow(x.values[i], j, x.names, x.fields, x.lookup), (i, j + 1, k + 1)
    end
end

Tables.isrowtable(::Type{<:DeconstructedRowsIterator}) = true
Tables.schema(x::DeconstructedRowsIterator) = Tables.Schema(x.names, x.types)

struct DeconstructedRow{T} <: Tables.AbstractRow
    x::T # a single object we're deconstructing
    index::Int # index of this specific row (may be > 1 for objects w/ collection fields)
    names::Vector{Symbol}
    fields::Vector{FieldNode}
    lookup::Dict{Symbol, Int}
end

obj(x::DeconstructedRow) = getfield(x, :x)
ind(x::DeconstructedRow) = getfield(x, :index)
names(x::DeconstructedRow) = getfield(x, :names)
fields(x::DeconstructedRow) = getfield(x, :fields)
lookup(x::DeconstructedRow) = getfield(x, :lookup)

Tables.columnnames(row::DeconstructedRow) = names(row)
Tables.getcolumn(row::DeconstructedRow, ::Type{T}, i::Int, nm::Symbol) where {T} =
    getfieldvalue(obj(row), ind(row), fields(row)[i])
Tables.getcolumn(row::DeconstructedRow, i::Int) =
    getfieldvalue(obj(row), ind(row), fields(row)[i])
Tables.getcolumn(row::DeconstructedRow, nm::Symbol) = Tables.getcolumn(row, lookup(row)[nm])

mutable struct DeconstructClosure{}
    len::Int
    names::Vector{Symbol}
    types::Vector{Type}
    fields::Vector{FieldNode}
    prefix::Symbol
    fieldnode::Union{FieldNode, Nothing}
    i::Int
    processedArrayType::Bool
    emptyArrayType::Bool
    parentType::Type
end

DeconstructClosure() = DeconstructClosure(1, Symbol[], Type[], FieldNode[], Symbol(), nothing, 1, false, false, Any)

# if we're direct root child, sub is full field ancestry (i.e. name/index are name/index of root object)
# if we're further nested than direct root child,
# we make a copy of the fieldnode ancestry and set
# our new node as new leaf node
function getfieldnode(f::DeconstructClosure, sub::FieldNode)
    if f.fieldnode === nothing
        return sub
    else
        fn = copy(f.fieldnode)
        setsubfield!(fn, sub)
        return fn
    end
end

function reset!(f::DeconstructClosure)
    f.i = 1
    f.len = 1
    f.processedArrayType = false
    f.emptyArrayType = false
    return
end

(f::DeconstructClosure)(x::T) where {T} = f(StructTypes.StructType(T), x)

struct EmptyArrayTypeValue end

function (f::DeconstructClosure)(::Union{StructTypes.Struct, StructTypes.Mutable}, x::T) where {T}
    # x is root object
    reset!(f)
    f.parentType = T
    StructTypes.foreachfield(x) do i, nm, TT, v
        # each root field is a separate prefix/fieldnode ancestry branch
        f.prefix = Symbol()
        f.fieldnode = nothing
        f(i, nm, TT, v)
    end
    return
end

function (f::DeconstructClosure)(::StructTypes.CustomStruct, x::T) where {T}
    y = StructTypes.lower(x)
    return f(StructTypes.StructType(y), y)
end

function (f::DeconstructClosure)(::StructTypes.DictType, x::T) where {T}
    # x is root object
    reset!(f)
    f.parentType = T
    i = 1
    for (k, v) in StructTypes.keyvaluepairs(x)
        # each root field is a separate prefix/fieldnode ancestry branch
        f.prefix = Symbol()
        f.fieldnode = nothing
        f(i, k, typeof(v), v)
        i += 1
    end
    return
end

(f::DeconstructClosure)(i, nm, TT, v; kw...) = f(StructTypes.StructType(TT), i, nm, TT, v)
(f::DeconstructClosure)(i, nm, TT; kw...) = f(StructTypes.StructType(TT), i, nm, TT, EmptyArrayTypeValue())
function (f::DeconstructClosure)(::Union{StructTypes.Struct, StructTypes.Mutable}, i, nm, TT, v)
    f.prefix = Symbol(f.prefix, StructTypes.fieldprefix(f.parentType, nm))
    f.parentType = TT
    f.fieldnode = getfieldnode(f, FieldNode(i, nm, nothing))
    StructTypes.foreachfield(f, v)
    return
end

function (f::DeconstructClosure)(::Union{StructTypes.Struct, StructTypes.Mutable}, i, nm, TT, ::EmptyArrayTypeValue)
    f.prefix = Symbol(f.prefix, StructTypes.fieldprefix(f.parentType, nm))
    f.parentType = TT
    f.fieldnode = getfieldnode(f, FieldNode(i, nm, nothing))
    StructTypes.foreachfield(f, TT)
    return
end

function (f::DeconstructClosure)(::StructTypes.DictType, i, nm, TT, x)
    f.prefix = Symbol(f.prefix, StructTypes.fieldprefix(f.parentType, nm))
    f.parentType = TT
    f.fieldnode = getfieldnode(f, FieldNode(i, nm, nothing))
    i = 1
    for (k, v) in StructTypes.keyvaluepairs(x)
        f(i, k, typeof(v), v)
        i += 1
    end
    return
end

(f::DeconstructClosure)(::StructTypes.DictType, i, nm, TT, ::EmptyArrayTypeValue) =
    error("unable to detect field types from empty StructTypes.ArrayType object with StructType.DictType eltype (since DictType keys/values are flattened each as fields and should always contain the same keys from element to element)")

@noinline multiplearraytype(TT) = error("multiple StructTypes.ArrayType objects detected during Strapping.deconstruct; only one is allowed for determining repeated deconstructed rows; see StructTypes.excludes(T) = (:nm,) for excluding fields from deconstructing consideration")
@noinline lengthrequired(T) = error("length(::$T) is required for Strapping.deconstruct")
@noinline eltyperequired(T) = erorr("eltype($T) is required for Strapping.deconstruct")

function (f::DeconstructClosure)(::StructTypes.ArrayType, i, nm, TT, v)
    f.processedArrayType && multiplearraytype(TT)
    f.processedArrayType = true
    Base.haslength(v) || lengthrequired(TT)
    f.len = max(1, length(v))
    T = eltype(v)
    st = iterate(v)
    f.emptyArrayType = st === nothing
    f(i, nm, T, st === nothing ? EmptyArrayTypeValue() : st[1])
    return
end

function (f::DeconstructClosure)(ST, i, nm, TT, v)
    if f.i > length(f.names)
        # first time deconstructing obj
        push!(f.names, Symbol(f.prefix, nm))
        push!(f.types, f.emptyArrayType ? Union{Missing, TT} : TT)
        push!(f.fields, getfieldnode(f, FieldNode(i, nm, nothing)))
    else
        # promote field type if needed
        T = f.types[f.i]
        TTT = promote_type(TT, T)
        if !(TTT <: T)
            f.types[f.i] = f.emptyArrayType ? Union{Missing, TTT} : TTT
        end
    end
    f.i += 1
    return
end

deconstruct(x) = deconstruct([x])

function deconstruct(values::A) where {T, A <: AbstractVector{T}}
    c = DeconstructClosure()
    lens = Int[]
    len = 0
    i = 0
    for x in values
        c(x)
        push!(lens, c.len)
        len += c.len
        i += 1
    end
    lookup = Dict(nm => i for (i, nm) in enumerate(c.names))
    return DeconstructedRowsIterator{T, A}(values, lens, len, c.names, c.types, c.fields, lookup)
end

end # module


text:
using DataFrames
include("RVM.jl")
using .RVM

function standard_score(Obs::AbstractMatrix)
    μ = mean(Obs, 1)
    σ = var(Obs, 1)

    f(M) = (M .- μ) ./ σ
end

function identitize(Obs::AbstractMatrix)
    identity
end

function get_data(ifname::AbstractString, target::Symbol)
    data = DataFrames.readtable(ifname)
    split_target(data, target)
end

function split_target(data::DataFrame, target::Symbol)
    Obs  = convert(Matrix{Float64},
                   data[setdiff(names(data), [target])])
    ObsT = convert(Vector, data[target])
    Obs, ObsT
end

SIZE = 100, 1000 , 2300, 4300
Train, TrainT = get_data("../data/training.csv", :class)

for K in (MLKernels.LinearKernel(), MLKernels.RadialBasisKernel(), MLKernels.GaussianKernel())
    @show "******************"
    @show K
    spec  = RVM.RVMSpec(K, 50, standard_score, 1e5, 1e6)
    for S in SIZE
        @show S
        @time model = RVM.fit(spec, Train[1:S, :], TrainT[1:S])
        @show size(model.RV)
        TrainP = RVM.predict(model, Train)
#=
        @show confusmat(2,
                        labelencode(model.labelmap, TrainP),
                        labelencode(model.labelmap, TrainT))
        @show mean(TrainP .== TrainT)

        Test, TestT = get_data("../data/testing.csv", :class)
        TestP = predict(model, Test)
        @show confusmat(2,
                        labelencode(model.labelmap, TestP),
                        labelencode(model.labelmap, TestT))
        @show mean(TestP .== TestT)
=#
    end
end


text:
# Julia wrapper for header: /usr/include/scip/scip_reader.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeReader(scip, name, desc, extension, readercopy, readerfree, readerread, readerwrite, readerdata)
    ccall((:SCIPincludeReader, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Cstring, Cstring, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{SCIP_READERDATA}), scip, name, desc, extension, readercopy, readerfree, readerread, readerwrite, readerdata)
end

function SCIPincludeReaderBasic(scip, readerptr, name, desc, extension, readerdata)
    ccall((:SCIPincludeReaderBasic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_READER}}, Cstring, Cstring, Cstring, Ptr{SCIP_READERDATA}), scip, readerptr, name, desc, extension, readerdata)
end

function SCIPsetReaderCopy(scip, reader, readercopy)
    ccall((:SCIPsetReaderCopy, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_READER}, Ptr{Cvoid}), scip, reader, readercopy)
end

function SCIPsetReaderFree(scip, reader, readerfree)
    ccall((:SCIPsetReaderFree, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_READER}, Ptr{Cvoid}), scip, reader, readerfree)
end

function SCIPsetReaderRead(scip, reader, readerread)
    ccall((:SCIPsetReaderRead, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_READER}, Ptr{Cvoid}), scip, reader, readerread)
end

function SCIPsetReaderWrite(scip, reader, readerwrite)
    ccall((:SCIPsetReaderWrite, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_READER}, Ptr{Cvoid}), scip, reader, readerwrite)
end

function SCIPfindReader(scip, name)
    ccall((:SCIPfindReader, libscip), Ptr{SCIP_READER}, (Ptr{SCIP_}, Cstring), scip, name)
end

function SCIPgetReaders(scip)
    ccall((:SCIPgetReaders, libscip), Ptr{Ptr{SCIP_READER}}, (Ptr{SCIP_},), scip)
end

function SCIPgetNReaders(scip)
    ccall((:SCIPgetNReaders, libscip), Cint, (Ptr{SCIP_},), scip)
end


text:
#__precompile__()

module MaskRCNN

using Knet, ArgCheck, ArgParse, Images, FileIO, MAT
using Pkg; Pkg.build("Knet")

#export
#    MaskRCNN,
#    config,
#    predict,
#    update

export
    smooth_l1_loss,
    clip_to_window,
    clip_boxes,
    resize_image,
    mold_image,
    compose_image_meta,
    parse_image_meta,
    mold_inputs,
    compute_iou,
    compute_overlaps_bbox,
    extract_bboxes,
    box_refinement_graph,
    matconvnet,
    imgdata,
    resnet50,
    resnet101,
    resnet152,
    get_params

include("utils.jl")
include("mrcnn.jl")
include("resnet.jl")
include("rpn.jl")
include("fpn.jl")
include("roi_align.jl")

end # module


text:
module LookUpTable

export SubstatHasKey,SubstatLookUp,MainstatHasKey,MainstatLookUp

#=================user configured section===================#

# Please add your output substat Sequence that you wish
substatProject=Dict(
    "hp"=>1,
    "hp_percent"=>2,
    "attack"=>3,
    "attack_percent"=>4,
    "defence"=>5,
    "defence_percent"=>6,
    "charge_efficiency"=>7,
    "element_mastery"=>8,
    "critical"=>9,
    "critical_hurt"=>10,
)

substatProjectRev=

# Please add your output mainstat Sequence that you wish
mainstatProject=Dict(
    "hp"=>1,
    "hp_percent"=>2,
    "attack"=>3,
    "attack_percent"=>4,
    "defence"=>5,
    "defence_percent"=>6,
    "critical"=>7,
    "critical_hurt"=>8,
    "charge_efficiency"=>9,
    "heal_add"=>10,
    "element_mastery"=>11,
    "fire_add_hurt"=>12,
    "elec_add_hurt"=>13,
    "water_add_hurt"=>14,
    "wind_add_hurt"=>15,
    "rock_add_hurt"=>16,
    "ice_add_hurt"=>17,
    "physical_add_hurt"=>18,
)

mainstatProjectRev=

# Right side is your variable name for those stats. Also used in comment.
statName=Dict(
    "hp"=>"HP",
    "hp_percent"=>"HP%",
    "attack"=>"ATK",
    "attack_percent"=>"ATK%",
    "defence"=>"DEF",
    "defence_percent"=>"DEF%",
    "charge_efficiency"=>"ELEMENT_CHARGE",
    "element_mastery"=>"ELEMENT_MASTERY",
    "critical"=>"CRIT_RATE",
    "critical_hurt"=>"CRIT_DMG",
    "heal_add"=>10,
    "fire_add_hurt"=>12,
    "elec_add_hurt"=>13,
    "water_add_hurt"=>14,
    "wind_add_hurt"=>15,
    "rock_add_hurt"=>16,
    "ice_add_hurt"=>17,
    "physical_add_hurt"=>18,
)

#==============end of user configured section===============#

# constant look-up table

statIntern=Dict(
    "FIGHT_PROP_HP"=>"hp",
    "FIGHT_PROP_HP_PERCENT"=>"hp_percent",
    "FIGHT_PROP_ATTACK"=>"attack",
    "FIGHT_PROP_ATTACK_PERCENT"=>"attack_percent",
    "FIGHT_PROP_DEFENSE"=>"defence",
    "FIGHT_PROP_DEFENSE_PERCENT"=>"defence_percent",
    "FIGHT_PROP_CRITICAL"=>"critical",
    "FIGHT_PROP_CRITICAL_HURT"=>"critical_hurt",
    "FIGHT_PROP_CHARGE_EFFICIENCY"=>"charge_efficiency",
    "FIGHT_PROP_HEAL_ADD"=>"heal_add",
    "FIGHT_PROP_ELEMENT_MASTERY"=>"element_mastery",
    "FIGHT_PROP_FIRE_ADD_HURT"=>"fire_add_hurt",
    "FIGHT_PROP_ELEC_ADD_HURT"=>"elec_add_hurt",
    "FIGHT_PROP_WATER_ADD_HURT"=>"water_add_hurt",
    "FIGHT_PROP_WIND_ADD_HURT"=>"wind_add_hurt",
    "FIGHT_PROP_ROCK_ADD_HURT"=>"rock_add_hurt",
    "FIGHT_PROP_GRASS_ADD_HURT"=>"grass_add_hurt",
    "FIGHT_PROP_ICE_ADD_HURT"=>"ice_add_hurt",
    "FIGHT_PROP_PHYSICAL_ADD_HURT"=>"physical_add_hurt",
    "FIGHT_PROP_FIRE_SUB_HURT"=>"fire_sub_hurt",
);

function SubstatHasKey(s::String)
    return haskey(substatProject,statIntern[s])
end

function SubstatLookUp(s::String)
    return substatProject[statIntern[s]]
end

function MainstatHasKey(s::String)
    return haskey(mainstatProject,statIntern[s])
end

function MainstatLookUp(s::String)
    return mainstatProject[statIntern[s]]
end

end

max_stars_repo_name: hypergraphUniverse/GenshinDataConverter

max_issues_repo_name: hypergraphUniverse/GenshinDataConverter

max_forks_repo_name: hypergraphUniverse/GenshinDataConverter


text:
# MDP parameters, ValueIteration minimizes the cost => cost is positive, reward is negative
no_states = 10
_horizon = 5
actions = [:l, :r]
actionCost = 1.
actionsImpact = Base.ImmutableDict(:l => -1, :r => 1)
reward_states = [1, no_states]
reward = -10.
discount_factor = 1.
noise = .6

mdp = CustomFHExample(no_states, _horizon, actions, actionCost, actionsImpact, reward_states, reward, discount_factor, noise)

println("Verbose Check")
fhsolver = FiniteHorizonSolver(verbose=true)
FHPolicy = FiniteHorizonPOMDPs.solve(fhsolver, mdp);


text:
#
# The mul! methods that contract derivative operators with arrays.
#
# These methods handle the {N} dimension type parameter.
#
# There are four methods.  One is a general method for multidimensional
# arrays, and the others are short cuts for 1, 2 and 3 dimensions.  The 1 dimensional
# method for AbstractVector is defined in convolutions.jl, and the
# others are defined here.
#
# At the interior points of an evenly spaced grid, a finite difference
# operator is a convolution.  Where possible, the heavy lifting is
# done by the efficient conv!  routine from NNlib.  The routines
# defined in convolutions.jl cover the cases where that isn't possible.
#

# Fallback mul! implementation for a single DerivativeOperator operating on an AbstractArray
function LinearAlgebra.mul!(x_temp::AbstractArray{T}, A::DerivativeOperator{T,N}, M::AbstractArray{T}; overwrite = true) where {T,N}

    # Check that x_temp has valid dimensions, allowing unnecessary padding in M
    v = zeros(ndims(x_temp))
    v .= 2
    @assert all(([size(x_temp)...] .== [size(M)...])
        .| (([size(x_temp)...] .+ v) .== [size(M)...])
        )

    # Check that axis of differentiation is in the dimensions of M and x_temp
    ndims_M = ndims(M)
    @assert N <= ndims_M
    @assert size(x_temp, N) + 2 == size(M, N) # differentiated dimension must be padded

    alldims = [1:ndims(M);]
    otherdims = setdiff(alldims, N)

    idx = Any[first(ind) for ind in axes(M)]
    nidx = length(otherdims)

    dims_M = [axes(M)...]
    dims_x_temp = [axes(x_temp)...]
    minimal_padding_indices = map(enumerate(dims_x_temp)) do (dim, val)
        if dim == N || length(dims_x_temp[dim]) == length(dims_M[dim])
            Colon()
        else
            dims_M[dim][begin+1:end-1]
        end
    end
    minimally_padded_M = view(M, minimal_padding_indices...)

    itershape = tuple(dims_x_temp[otherdims]...)
    indices = Iterators.drop(CartesianIndices(itershape), 0)

    setindex!(idx, :, N)
    for I in indices
        # replace all elements of idx with corresponding elt of I, except at index N
        Base.replace_tuples!(nidx, idx, idx, otherdims, I)
        mul!(view(x_temp, idx...), A, view(minimally_padded_M, idx...), overwrite = overwrite)
    end
end

# A more efficient mul! implementation for a single, regular-grid, centered difference,
# scalar coefficient, non-winding, DerivativeOperator operating on a 2D or 3D AbstractArray
for MT in [2,3]
    @eval begin
        function LinearAlgebra.mul!(x_temp::AbstractArray{T,$MT}, A::DerivativeOperator{T,N,false,T2,S1,S2,T3}, M::AbstractArray{T,$MT}) where
                                                                            {T,N,T2,SL,S1<:SArray{Tuple{SL},T,1,SL},S2,T3<:Union{Nothing,Number}}
            # Check that x_temp has valid dimensions, allowing unnecessary padding in M
            v = zeros(ndims(x_temp))
            v .= 2
            @assert all(([size(x_temp)...] .== [size(M)...])
                .| (([size(x_temp)...] .+ v) .== [size(M)...])
                )

            # Check that axis of differentiation is in the dimensions of M and x_temp
            ndims_x_temp = ndims(x_temp)
            @assert N <= ndims_x_temp
            @assert size(x_temp, N) + 2 == size(M, N) # differentiated dimension must be padded

            # Determine padding for NNlib.conv!
            bpc = A.boundary_point_count
            pad = zeros(Int64, ndims_x_temp)
            pad[N] = bpc

            # Reshape x_temp for NNlib.conv!
            _x_temp = reshape(x_temp, (size(x_temp)...,1,1))

            # Reshape M for NNlib.conv!
            dims_M = [axes(M)...]
            dims_x_temp = [axes(x_temp)...]
            minimal_padding_indices = map(enumerate(dims_x_temp)) do (dim, val)
                if dim == N || length(dims_x_temp[dim]) == length(dims_M[dim])
                    Colon()
                else
                    dims_M[dim][begin+1:end-1]
                end
            end
            minimally_padded_M = view(M, minimal_padding_indices...)
            _M = reshape(minimally_padded_M, (size(minimally_padded_M)...,1,1))

            # Setup W, the kernel for NNlib.conv!
            s = A.stencil_coefs
            sl = A.stencil_length
            Wdims = ones(Int64, ndims(_x_temp))
            Wdims[N] = sl
            W = zeros(Wdims...)
            Widx = Any[Wdims...]
            setindex!(Widx,:,N)
            coeff = A.coefficients === nothing ? true : A.coefficients
            W[Widx...] = coeff*s

            cv = DenseConvDims(_M, W, padding=pad, flipkernel=true)
            conv!(_x_temp, _M, W, cv)

            # Now deal with boundaries
            if bpc > 0
                dims_minimal_M = [axes(minimally_padded_M)...]
                alldims = [1:ndims(minimally_padded_M);]
                otherdims = setdiff(alldims, N)

                idx = Any[first(ind) for ind in axes(minimally_padded_M)]
                itershape = tuple(dims_minimal_M[otherdims]...)
                nidx = length(otherdims)
                indices = Iterators.drop(CartesianIndices(itershape), 0)

                setindex!(idx, :, N)
                for I in indices
                    Base.replace_tuples!(nidx, idx, idx, otherdims, I)
                    convolve_BC_left!(view(x_temp, idx...), view(minimally_padded_M, idx...), A)
                    convolve_BC_right!(view(x_temp, idx...), view(minimally_padded_M, idx...), A)
                end
            end
        end
    end
end

###########################################

function *(A::DerivativeOperator{T,N},M::AbstractArray{T}) where {T<:Real,N}
    size_x_temp = [size(M)...]
    size_x_temp[N] -= 2
    x_temp = zeros(promote_type(eltype(A),eltype(M)), size_x_temp...)
    LinearAlgebra.mul!(x_temp, A, M)
    return x_temp
end

function *(c::Number, A::DerivativeOperator{T,N,Wind}) where {T,N,Wind}
    coefficients = A.coefficients === nothing ? oneunit(T) .* c : c .* A.coefficients
    DerivativeOperator{T,N,Wind,typeof(A.dx),typeof(A.stencil_coefs),
                       typeof(A.low_boundary_coefs),typeof(coefficients),
                       typeof(A.coeff_func)}(
        A.derivative_order, A.approximation_order,
        A.dx, A.len, A.stencil_length,
        A.stencil_coefs,
        A.boundary_stencil_length,
        A.boundary_point_count,
        A.low_boundary_coefs,
        A.high_boundary_coefs,coefficients,A.coeff_func)
end

function *(c::AbstractVector{<:Number}, A::DerivativeOperator{T,N,Wind}) where {T,N,Wind}
    if length(c) != A.len
        throw(DimensionMismatch("length of c ($(length(c))) must match length of A ($A.len)"))
    end
    coefficients = A.coefficients === nothing ? c : c .* A.coefficients
    DerivativeOperator{T,N,Wind,typeof(A.dx),typeof(A.stencil_coefs),
                       typeof(A.low_boundary_coefs),typeof(coefficients),
                       typeof(A.coeff_func)}(
        A.derivative_order, A.approximation_order,
        A.dx, A.len, A.stencil_length,
        A.stencil_coefs,
        A.boundary_stencil_length,
        A.boundary_point_count,
        A.low_boundary_coefs,
        A.high_boundary_coefs,coefficients,A.coeff_func)
end

###########################################

# A more efficient mul! implementation for compositions of operators which may include regular-grid, centered difference,
# scalar coefficient, non-winding, DerivativeOperator, operating on a 2D or 3D AbstractArray
function LinearAlgebra.mul!(x_temp::AbstractArray{T,2}, A::AbstractDiffEqCompositeOperator, M::AbstractArray{T,2}) where {T}

    # opsA operators satisfy conditions for NNlib.conv! call, opsB operators do not
    opsA = DerivativeOperator[]
    opsB = DerivativeOperator[]
    for L in A.ops
        if (L.coefficients isa Number || L.coefficients === nothing) && use_winding(L) == false && L.dx isa Number
            push!(opsA, L)
        else
            push!(opsB,L)
        end
    end

    # Check that we can make at least one NNlib.conv! call
    if !isempty(opsA)
        ndimsM = ndims(M)
        Wdims = ones(Int64,ndimsM)
        pad = zeros(Int64, ndimsM)

        # compute dimensions of interior kernel W
        # Here we still use A.ops since operators in opsB may indicate that
        # we have more padding to account for
        for L in A.ops
            axis = typeof(L).parameters[2]
            @assert axis <= ndimsM
            Wdims[axis] = max(Wdims[axis],L.stencil_length)
            pad[axis] = max(pad[axis], L.boundary_point_count)
        end

        # create zero-valued kernel
        W = zeros(T, Wdims...)
        mid_Wdims = div.(Wdims,2).+1
        idx = div.(Wdims,2).+1

        # add to kernel each stencil
        for L in opsA
            s = L.stencil_coefs
            sl = L.stencil_length
            axis = typeof(L).parameters[2]
            offset = convert(Int64,(Wdims[axis] - sl)/2)
            coeff = L.coefficients isa Number ? L.coefficients : true
            for i in offset+1:Wdims[axis]-offset
                idx[axis]=i
                W[idx...] += coeff*s[i-offset]
                idx[axis] = mid_Wdims[axis]
            end
        end

        # Reshape x_temp for NNlib.conv!
        _x_temp = reshape(x_temp, (size(x_temp)...,1,1))

        # Reshape M for NNlib.conv!
        _M = reshape(M, (size(M)...,1,1))

        _W = reshape(W, (size(W)...,1,1))

        # Call NNlib.conv!
        cv = DenseConvDims(_M, _W, padding=pad, flipkernel=true)
        conv!(_x_temp, _M, _W, cv)

        # convolve boundary and interior points near boundary
        # partition operator indices along axis of differentiation
        if pad[1] > 0 || pad[2] > 0
            ops_1 = Int64[]
            ops_1_max_bpc_idx = [0]
            ops_2 = Int64[]
            ops_2_max_bpc_idx = [0]
            for i in 1:length(opsA)
                L = opsA[i]
                if typeof(L).parameters[2] == 1
                    push!(ops_1,i)
                    if L.boundary_point_count == pad[1]
                        ops_1_max_bpc_idx[1] = i
                    end
                else
                    push!(ops_2,i)
                    if L.boundary_point_count == pad[2]
                        ops_2_max_bpc_idx[1]= i
                    end
                end
            end

            # need offsets since some axis may have ghost nodes and some may not
            offset_x = 0
            offset_y = 0

            if length(ops_2) > 0
                offset_x = 1
            end
            if length(ops_1) > 0
                offset_y = 1
            end

            # convolve boundaries and unaccounted for interior in axis 1
            if length(ops_1) > 0
                for i in 1:size(x_temp)[2]
                    convolve_BC_left!(view(x_temp,:,i), view(M,:,i+offset_x), opsA[ops_1_max_bpc_idx...])
                    convolve_BC_right!(view(x_temp,:,i), view(M,:,i+offset_x), opsA[ops_1_max_bpc_idx...])
                    if i <= pad[2] || i > size(x_temp)[2]-pad[2]
                        convolve_interior!(view(x_temp,:,i), view(M,:,i+offset_x), opsA[ops_1_max_bpc_idx...])
                    end

                    for Lidx in ops_1
                        if Lidx != ops_1_max_bpc_idx[1]
                            convolve_BC_left!(view(x_temp,:,i), view(M,:,i+offset_x), opsA[Lidx], overwrite = false)
                            convolve_BC_right!(view(x_temp,:,i), view(M,:,i+offset_x), opsA[Lidx], overwrite = false)
                            if i <= pad[2] || i > size(x_temp)[2]-pad[2]
                                convolve_interior!(view(x_temp,:,i), view(M,:,i+offset_x), opsA[Lidx], overwrite = false)
                            elseif pad[1] - opsA[Lidx].boundary_point_count > 0
                                convolve_interior!(view(x_temp,:,i), view(M,:,i+offset_x), opsA[Lidx], overwrite = false, add_range = true, offset = pad[1] - opsA[Lidx].boundary_point_count)
                            end
                        end
                    end
                end
            end
            # convolve boundaries and unaccounted for interior in axis 2
            if length(ops_2) > 0
                for i in 1:size(x_temp)[1]
                    # in the case of no axis 1 operators, we need to overwrite x_temp
                    if length(ops_1) == 0
                        convolve_BC_left!(view(x_temp,i,:), view(M,i+offset_y,:), opsA[ops_2_max_bpc_idx...])
                        convolve_BC_right!(view(x_temp,i,:), view(M,i+offset_y,:), opsA[ops_2_max_bpc_idx...])
                        if i <= pad[1] || i > size(x_temp)[1]-pad[1]
                            convolve_interior!(view(x_temp,i,:), view(M,i+offset_y,:), opsA[ops_2_max_bpc_idx...])
                        end

                    else
                        convolve_BC_left!(view(x_temp,i,:), view(M,i+offset_y,:), opsA[ops_2_max_bpc_idx...], overwrite = false)
                        convolve_BC_right!(view(x_temp,i,:), view(M,i+offset_y,:), opsA[ops_2_max_bpc_idx...], overwrite = false)
                        if i <= pad[1] || i > size(x_temp)[1]-pad[1]
                            convolve_interior!(view(x_temp,i,:), view(M,i+offset_y,:), opsA[ops_2_max_bpc_idx...], overwrite = false)
                        end

                    end
                    for Lidx in ops_2
                        if Lidx != ops_2_max_bpc_idx[1]
                            convolve_BC_left!(view(x_temp,i,:), view(M,i+offset_y,:), opsA[Lidx], overwrite = false)
                            convolve_BC_right!(view(x_temp,i,:), view(M,i+offset_y,:), opsA[Lidx], overwrite = false)
                            if i <= pad[1] || i > size(x_temp)[1]-pad[1]
                                convolve_interior!(view(x_temp,i,:), view(M,i+offset_y,:), opsA[Lidx], overwrite = false)
                            elseif pad[2] - opsA[Lidx].boundary_point_count > 0
                                convolve_interior!(view(x_temp,i,:), view(M,i+offset_y,:), opsA[Lidx], overwrite = false, add_range = true, offset = pad[2] - opsA[Lidx].boundary_point_count)
                            end
                        end
                    end
                end
            end
        end

        # Here we compute mul! (additively) for every operator in opsB

        operating_dims = zeros(Int64,2)
        # need to consider all dimensions and operators to determine the truncation
        # of M to x_temp
        for L in A.ops
            if diff_axis(L) == 1
                operating_dims[1] = 1
            else
                operating_dims[2] = 1
            end
        end

        x_temp_1, x_temp_2 = size(x_temp)

        for L in opsB
            N = diff_axis(L)
            if N == 1
                if operating_dims[2] == 1
                    mul!(x_temp,L,view(M,1:x_temp_1+2,1:x_temp_2), overwrite = false)
                else
                    mul!(x_temp,L,M, overwrite = false)
                end
            else
                if operating_dims[1] == 1
                    mul!(x_temp,L,view(M,1:x_temp_1,1:x_temp_2+2), overwrite = false)
                else
                    mul!(x_temp,L,M, overwrite = false)
                end
            end
        end

    # The case where we call everything in A.ops using the fallback mul!
    else
        #operating_dims
        operating_dims = zeros(Int64,2)
        for L in A.ops
            if diff_axis(L) == 1
                operating_dims[1] = 1
            else
                operating_dims[2] = 1
            end
        end

        x_temp_1, x_temp_2 = size(x_temp)

        # Handle first case non-additively
        N = diff_axis(A.ops[1])
        if N == 1
            if operating_dims[2] == 1
                mul!(x_temp,A.ops[1],view(M,1:x_temp_1+2,1:x_temp_2))
            else
                mul!(x_temp,A.ops[1],M)
            end
        else
            if operating_dims[1] == 1
                mul!(x_temp,A.ops[1],view(M,1:x_temp_1,1:x_temp_2+2))
            else
                mul!(x_temp,A.ops[1],M)
            end
        end

        for L in A.ops[2:end]
            N = diff_axis(L)
            if N == 1
                if operating_dims[2] == 1
                    mul!(x_temp,L,view(M,1:x_temp_1+2,1:x_temp_2), overwrite = false)
                else
                    mul!(x_temp,L,M, overwrite = false)
                end
            else
                if operating_dims[1] == 1
                    mul!(x_temp,L,view(M,1:x_temp_1,1:x_temp_2+2), overwrite = false)
                else
                    mul!(x_temp,L,M, overwrite = false)
                end
            end
        end
    end
end

# A more efficient mul! implementation for compositions of operators which may include regular-grid, centered difference,
# scalar coefficient, non-winding, DerivativeOperator, operating on a 2D or 3D AbstractArray
function LinearAlgebra.mul!(x_temp::AbstractArray{T,3}, A::AbstractDiffEqCompositeOperator, M::AbstractArray{T,3}) where {T}

    # opsA operators satisfy conditions for NNlib.conv! call, opsB operators do not
    opsA = DerivativeOperator[]
    opsB = DerivativeOperator[]
    for L in A.ops
        if (L.coefficients isa Number || L.coefficients === nothing) && use_winding(L) == false && L.dx isa Number
            push!(opsA, L)
        else
            push!(opsB,L)
        end
    end

    # Check that we can make at least one NNlib.conv! call
    if !isempty(opsA)
        ndimsM = ndims(M)
        Wdims = ones(Int64,ndimsM)
        pad = zeros(Int64, ndimsM)

        # compute dimensions of interior kernel W
        # Here we still use A.ops since operators in opsB may indicate that
        # we have more padding to account for
        for L in A.ops
            axis = typeof(L).parameters[2]
            @assert axis <= ndimsM
            Wdims[axis] = max(Wdims[axis],L.stencil_length)
            pad[axis] = max(pad[axis], L.boundary_point_count)
        end

        # create zero-valued kernel
        W = zeros(T, Wdims...)
        mid_Wdims = div.(Wdims,2).+1
        idx = div.(Wdims,2).+1

        # add to kernel each stencil
        for L in opsA
            s = L.stencil_coefs
            sl = L.stencil_length
            axis = typeof(L).parameters[2]
            offset = convert(Int64,(Wdims[axis] - sl)/2)
            coeff = L.coefficients isa Number ? L.coefficients : true
            for i in offset+1:Wdims[axis]-offset
                idx[axis]=i
                W[idx...] += coeff*s[i-offset]
                idx[axis] = mid_Wdims[axis]
            end
        end

        # Reshape x_temp for NNlib.conv!
        _x_temp = reshape(x_temp, (size(x_temp)...,1,1))

        # Reshape M for NNlib.conv!
        _M = reshape(M, (size(M)...,1,1))

        _W = reshape(W, (size(W)...,1,1))

        # Call NNlib.conv!
        cv = DenseConvDims(_M, _W, padding=pad, flipkernel=true)
        conv!(_x_temp, _M, _W, cv)

        # convolve boundary and interior points near boundary
        # partition operator indices along axis of differentiation
        if pad[1] > 0 || pad[2] > 0 || pad[3] > 0
            ops_1 = Int64[]
            ops_1_max_bpc_idx = [0]
            ops_2 = Int64[]
            ops_2_max_bpc_idx = [0]
            ops_3 = Int64[]
            ops_3_max_bpc_idx = [0]

            for i in 1:length(opsA)
                L = opsA[i]
                if typeof(L).parameters[2] == 1
                    push!(ops_1,i)
                    if L.boundary_point_count == pad[1]
                        ops_1_max_bpc_idx[1] = i
                    end
                elseif typeof(L).parameters[2] == 2
                    push!(ops_2,i)
                    if L.boundary_point_count == pad[2]
                        ops_2_max_bpc_idx[1]= i
                    end
                else
                    push!(ops_3,i)
                    if L.boundary_point_count == pad[3]
                        ops_3_max_bpc_idx[1]= i
                    end
                end
            end

            # need offsets since some axis may have ghost nodes and some may not
            offset_x = 0
            offset_y = 0
            offset_z = 0

            if length(ops_1) > 0
                offset_x = 1
            end
            if length(ops_2) > 0
                offset_y = 1
            end
            if length(ops_3) > 0
                offset_z = 1
            end

            # convolve boundaries and unaccounted for interior in axis 1
            if length(ops_1) > 0
                for i in 1:size(x_temp)[2]
                    for j in 1:size(x_temp)[3]
                        convolve_BC_left!(view(x_temp,:,i,j), view(M,:,i+offset_y,j+offset_z), opsA[ops_1_max_bpc_idx...])
                        convolve_BC_right!(view(x_temp,:,i,j), view(M,:,i+offset_y,j+offset_z), opsA[ops_1_max_bpc_idx...])
                        if i <= pad[2] || i > size(x_temp)[2]-pad[2] || j <= pad[3] || j > size(x_temp)[3]-pad[3]
                            convolve_interior!(view(x_temp,:,i,j), view(M,:,i+offset_y,j+offset_z), opsA[ops_1_max_bpc_idx...])
                        end

                        for Lidx in ops_1
                            if Lidx != ops_1_max_bpc_idx[1]
                                convolve_BC_left!(view(x_temp,:,i,j), view(M,:,i+offset_y,j+offset_z), opsA[Lidx], overwrite = false)
                                convolve_BC_right!(view(x_temp,:,i,j), view(M,:,i+offset_y,j+offset_z), opsA[Lidx], overwrite = false)
                                if i <= pad[2] || i > size(x_temp)[2]-pad[2] || j <= pad[3] || j > size(x_temp)[3]-pad[3]
                                    convolve_interior!(view(x_temp,:,i,j), view(M,:,i+offset_y,j+offset_z), opsA[Lidx], overwrite = false)
                                elseif pad[1] - opsA[Lidx].boundary_point_count > 0
                                    convolve_interior!(view(x_temp,:,i,j), view(M,:,i+offset_y,j+offset_z), opsA[Lidx], overwrite = false, add_range = true, offset = pad[1] - opsA[Lidx].boundary_point_count)
                                end
                            end
                        end
                    end
                end
            end
            # convolve boundaries and unaccounted for interior in axis 2
            if length(ops_2) > 0
                for i in 1:size(x_temp)[1]
                    for j in 1:size(x_temp)[3]
                        # in the case of no axis 1 operators, we need to overwrite x_temp
                        if length(ops_1) == 0
                            convolve_BC_left!(view(x_temp,i,:,j), view(M,i+offset_x,:,j+offset_z), opsA[ops_2_max_bpc_idx...])
                            convolve_BC_right!(view(x_temp,i,:,j), view(M,i+offset_x,:,j+offset_z), opsA[ops_2_max_bpc_idx...])
                            if i <= pad[1] || i > size(x_temp)[1]-pad[1] || j <= pad[3] || j > size(x_temp)[3]-pad[3]
                                convolve_interior!(view(x_temp,i,:,j), view(M,i+offset_x,:,j+offset_z), opsA[ops_2_max_bpc_idx...])
                            end

                        else
                            convolve_BC_left!(view(x_temp,i,:,j), view(M,i+offset_x,:,j+offset_z), opsA[ops_2_max_bpc_idx...], overwrite = false)
                            convolve_BC_right!(view(x_temp,i,:,j), view(M,i+offset_x,:,j+offset_z), opsA[ops_2_max_bpc_idx...], overwrite = false)
                            if i <= pad[1] || i > size(x_temp)[1]-pad[1] || j <= pad[3] || j > size(x_temp)[3]-pad[3]
                                convolve_interior!(view(x_temp,i,:,j), view(M,i+offset_y,:,j+offset_z), opsA[ops_2_max_bpc_idx...], overwrite = false)
                            end

                        end
                        for Lidx in ops_2
                            if Lidx != ops_2_max_bpc_idx[1]
                                convolve_BC_left!(view(x_temp,i,:,j), view(M,i+offset_x,:,j+offset_z), opsA[Lidx], overwrite = false)
                                convolve_BC_right!(view(x_temp,i,:,j), view(M,i+offset_x,:,j+offset_z), opsA[Lidx], overwrite = false)
                                if i <= pad[1] || i > size(x_temp)[1]-pad[1] || j <= pad[3] || j > size(x_temp)[3]-pad[3]
                                    convolve_interior!(view(x_temp,i,:,j), view(M,i+offset_x,:,j+offset_z), opsA[Lidx], overwrite = false)
                                elseif pad[2] - opsA[Lidx].boundary_point_count > 0
                                    convolve_interior!(view(x_temp,i,:,j), view(M,i+offset_x,:,j+offset_z), opsA[Lidx], overwrite = false, add_range = true, offset = pad[2] - opsA[Lidx].boundary_point_count)
                                end
                            end
                        end
                    end
                end
            end
            # convolve boundaries and unaccounted for interior in axis 3
            if length(ops_3) > 0
                for i in 1:size(x_temp)[1]
                    for j in 1:size(x_temp)[2]
                        # in the case of no axis 1 and 2 operators, we need to overwrite x_temp
                        if length(ops_1) == 0 && length(ops_2) == 0
                            convolve_BC_left!(view(x_temp,i,j,:), view(M,i+offset_x,j+offset_y,:), opsA[ops_3_max_bpc_idx...])
                            convolve_BC_right!(view(x_temp,i,j,:), view(M,i+offset_x,j+offset_y,:), opsA[ops_3_max_bpc_idx...])
                            if i <= pad[1] || i > size(x_temp)[1]-pad[1]
                                convolve_interior!(view(x_temp,i,j,:), view(M,i+offset_x,j+offset_y,:), opsA[ops_3_max_bpc_idx...])
                            end

                        else
                            convolve_BC_left!(view(x_temp,i,j,:), view(M,i+offset_x,j+offset_y,:), opsA[ops_3_max_bpc_idx...], overwrite = false)
                            convolve_BC_right!(view(x_temp,i,j,:), view(M,i+offset_x,j+offset_y,:), opsA[ops_3_max_bpc_idx...], overwrite = false)
                            if i <= pad[1] || i > size(x_temp)[1]-pad[1] || j <= pad[2] || j > size(x_temp)[2]-pad[2]
                                convolve_interior!(view(x_temp,i,j,:), view(M,i+offset_x,j+offset_y,:), opsA[ops_3_max_bpc_idx...], overwrite = false)
                            end

                        end
                        for Lidx in ops_3
                            if Lidx != ops_3_max_bpc_idx[1]
                                convolve_BC_left!(view(x_temp,i,j,:), view(M,i+offset_x,j+offset_y,:), opsA[Lidx], overwrite = false)
                                convolve_BC_right!(view(x_temp,i,j,:), view(M,i+offset_x,j+offset_y,:), opsA[Lidx], overwrite = false)
                                if i <= pad[1] || i > size(x_temp)[1]-pad[1] || j <= pad[2] || j > size(x_temp)[2]-pad[2]
                                    convolve_interior!(view(x_temp,i,j,:), view(M,i+offset_x,j+offset_y,:), opsA[Lidx], overwrite = false)
                                elseif pad[3] - opsA[Lidx].boundary_point_count > 0
                                    convolve_interior!(view(x_temp,i,j,:), view(M,i+offset_x,j+offset_y,:), opsA[Lidx], overwrite = false, add_range = true, offset = pad[3] - opsA[Lidx].boundary_point_count)
                                end
                            end
                        end
                    end
                end
            end
        end

        # Here we compute mul! (additively) for every operator in opsB

        operating_dims = zeros(Int64,3)
        # need to consider all dimensions and operators to determine the truncation
        # of M to x_temp
        for L in A.ops
            operating_dims[diff_axis(L)] = 1
        end

        x_temp_1, x_temp_2, x_temp_3 = size(x_temp)

        for L in opsB
            N = diff_axis(L)
            if N == 1
                mul!(x_temp, L, view(M,1:x_temp_1+2,1:x_temp_2,1:x_temp_3), overwrite = false)
            elseif N == 2
                 mul!(x_temp, L, view(M,1:x_temp_1,1:x_temp_2+2,1:x_temp_3), overwrite = false)
            else
                mul!(x_temp, L, view(M,1:x_temp_1,1:x_temp_2,1:x_temp_3+2), overwrite = false)
            end
        end

    # The case where we call everything in A.ops using the fallback mul!
    else
        # operating_dims indicates which dimensions we are multiplying along
        operating_dims = zeros(Int64,3)
        for L in A.ops
            operating_dims[diff_axis(L)] = 1
        end

        x_temp_1, x_temp_2, x_temp_3 = size(x_temp)

        # Handle first case non-additively
        N = diff_axis(A.ops[1])

        if N == 1
            mul!(x_temp, A.ops[1], view(M,1:x_temp_1+2,1:x_temp_2,1:x_temp_3))
        elseif N == 2
             mul!(x_temp, A.ops[1], view(M,1:x_temp_1,1:x_temp_2+2,1:x_temp_3))
        else
            mul!(x_temp, A.ops[1], view(M,1:x_temp_1,1:x_temp_2,1:x_temp_3+2))
        end

        for L in A.ops[2:end]
            N = diff_axis(L)
            if N == 1
                mul!(x_temp, L, view(M,1:x_temp_1+2,1:x_temp_2,1:x_temp_3), overwrite = false)
            elseif N == 2
                 mul!(x_temp, L, view(M,1:x_temp_1,1:x_temp_2+2,1:x_temp_3), overwrite = false)
            else
                mul!(x_temp, L, view(M,1:x_temp_1,1:x_temp_2,1:x_temp_3+2), overwrite = false)
            end
        end
    end
end


text:
## Extra features not in base. Present for backwards compatibility.

export cartesian_linear, linear, @forcartesian, @indexedvariable, @nlinear, @nlookup, @nrefshift

macro forcartesian(sym, sz, ex)
    idim = gensym()
    N = gensym()
    sz1 = gensym()
    isdone = gensym()
    quote
        if !(isempty($(esc(sz))) || prod($(esc(sz))) == 0)
            $N = length($(esc(sz)))
            $sz1 = $(esc(sz))[1]
            $isdone = false
            $(esc(sym)) = ones(Int, $N)
            while !$isdone
                $(esc(ex))
                if ($(esc(sym))[1] += 1) > $sz1
                    $idim = 1
                    while $(esc(sym))[$idim] > $(esc(sz))[$idim] && $idim < $N
                        $(esc(sym))[$idim] = 1
                        $idim += 1
                        $(esc(sym))[$idim] += 1
                    end
                    $isdone = $(esc(sym))[$N] > $(esc(sz))[$N]
                end
            end
        end
    end
end

function cartesian_linear(A::AbstractArray, I::Vector{Int})
    k = 0
    for j = length(I):-1:1
        k = size(A, j)*k + I[j]-1
    end
    k += 1
end

# Generate expression A[i1+j1, i2+j2, ...]
macro nrefshift(N, A, sym, shiftexpr)
    _nrefshift(N, A, sym, shiftexpr)
end

# ... using an offset symbol. This is useful with nested @nloops
function _nrefshift(N::Int, A::Symbol, iter1::Symbol, iter2::Symbol)
    vars = [ :($(namedvar(iter1, i))+$(namedvar(iter2, i))) for i = 1:N ]
    Expr(:escape, Expr(:ref, A, vars...))
end

# ... using a shiftexpr, e.g., A[i1, i2+1, ...] with d->(d==2)?1:0
function _nrefshift(N::Int, A::Symbol, sym::Symbol, shiftexpr::Expr)
    vars = [ exprresolve(:($(namedvar(sym, i))+$(inlineanonymous(shiftexpr, i)))) for i = 1:N ]
    Expr(:escape, Expr(:ref, A, vars...))
end

# Generate expression A[ I1[i1], I2[i2], ... ]
macro nlookup(N, A, indexes, itersym)
    _nlookup(N, A, indexes, itersym)
end

function _nlookup(N::Int, A::Symbol, indexes::Symbol, itersym::Symbol)
    vars = [ :($(namedvar(indexes, i))[$(namedvar(itersym, i))]) for i = 1:N ]
    Expr(:escape, Expr(:ref, A, vars...))
end

# Convert to a linear index
macro nlinear(N, A, itersym)
    _nlinear(N, A, itersym)
end

function _nlinear(N::Int, A::Symbol, itersym::Symbol)
    Expr(:call, :linear, :($(esc(A))), [Expr(:escape, namedvar(itersym, i)) for i = 1:N]...)
end

function _nlinear(N::Int, A::Symbol, ex::Expr)
    Expr(:call, :linear, :($(esc(A))), [Expr(:escape, inlineanonymous(ex, i)) for i = 1:N]...)
end

namedvar(base::Symbol, ext) = symbol(string(base)*"_"*string(ext))

macro indexedvariable(N, sym)
    _indexedvariable(N, sym)
end

# _indexedvariable(ex::Expr, sym::Symbol) = _indexedvariable(eval(ex), sym)

function _indexedvariable(i::Integer, sym::Symbol)
    nv = namedvar(sym, i)
    :($(esc(nv)))
end

linear(A::Array, i1::Integer) = A, i1
linear(A::Array, i1::Integer, i2::Integer) = A, i1+size(A,1)*(i2-1)
linear(A::Array, i1::Integer, i2::Integer, i3::Integer) = A, i1+size(A,1)*(i2-1+size(A,2)*(i3-1))
linear(A::Array, i1::Integer, i2::Integer, i3::Integer, i4::Integer) = A, i1+size(A,1)*(i2-1+size(A,2)*(i3-1+size(A,3)*(i4-1)))
linear(A::Array, i1::Integer, i2::Integer, i3::Integer, i4::Integer, i5::Integer) = A, i1+size(A,1)*(i2-1+size(A,2)*(i3-1+size(A,3)*(i4-1+size(A,4)*(i5-1))))

linear{T}(s::SubArray{T,1}, i::Integer) =
    s.parent, s.first_index + (i-1)*s.strides[1]
linear{T}(s::SubArray{T,2}, i::Integer, j::Integer) =
    s.parent, s.first_index + (i-1)*s.strides[1] + (j-1)*s.strides[2]
linear{T}(s::SubArray{T,3}, i::Integer, j::Integer, k::Integer) =
    s.parent, s.first_index + (i-1)*s.strides[1] + (j-1)*s.strides[2] + (k-1)*s.strides[3]
linear{T}(s::SubArray{T,4}, i::Integer, j::Integer, k::Integer, l::Integer) =
    s.parent, s.first_index + (i-1)*s.strides[1] + (j-1)*s.strides[2] + (k-1)*s.strides[3] + (l-1)*s.strides[4]
linear{T}(s::SubArray{T,5}, i::Integer, j::Integer, k::Integer, l::Integer, m::Integer) =
    s.parent, s.first_index + (i-1)*s.strides[1] + (j-1)*s.strides[2] + (k-1)*s.strides[3] + (l-1)*s.strides[4] + (m-1)*s.strides[5]


text:
using JSON

const BEGIN_MARKER = "###JULIA-WEBIO-CONFIG-BEGIN"
const END_MARKER = "###JULIA-WEBIO-CONFIG-END"
function install_ijulia_config()
    config_file = joinpath(homedir(), ".jupyter", "jupyter_notebook_config.py")
    if isfile(config_file)
        config_str = String(read(config_file))
    else
        mkpath(dirname(config_file))
        config_str = ""
    end

    # remove previous config
    config_str = replace(config_str, Regex("\n?" * BEGIN_MARKER * ".*" * END_MARKER * "\n?", "s") => "")

    config_str *= """

    $BEGIN_MARKER
    import sys, os
    if os.path.isfile($(repr(joinpath(dirname(@__FILE__), "jlstaticserve.py")))):
        sys.path.append($(repr(dirname(@__FILE__))))
        c = get_config()
        c.NotebookApp.nbserver_extensions = {
            "jlstaticserve": True
        }
    else:
        print("WebIO config in ~/.jupyter/jupyter_notebook_config.py but WebIO plugin not found")
    $END_MARKER
    """
    write(config_file, config_str)
    config_file_json = joinpath(homedir(), ".jupyter", "jupyter_notebook_config.json")

    if isfile(config_file_json)
        dict = try
            JSON.parse(read(config_file_json, String))
        catch err
            println(stderr, "Error parsing Jupyter config file $config_file_json - fix it and build again or delete it to enable WebIO")
            @goto jsondone
        end
        app = Base.@get! dict "NotebookApp" Dict()
        nbext = Base.@get! app "nbserver_extensions" Dict()
        nbext["jlstaticserve"] = true
        open(config_file_json, "w") do io
            prettyio = JSON.Writer.PrettyContext(io, 4)
            JSON.print(prettyio, dict)
        end
    end
    @label jsondone
end

function get_jupyter_datadir()
    try
        return readline(open(`jupyter --data-dir`))
    catch (e)
        # Is there a way to use Conda.jl if it's installed?
        @warn "Didn't detect Jupyter."
    end

    # Try guessing based on OS
    # https://jupyter.readthedocs.io/en/latest/projects/jupyter-directories.html
    if Sys.iswindows()
        # Is this right?
        return joinpath(homedir(), "%APPDATA%", "jupyter")
    elseif Sys.isapple()
        return joinpath(homedir(), "Library", "Jupyter")
    else
        # Maybe need to check XDG_DATA_HOME environment variable?
        return joinpath(homedir(), ".local", "share", "jupyter")
    end
end

"""
Install the Jupyter WebIO notebook extension.
"""
function install_webio_nbextension()
    extension_dir = joinpath(get_jupyter_datadir(), "nbextensions")
    mkpath(extension_dir)

    # I think the config dir is always ~/.jupyter, even on Windows.
    config_dir = joinpath(homedir(), ".jupyter", "nbconfig")
    mkpath(config_dir)
    config_file_json = joinpath(config_dir, "notebook.json")

    # Copy the nbextension files.
    @info "Copying WebIO nbextension files to $(extension_dir)."
    cp(
        joinpath(@__DIR__, "../packages/jupyter-notebook-provider/dist"),
        joinpath(extension_dir, "webio"),
        ; force=true
    )

    # Enable the notebook extension.
    config_data = Dict()
    if isfile(config_file_json)
        config_data = try
            JSON.parse(read(config_file_json, String))
        catch err
            println(stderr, "Error parsing Jupyter config file $config_file_json - fix it and build again or delete it to enable WebIO")
            return
        end
    end
    config_data["load_extensions"] = get(config_data, "load_extensions", Dict())
    config_data["load_extensions"]["webio/main"] = true
    open(config_file_json, "w") do io
        prettyio = JSON.Writer.PrettyContext(io, 4)
        JSON.print(prettyio, config_data)
    end
end

install_ijulia_config()
install_webio_nbextension()


text:
function node( obj::tikzObj, x::Float64, y::Float64; lw = 0.0, color = "",
               linestyle="", style=tikzStyle(), extracommands="", name="", text="" )

    n = length(x)

    # print start
    @printf( obj.fID, "\\node[" )
    if length(color) > 0
        @printf( obj.fID, "draw=%s, ", color )
    end
    if lw > 0.0
        @printf( obj.fID, "line width=%3.2fpt, ", lw )
    end
    if style.set
        @printf( obj.fID, "%s, ", style.name )
    end
    @printf( obj.fID, "%s] ", extracommands )

    if length(name) > 0
        @printf( obj.fID, " (%s) ", name )
    end
    @printf( obj.fID, " at ( %6.4f, %6.4f ) {%s};\n\n", x, y, text )

end

function node( obj::tikzObj, objAx::axisObj, x::Float64, y::Float64; lw = 0.0, color = "",
               linestyle="", style=tikzStyle(), extracommands="", name="", text="" )

    n = length(x)

    # print start
    @printf( obj.fID, "\\node[" )
    if length(color) > 0
        @printf( obj.fID, "draw=%s, ", color )
    end
    if lw > 0.0
        @printf( obj.fID, "line width=%3.2fpt, ", lw )
    end
    if style.set
        @printf( obj.fID, "%s, ", style.name )
    end
    @printf( obj.fID, "%s] ", extracommands )

    if length(name) > 0
        @printf( obj.fID, " (%s) ", name )
    end
    @printf( obj.fID, " at (axis cs: %6.4f, %6.4f ) {%s};\n\n", x, y, text )

end

max_stars_repo_name: mopg/julia-tikz

max_issues_repo_name: mopg/julia-tikz

max_forks_repo_name: mopg/julia-tikz


text:
ipft = [2 1 1]
L = [1 1.5 2]

ncoh = length(ipft)

lr = [0.045 0.056 0.056]
lt = [0.025 0.046 0.046]
orient = [0.2 0.1][ipft]

theta = 15
S0_beam = 0.8
S0_diff = 1 - S0_beam
alb_ground = 0.1

nu = lr + lt
omega = nu .+ 0.25 .* (lr .- lt) .* (1 .- orient) .^ 2 ./ (2 .* nu)

phi1 = ...
phi2 = ...

a = [1 - alb_ground ]

max_stars_repo_name: femeunier/fortebaseline

max_issues_repo_name: femeunier/fortebaseline

max_forks_repo_name: femeunier/fortebaseline


text:
include("AirtableTests.jl")
AirtableTests.retest()

using Airtable
Airtable.retest()


text:
using SweepContractor
using ..ITensorAutoHOOT
using SweepContractor: TensorNetwork, LabelledTensorNetwork
using ..ITensorAutoHOOT: SubNetwork

function ITensor_networks(TN::TensorNetwork)
  SweepContractor.sort!(TN)
  index_dict = Dict()
  function itensor(i, t)
    # Construct indices
    inds = []
    for (dim, j) in enumerate(t.adj)
      label = sort([i, j])
      if !haskey(index_dict, label)
        s = size(t.arr)[dim]
        index_dict[label] = Index(s, string(label))
      end
      push!(inds, index_dict[label])
    end
    # build the tensor
    return ITensor(t.arr, inds...)
  end
  inetwork = [itensor(i, t) for (i, t) in enumerate(TN)]
  return inetwork
end

function ITensor_networks(LTN::LabelledTensorNetwork)
  return ITensor_networks(SweepContractor.delabel(LTN))
end

function line_network(network::Vector)
  if length(network) <= 2
    return SubNetwork(network)
  end
  return SubNetwork(line_network(network[1:(end - 1)]), network[end])
end


text:
function getCoefficientsRadIIA2(T=Float64)
    a = [[5//12  -1//12]
         [3//4    1//4 ]]
    b = [3//4, 1//4]
    c = [1//3, 1//1]

    CoefficientsRK(T, :RadIIA2, 3, a, b, c)
end

function getCoefficientsRadIIA3(T=Float64)
    a = @dec128 [
          [11/45 -  7*√6/360    37/225-169*√6/1800   -2/225+√6/75]
          [37/225+169*√6/1800   11/45 +  7*√6/360    -2/225-√6/75]
          [ 4/9  -    √6/36      4/9  +    √6/36      1/9        ]
        ]
    b = @dec128 [4/9-√6/36,   4/9+√6/36,  1/9 ]
    c = @dec128 [2/5-√6/10,   2/5+√6/10,  1   ]

    CoefficientsRK(T, :RadIIA3, 5, a, b, c)
end