text:
module HybridArrays

import Base: convert,
    copy,
    copyto!,
    hcat,
    vcat,
    similar,
    axes,
    getindex,
    dataids,
    promote_rule,
    setindex!,
    size,
    length,
    +,
    -,
    *

import Base.Array

using StaticArrays
using StaticArrays: Dynamic
import StaticArrays: _setindex!_scalar, Size

using LinearAlgebra

@generated function hasdynamic(::Type{Size}) where Size<:Tuple
    for s ∈ Size.parameters
        if isa(s, Dynamic)
            return true
        end
    end
    return false
end

function all_dynamic_fixed_val(::Type{Size}, inds::StaticArrays.StaticIndexing{<:Union{Int, AbstractArray, Colon}}...) where Size<:Tuple
    return all_dynamic_fixed_val(Size, map(StaticArrays.unwrap, inds)...)
end

@generated function all_dynamic_fixed_val(::Type{Size}, inds::Union{Int, AbstractArray, Colon}...) where Size<:Tuple
    all_fixed = true
    for (i, param) in enumerate(Size.parameters)

        destatizing = (inds[i] <: AbstractArray && !(
            inds[i] <: StaticArray ||
            inds[i] <: Base.Slice{<:StaticArray} ||
            inds[i] <: SOneTo ||
            inds[i] <: Base.Slice{<:SOneTo}))

        nonstatizing = inds[i] == Colon || destatizing

        if destatizing || (isa(param, Dynamic) && nonstatizing)
            all_fixed = false
            break
        end

    end

    if all_fixed
        return Val(:dynamic_fixed_true)
    else
        return Val(:dynamic_fixed_false)
    end
end

@generated function tuple_nodynamic_prod(::Type{Size}) where Size<:Tuple
    i = 1
    for s ∈ Size.parameters
        if !isa(s, Dynamic)
            i *= s
        end
    end
    return i
end

# conversion of SizedArray from StaticArrays.jl
"""
    HybridArray{Tuple{dims...}}(array)

Wraps an `AbstractArray` with a combination of static and dynamic sizes,
so to take advantage of the (faster) methods defined by the StaticArrays
package. The size is checked once upon construction to determine if the
number of elements (`length`) match, but the array may be reshaped.
"""
struct HybridArray{S<:Tuple, T, N, M, TData<:AbstractArray{T,M}} <: AbstractArray{T, N}
    data::TData

    function HybridArray{S, T, N, M, TData}(a::TData) where {S, T, N, M, TData<:AbstractArray{T,M}}
        tnp = tuple_nodynamic_prod(S)
        lena = length(a)
        dynamic_nodivisible = hasdynamic(S) && tnp != 0 && mod(lena, tnp) != 0
        nodynamic_notequal = !hasdynamic(S) && lena != StaticArrays.tuple_prod(S)
        if nodynamic_notequal || dynamic_nodivisible || (tnp == 0 && lena != 0)
            error("Dimensions $(size(a)) don't match static size $S")
        end
        new{S,T,N,M,TData}(a)
    end

    function HybridArray{S, T, N, 1}(::UndefInitializer) where {S, T, N}
        new{S, T, N, 1, Array{T, 1}}(Array{T, 1}(undef, StaticArrays.tuple_prod(S)))
    end
    function HybridArray{S, T, N, N}(::UndefInitializer) where {S, T, N}
        new{S, T, N, N, Array{T, N}}(Array{T, N}(undef, size_to_tuple(S)...))
    end
end

@inline HybridArray{S,T,N}(a::TData) where {S,T,N,M,TData<:AbstractArray{T,M}} = HybridArray{S,T,N,M,TData}(a)
@inline HybridArray{S,T}(a::TData) where {S,T,M,TData<:AbstractArray{T,M}} = HybridArray{S,T,StaticArrays.tuple_length(S),M,TData}(a)
@inline HybridArray{S}(a::TData) where {S,T,M,TData<:AbstractArray{T,M}} = HybridArray{S,T,StaticArrays.tuple_length(S),M,TData}(a)

@inline HybridArray{S,T,N}(::UndefInitializer) where {S,T,N} = HybridArray{S,T,N,N}(undef)
@inline HybridArray{S,T}(::UndefInitializer) where {S,T} = HybridArray{S,T,StaticArrays.tuple_length(S),StaticArrays.tuple_length(S)}(undef)
@generated function (::Type{HybridArray{S,T,N,M,TData}})(x::NTuple{L,Any}) where {S,T,N,M,TData,L}
    if L != StaticArrays.tuple_prod(S)
        error("Dimension mismatch")
    end
    exprs = [:(a[$i] = x[$i]) for i = 1:L]
    return quote
        $(Expr(:meta, :inline))
        a = HybridArray{S,T,N,M}(undef)
        @inbounds $(Expr(:block, exprs...))
        return a
    end
end

@inline HybridArray{S,T,N}(x::Tuple) where {S,T,N} = HybridArray{S,T,N,N,Array{T,N}}(x)
@inline HybridArray{S,T}(x::Tuple) where {S,T} = HybridArray{S,T,StaticArrays.tuple_length(S)}(x)
@inline HybridArray{S}(x::NTuple{L,T}) where {S,T,L} = HybridArray{S,T}(x)

HybridVector{S,T,M} = HybridArray{Tuple{S},T,1,M}
@inline HybridVector{S}(a::TData) where {S,T,M,TData<:AbstractArray{T,M}} = HybridArray{Tuple{S},T,1,M,TData}(a)
@inline HybridVector{S}(x::NTuple{L,T}) where {S,T,L} = HybridArray{Tuple{S},T,1,1,Vector{T}}(x)

HybridMatrix{S1,S2,T,M} = HybridArray{Tuple{S1,S2},T,2,M}
@inline HybridMatrix{S1,S2}(a::TData) where {S1,S2,T,M,TData<:AbstractArray{T,M}} = HybridArray{Tuple{S1,S2},T,2,M,TData}(a)
@inline HybridMatrix{S1,S2}(x::NTuple{L,T}) where {S1,S2,T,L} = HybridArray{Tuple{S1,S2},T,2,2,Matrix{T}}(x)

export HybridArray, HybridMatrix, HybridVector

include("ssubarray.jl")

include("abstractarray.jl")
include("arraymath.jl")
include("broadcast.jl")
include("convert.jl")
include("indexing.jl")
include("linalg.jl")
include("utils.jl")

end # module


text:
import Base: unbox, box, reinterpret

for (utype,itype,ftype) in ((:UInt32, :Int32, :Float32), (:UInt64, :Int64, :Float64))
    @eval begin
        reinterpret(::Type{Signed},   ::Type{$utype}) = ($itype)
        reinterpret(::Type{Signed},   ::Type{$ftype}) = ($itype)
        reinterpret(::Type{Unsigned}, ::Type{$itype}) = ($utype)
        reinterpret(::Type{Unsigned}, ::Type{$ftype}) = ($utype)
        reinterpret(::Type{AbstractFloat}, ::Type{$utype}) = ($ftype)
        reinterpret(::Type{AbstractFloat}, ::Type{$itype}) = ($ftype)

        reinterpret(::Type{Signed}, x::($utype))   = box($itype,unbox($utype,x))
        reinterpret(::Type{Signed}, x::($itype))   = x
        reinterpret(::Type{Signed}, x::($ftype))   = box($itype,unbox($ftype,x))

        reinterpret(::Type{Unsigned}, x::($utype)) = x
        reinterpret(::Type{Unsigned}, x::($itype)) = box($utype,unbox($itype,x))

        reinterpret(::Type{AbstractFloat}, x::($utype)) = box($ftype,unbox($utype,x))
        reinterpret(::Type{AbstractFloat}, x::($itype)) = box($ftype,unbox($itype,x))
        reinterpret(::Type{AbstractFloat}, x::($ftype)) = x
    end
end

ext: jl


text:
using Dash
using DashHtmlComponents
using DashCoreComponents

app = dash()

app.layout = html_div() do
    dcc_dropdown(id = "input",
        options = [
            (label="regular", value="regular"),
            (label="prevent", value="prevent"),
            (label="no_update1", value="no_update1"),
            (label="no_update2", value="no_update2"),
        ]
    ),
    html_div(id="output1"),
    html_div(id="output2"),
    html_div(id="regular_output")
end
callback!(app, callid"input.value => regular_output.children") do input
    return input
end
callback!(app, callid"input.value => output1.children, output2.children") do input
    input == "prevent" && throw(PreventUpdate())
    if input == "no_update1"
        return (no_update(), input)
    end
    if input == "no_update2"
        return (input, no_update())
    end
    return (input, input)
end

run_server(app)


text:
@testset "Logarith of power series matrices" begin

    T, t = Nemo.PowerSeriesRing(Nemo.QQ, 300, "t"; model=:capped_absolute)

    for d in 1:5
        diag_elements = [1 + t * random_ps(T) for i in 1:d]
        S = Nemo.MatrixSpace(T, d, d)
        M = S([(i == j ? diag_elements[i] : T(0)) for i = 1:d, j = 1:d])
        result = ps_matrix_log(M)
        correct = S([(i == j ? log(diag_elements[i]) : T(0)) for i = 1:d, j = 1:d])
        @test result == correct
    end

    for c in 1:5
        f = random_ps(T) * t
        S = Nemo.MatrixSpace(T, 2, 2)
        m = S([
            T(1) f;
            T(0) T(1)
        ])
        correct = S([
            T(0) f;
            T(0) T(0)
        ])
        @test ps_matrix_log(m) == correct
    end

    for c in 1:5
        f, g = random_ps(T) * t, random_ps(T) * t
        S = Nemo.MatrixSpace(T, 3, 3)
        m = S([
            T(1) f    g;
            T(0) T(1) f;
            T(0) T(0) T(1)
        ])
        correct = S([
            T(0) f    (-1 // 2) * f^2 + g;
            T(0) T(0) f;
            T(0) T(0) T(0)
        ])
        @test ps_matrix_log(m) == correct
    end

end


text:
using ThreadedIterables;
using Test;

if VERSION >= v"1.2"
	using OffsetArrays;
end

using Random;

include("utils.jl");

include("testForeach.jl");
include("testMap.jl");
include("testReduce.jl");
include("testMapreduce.jl");
include("testMapfoldl.jl");
include("testMapfoldr.jl");

include("testThreaded.jl");

@testset "ThreadedIterables.jl" begin
	@testset "foreach" begin
		@testset "tforeach(f, c)" begin
			for i in 1:1000
				@test testTforeach();
			end
		end
		if VERSION >= v"1.4"
			@testset "tforeach(f, c::OffsetArray)" begin
				@test testTforeachOffset();
			end
		end
		@testset "tforeach(f, [])" begin
			@test testTforeachEmpty();
		end
		@testset "tforeach(f, c1, c2)" begin
			for i in 1:1000
				@test testTforeachMultiple();
			end
		end
		if VERSION >= v"1.2"
			@testset "tforeach(f, c1::OffsetArray, c2::OffsetArray)" begin
				@test testTforeachMultipleOffset();
			end
		end
		@testset "tforeach(f, [], [])" begin
			@test testTforeachMultipleEmpty();
		end
	end
	@testset "map" begin
		@testset "tmap(f, c)" begin
			@testset "Type stable" begin
				for i in 1:1000
					@test testTmapStable();
				end
			end
			@testset "Type unstable" begin
				for i in 1:1000
					@test testTmapUnstable();
				end
			end
		end
		if VERSION >= v"1.2"
			@testset "tmap(f, c::OffsetArray)" begin
				@test testTmapOffset();
			end
		end
		@testset "tmap(f, [])" begin
			@test testTmapEmpty();
		end
		@testset "tmap(f, c1, c2)" begin
			@testset "Type stable" begin
				for i in 1:1000
					@test testTmapStableMultiple();
				end
			end
			@testset "Type unstable" begin
				for i in 1:1000
					@test testTmapUnstableMultiple();
				end
			end
		end
		if VERSION >= v"1.2"
			@testset "tmap(f, c1::OffsetArray, c2::OffsetArray)" begin
				@test testTmapMultipleOffset();
			end
		end
		@testset "tmap(f, [], [])" begin
			@test testTmapMultipleEmpty();
		end
		@testset "tmap!(f, destination, collection)" begin
			for i in 1:1000
				@test testTmap!();
			end
		end
		@testset "tmap!(f, destination, [])" begin
			@test testTmap!Empty();
		end
		@testset "tmap!(f, destination, collection1, collection2)" begin
			for i in 1:1000
				@test testTmap!Multiple();
			end
		end
		@testset "tmap!(f, destination, [], [])" begin
			@test testTmap!MultipleEmpty();
		end
		@testset "tmap(f, c1, c2) where length(c1) != length(c2)" begin
			@test_throws DimensionMismatch testTmapMultipleDimensionMismatch();
		end
	end
	@testset "reduce" begin
		@testset "reduce(op, v0, itr)" begin
			for i in 1:1000
				@test testTreduce1();
			end
		end
		if VERSION >= v"1.2"
			@testset "reduce(op, v0, itr::OffsetArray)" begin
				@test testTreduce1Offset();
			end
		end
		@testset "reduce(op, v0, [])" begin
			@test testTreduce1Empty();
		end
		@testset "reduce(op, itr)" begin
			for i in 1:1000
				@test testTreduce2();
			end
		end
		if VERSION >= v"1.2"
			@testset "reduce(op, itr::OffsetArray)" begin
				@test testTreduce2Offset();
			end
		end
		@testset "reduce(op, [])" begin
			@test_throws ArgumentError testTreduce2Empty();
		end
	end
	@testset "mapreduce" begin
		@testset "tmapreduce(f, op, itr; init)" begin
			@testset "Type stable" begin
				for i in 1:1000
					@test testTmapreduce1Stable();
				end
			end
			@testset "Type unstable" begin
				for i in 1:1000
					@test testTmapreduce1Unstable();
				end
			end
		end
		if VERSION >= v"1.2"
			@testset "tmapreduce(f, op, itr::OffsetArray; init)" begin
				@test testTmapreduce1Offset();
			end
		end
		@testset "tmapreduce(f, op, []; init)" begin
			@test testTmapreduce1Empty();
		end
		if VERSION >= v"1.2"
			@testset "tmapreduce(f, op, itr1, itr2; init)" begin
				@testset "Type stable" begin
					for i in 1:1000
						@test testTmapreduce1StableMultiple();
					end
				end
				@testset "Type unstable" begin
					for i in 1:1000
						@test testTmapreduce1UnstableMultiple();
					end
				end
			end
			@testset "tmapreduce(f, op, itr1::OffsetArray, itr2::OffsetArray; init)" begin
				@test testTmapreduce1MultipleOffset();
			end
			@testset "tmapreduce(f, op, [], []; init)" begin
				@test testTmapreduce1MultipleEmpty();
			end
		end
		@testset "tmapreduce(f, op, itr)" begin
			@testset "Type stable" begin
				for i in 1:1000
					@test testTmapreduce2Stable();
				end
			end
			@testset "Type unstable" begin
				for i in 1:1000
					@test testTmapreduce2Unstable();
				end
			end
		end
		if VERSION >= v"1.2"
			@testset "tmapreduce(f, op, itr::OffsetArray)" begin
				@test testTmapreduce2Offset();
			end
		end
		@testset "tmapreduce(f, op, [])" begin
			@test_throws Union{ArgumentError,MethodError} testTmapreduce2Empty();
		end
		if VERSION >= v"1.2"
			@testset "tmapreduce(f, op, itr1, itr2)" begin
				@testset "Type stable" begin
					for i in 1:1000
						@test testTmapreduce2StableMultiple();
					end
				end
				@testset "Type unstable" begin
					for i in 1:1000
						@test testTmapreduce2UnstableMultiple();
					end
				end
			end
			@testset "tmapreduce(f, op, itr1::OffsetArray, itr2::OffsetArray)" begin
				@test testTmapreduce2MultipleOffset();
			end
			@testset "tmapreduce(f, op, [], [])" begin
				@test_throws Union{ArgumentError,MethodError} testTmapreduce2MultipleEmpty();
			end
		end
		@testset "maptreduce(f, op, itr; init)" begin
			for i in 1:1000
				@test testMaptreduce1();
			end
		end
		if VERSION >= v"1.2"
			@testset "maptreduce(f, op, itr::OffsetArray; init)" begin
				@test testMaptreduce1Offset();
			end
		end
		@testset "maptreduce(f, op, []; init)" begin
			@test testMaptreduce1Empty();
		end
		if VERSION >= v"1.2"
			@testset "maptreduce(f, op, itr1, itr2; init)" begin
				for i in 1:1000
					@test testMaptreduce1Multiple();
				end
			end
			@testset "maptreduce(f, op, itr1::OffsetArray, itr2::OffsetArray; init)" begin
				@test testMaptreduce1MultipleOffset();
			end
			@testset "maptreduce(f, op, [], []; init)" begin
				@test testMaptreduce1MultipleEmpty();
			end
		end
		@testset "maptreduce(f, op, itr)" begin
			for i in 1:1000
				@test testMaptreduce2();
			end
		end
		if VERSION >= v"1.2"
			@testset "maptreduce(f, op, itr::OffsetArray)" begin
				@test testMaptreduce2Offset();
			end
		end
		@testset "maptreduce(f, op, [])" begin
			@test_throws ArgumentError testMaptreduce2Empty();
		end
		if VERSION >= v"1.2"
			@testset "maptreduce(f, op, itr1, itr2)" begin
				for i in 1:1000
					@test testMaptreduce2Multiple();
				end
			end
			@testset "maptreduce(f, op, itr1::OffsetArray, itr2::OffsetArray)" begin
				@test testMaptreduce2MultipleOffset();
			end
			@testset "maptreduce(f, op, [], [])" begin
				@test_throws ArgumentError testMaptreduce2MultipleEmpty();
			end
		end
		@testset "tmaptreduce(f, op, itr; init)" begin
			for i in 1:1000
				@test testTmaptreduce1();
			end
		end
		if VERSION >= v"1.2"
			@testset "tmaptreduce(f, op, itr::OffsetArray; init)" begin
				@test testTmaptreduce1Offset();
			end
		end
		@testset "tmaptreduce(f, op, []; init)" begin
			@test testTmaptreduce1Empty();
		end
		if VERSION >= v"1.2"
			@testset "tmaptreduce(f, op, itr1, itr2; init)" begin
				for i in 1:1000
					@test testTmaptreduce1Multiple();
				end
			end
			@testset "tmaptreduce(f, op, itr1::OffsetArray, itr2::OffsetArray; init)" begin
				@test testTmaptreduce1MultipleOffset();
			end
			@testset "tmaptreduce(f, op, [], []; init)" begin
				@test testTmaptreduce1MultipleEmpty();
			end
		end
		@testset "tmaptpreduce(f, op, itr)" begin
			for i in 1:1000
				@test testTmaptreduce2();
			end
		end
		if VERSION >= v"1.2"
			@testset "tmatpreduce(f, op, itr::OffsetArray)" begin
				@test testTmaptreduce2Offset();
			end
		end
		@testset "tmaptpreduce(f, op, [])" begin
			@test_throws ArgumentError testTmaptreduce2Empty();
		end
		if VERSION >= v"1.2"
			@testset "tmaptpreduce(f, op, itr1, itr2)" begin
				for i in 1:1000
					@test testTmaptreduce2Multiple();
				end
			end
			@testset "tmatpreduce(f, op, itr1::OffsetArray, itr2::OffsetArray)" begin
				@test testTmaptreduce2MultipleOffset();
			end
			@testset "tmaptpreduce(f, op, [], [])" begin
				@test_throws ArgumentError testTmaptreduce2MultipleEmpty();
			end
		end
	end
	@testset "mapfoldl" begin
		@testset "tmapfoldl(f, op, itr; init)" begin
			@testset "Type stable" begin
				for i in 1:1000
					@test testTmapfoldl1Stable();
				end
			end
			@testset "Type unstable" begin
				for i in 1:1000
					@test testTmapfoldl1Unstable();
				end
			end
		end
		if VERSION >= v"1.2"
			@testset "tmapfoldl(f, op, itr::OffsetArray; init)" begin
				@test testTmapfoldl1Offset();
			end
		end
		@testset "tmapfoldl(f, op, []; init)" begin
			@test testTmapfoldl1Empty();
		end
		@testset "tmapfoldl(f, op, itr)" begin
			@testset "Type stable" begin
				for i in 1:1000
					@test testTmapfoldl2Stable();
				end
			end
			@testset "Type unstable" begin
				for i in 1:1000
					@test testTmapfoldl2Unstable();
				end
			end
		end
		if VERSION >= v"1.2"
			@testset "tmapfoldl(f, op, itr::OffsetArray)" begin
				@test testTmapfoldl2Offset();
			end
		end
		@testset "tmapfoldl(f, op, [])" begin
			@test_throws Union{ArgumentError,MethodError} testTmapfoldl2Empty();
		end
	end
	@testset "mapfoldr" begin
		@testset "tmapfoldr(f, op, itr; init)" begin
			@testset "Type stable" begin
				for i in 1:1000
					@test testTmapfoldr1Stable();
				end
			end
			@testset "Type unstable" begin
				for i in 1:1000
					@test testTmapfoldr1Unstable();
				end
			end
		end
		if VERSION >= v"1.2"
			@testset "tmapfoldr(f, op, itr::OffsetArray; init)" begin
				@test testTmapfoldr1Offset();
			end
		end
		@testset "tmapfoldr(f, op, []; init)" begin
			@test testTmapfoldr1Empty();
		end
		@testset "tmapfoldr(f, op, itr)" begin
			@testset "Type stable" begin
				for i in 1:1000
					@test testTmapfoldr2Stable();
				end
			end
			@testset "Type unstable" begin
				for i in 1:1000
					@test testTmapfoldr2Unstable();
				end
			end
		end
		if VERSION >= v"1.2"
			@testset "tmapfoldr(f, op, itr::OffsetArray)" begin
				@test testTmapfoldr2Offset();
			end
		end
		@testset "tmapfoldr(f, op, [])" begin
			@test_throws Union{ArgumentError,MethodError} testTmapfoldr2Empty();
		end
	end
	@testset "@threaded" begin
		@test testThreadedForeach();
		@test testThreadedMapfoldl();
		@test testThreadedMapfoldr();
		@test testThreadedMap();
		@test testThreadedMap!();
		@test testThreadedMapreduce();
		@test testThreadedReduce();
		@test testThreadedOther();
		@test testThreadedNested();
	end
end


text:
function GeometricDistribution(arg0::jdouble)
    return GeometricDistribution((jdouble,), arg0)
end

function cumulative_probability(obj::GeometricDistribution, arg0::jint)
    return jcall(obj, "cumulativeProbability", jdouble, (jint,), arg0)
end

function get_numerical_mean(obj::GeometricDistribution)
    return jcall(obj, "getNumericalMean", jdouble, ())
end

function get_numerical_variance(obj::GeometricDistribution)
    return jcall(obj, "getNumericalVariance", jdouble, ())
end

function get_probability_of_success(obj::GeometricDistribution)
    return jcall(obj, "getProbabilityOfSuccess", jdouble, ())
end

function get_support_lower_bound(obj::GeometricDistribution)
    return jcall(obj, "getSupportLowerBound", jint, ())
end

function get_support_upper_bound(obj::GeometricDistribution)
    return jcall(obj, "getSupportUpperBound", jint, ())
end

function inverse_cumulative_probability(obj::GeometricDistribution, arg0::jdouble)
    return jcall(obj, "inverseCumulativeProbability", jint, (jdouble,), arg0)
end

function is_support_connected(obj::GeometricDistribution)
    return jcall(obj, "isSupportConnected", jboolean, ())
end

function log_probability(obj::GeometricDistribution, arg0::jint)
    return jcall(obj, "logProbability", jdouble, (jint,), arg0)
end

function probability(obj::GeometricDistribution, arg0::jint)
    return jcall(obj, "probability", jdouble, (jint,), arg0)
end


text:
gensym_unhashed(s::Symbol) = Symbol(Base.replace(string(gensym(s)), "#"=>""))

function_walk(f, expr) = postwalk(expr) do x
    @capture(x, g_(a__)) && issymbol(g) ? f(g, a) : x
end

symbol_walk(f, ex) = postwalk(x -> issymbol(x) ? f(x) : x, ex)

atom_walk(f, x) = walk(x, x -> (@capture(x, y_(ys__)) && issymbol(y)) ? f(x) : atom_walk(f, x) , f)

issymbol(x) = x isa Symbol && Base.isidentifier(x)
isfunction(x) = @capture(x, s_(xs__)) && s isa Symbol && Base.isidentifier(s)

function symbols(x::Expr)
    ls = Symbol[]
    symbol_walk(x) do s
        push!(ls, s)
        s
    end
    Base.unique(ls)
end
symbols(x::Symbol) = [x]
symbols(::Number) = []

function (R::PolyRing)(x::Union{Symbol,Expr})
    vs = [:($(Symbol(string(g))) = $g) for g in gens(R)]
    qq = quote
        let $(vs...)
            $x
        end
    end
    eval(qq)
end

function (R::PolyRing)(x::Number)
    base_ring(R)(x)
end

mroots(p::PolyElem) = [-(Nemo.coeff(f, 0)/lead(f), m) for (f, m) in Nemo.factor(p)]
roots(p::PolyElem) = [-Nemo.coeff(f, 0)//lead(f) for (f, _) in Nemo.factor(p)]

fallfactorial(n, j) = j == 0 ? 1 : prod(n - i for i in 0:j-1)


text:
@testset "Polyhedral" begin
    @testset "affine + torus solutions" begin
        @var x y
        f = System([2y + 3 * y^2 - x * y^3, x + 4 * x^2 - 2 * x^3 * y])

        tracker, starts = polyhedral(f; only_torus = false)
        @test length(collect(starts)) == 8
        @test count(is_success, track.(tracker, starts)) == 6

        tracker, starts = polyhedral(f; only_torus = true)
        @test length(collect(starts)) == 3
        @test count(is_success, track.(tracker, starts)) == 3

        tracker, starts = polyhedral(
            f;
            only_torus = true,
            tracker_options = (automatic_differentiation = 3,),
        )
        @test length(collect(starts)) == 3
        @test count(is_success, track.(tracker, starts)) == 3
    end

    @testset "only torus" begin
        @var x₁ x₂ s

        HC_I = [
            x₁^3 * s^15 + x₁ * x₂ * s + x₂^3 + s^12,
            x₁^2 * s^9 + x₁ * x₂^2 + x₂ * s^3,
            x₁^2 * x₂ * s^5 + x₁ * s^8 + x₂^2,
        ]
        F = System(HC_I)
        _, starts = solver_startsolutions(F; start_system = :polyhedral, only_torus = false)
        @test length(collect(starts)) == paths_to_track(F; only_torus = false) == 92
        _, starts = solver_startsolutions(F; start_system = :polyhedral, only_torus = true)
        @test length(collect(starts)) == paths_to_track(F; only_torus = true) == 54
    end

    @testset "cyclic" begin
        tracker, starts = polyhedral(cyclic(5))
        res = track.(tracker, starts)
        @test count(is_success, res) == 70

        tracker, starts = polyhedral(cyclic(7))
        res = track.(tracker, starts)
        @test count(is_success, res) == 924
    end
end


text:
struct NNStatistics
    units::String
    arealUnits::String
    area::Polygon
    observedMeanDistance::Real
    expectedMeanDistance::Real
    nearestNeighbourIndex::Real
    pointsCount::Real
    zScore::Real
end

"""
    analysis(data::F, area::Union{P, Nothing}=nothing, units::String="kilometers") where {F <: AbstractFeatureCollection, P <: AbstractPolygon}

Nearest Neighbor Analysis calculates an index based the average distances
between points in the dataset, thereby providing inference as to whether the
data is clustered, dispersed, or randomly distributed within the study area.

It returns a Polygon of the study area, with the results of
the analysis attached as part of of the `nearestNeighborAnalysis` property
of the study area's `properties`. The attached
[_z_-score](https://en.wikipedia.org/wiki/Standard_score) indicates how many
standard deviations above or below the expected mean distance the data's
observed mean distance is. The more negative, the more clustered. The more
positive, the more evenly dispersed. A _z_-score between -2 and 2 indicates
a seemingly random distribution. That is, within _p_ of less than 0.05, the
distribution appears statistically significantly neither clustered nor
dispersed.

**Remarks**

- Though the analysis will work on any FeatureCollection type, it
works best with Point collections.

- This analysis is **very** sensitive to the study area provided.
If no Polygon is passed as the study area, the function draws a box
around the data, which may distort the findings. This analysis works best
with a bounded area of interest within with the data is either clustered,
dispersed, or randomly distributed. For example, a city's subway stops may
look extremely clustered if the study area is an entire state. On the other
hand, they may look rather evenly dispersed if the study area is limited to
the city's downtown.
"""
function analysis(data::F, area::Union{P, Nothing}=nothing, units::String="kilometers") where {F <: AbstractFeatureCollection, P <: AbstractPolygon}
    area == nothing && (area = bbox_polygon(bbox(data)[1]))

    feats = []

    for ft in data.features
        push!(feats, centroid(ft.geometry))
    end

    n = length(feats)

    obsDist = []
    for (index, p) in enumerate(feats)
        other = filter((f, i) -> i !== index, feats)
        push!(obsDist, distance(p, nearestpoint(p, other).coordinates, units))
    end

    reduce((a, b) -> a + b, obsDist)

    popDensity = n / convert_area(area(area), "meters", units)
    expDist = 1 / (2 * sqrt(popDensity))
    var = 0.26136 / (sqrt(n * popDensity))

    return NNStatistics(units, units * "²", area, obsDist, expDist, obsDist / expDist, n, (obsDist - expDist) / var)

end


text:
#
# Types
#

"BSONType mirrors C enum bson_type_t."
primitive type BSONType 8 end # 1 byte

Base.convert(::Type{T}, t::BSONType) where {T<:Number} = reinterpret(UInt8, t)
Base.convert(::Type{BSONType}, n::T) where {T<:Number} = reinterpret(BSONType, n)
BSONType(u::UInt8) = convert(BSONType, u)

#
# Constants for BSONType
#

const BSON_TYPE_EOD        = BSONType(0x00)
const BSON_TYPE_DOUBLE     = BSONType(0x01)
const BSON_TYPE_UTF8       = BSONType(0x02)
const BSON_TYPE_DOCUMENT   = BSONType(0x03)
const BSON_TYPE_ARRAY      = BSONType(0x04)
const BSON_TYPE_BINARY     = BSONType(0x05)
const BSON_TYPE_UNDEFINED  = BSONType(0x06)
const BSON_TYPE_OID        = BSONType(0x07)
const BSON_TYPE_BOOL       = BSONType(0x08)
const BSON_TYPE_DATE_TIME  = BSONType(0x09)
const BSON_TYPE_NULL       = BSONType(0x0A)
const BSON_TYPE_REGEX      = BSONType(0x0B)
const BSON_TYPE_DBPOINTER  = BSONType(0x0C)
const BSON_TYPE_CODE       = BSONType(0x0D)
const BSON_TYPE_SYMBOL     = BSONType(0x0E)
const BSON_TYPE_CODEWSCOPE = BSONType(0x0F)
const BSON_TYPE_INT32      = BSONType(0x10)
const BSON_TYPE_TIMESTAMP  = BSONType(0x11)
const BSON_TYPE_INT64      = BSONType(0x12)
const BSON_TYPE_DECIMAL128 = BSONType(0x13)
const BSON_TYPE_MAXKEY     = BSONType(0x7F)
const BSON_TYPE_MINKEY     = BSONType(0xFF)

"""
BSONIter mirrors C struct bson_iter_t and can be allocated in the stack.

According to [libbson documentation](http://mongoc.org/libbson/current/bson_iter_t.html),
it is meant to be used on the stack and can be discarded at any time
as it contains no external allocation.
The contents of the structure should be considered private
and may change between releases, however the structure size will not change.

Inspecting its size in C, we get:

```c
sizeof(bson_iter_t) == 80
```
"""
primitive type BSONIter 80 * 8 end # 80 bytes

"""
Mirrors C struct `bson_oid_t`.

`BSONObjectId` instances addresses are passed
to libbson/libmongoc API using `Ref(oid)`,
and are owned by the Julia process.

```c
typedef struct {
   uint8_t bytes[12];
} bson_oid_t;
```
"""
primitive type BSONObjectId 12 * 8 end # 12 bytes

function BSONObjectId()
    oid_ref = Ref{BSONObjectId}()
    bson_oid_init(oid_ref, C_NULL)
    return oid_ref[]
end

function BSONObjectId(oid_string::String)
    if !bson_oid_is_valid(oid_string)
        error("'$oid_string' is not a valid ObjectId.")
    end
    oid_ref = Ref{BSONObjectId}()
    bson_oid_init_from_string(oid_ref, oid_string)
    return oid_ref[]
end

"""
Mirrors C struct `bson_error_t` and can be allocated in the stack.

`BSONError` instances addresses are passed
to libbson/libmongoc API using `Ref(error)`,
and are owned by the Julia process.

```c
typedef struct {
   uint32_t domain;
   uint32_t code;
   char message[504];
} bson_error_t;
```
"""
mutable struct BSONError
    domain::UInt32
    code::UInt32
    message::NTuple{504, UInt8}

    BSONError() = new(0, 0, tuple(zeros(UInt8, 504)...))
end

"BSON element with JavaScript source code."
struct BSONCode
    code::String
end

"`BSON` is a wrapper for C struct `bson_t`."
mutable struct BSON
    handle::Ptr{Cvoid}

    function BSON(handle::Ptr{Cvoid})
        new_bson = new(handle)
        @compat finalizer(destroy!, new_bson)
        return new_bson
    end
end

function destroy!(bson::BSON)
    if bson.handle != C_NULL
        bson_destroy(bson.handle)
        bson.handle = C_NULL
    end
    nothing
end

function Base.deepcopy(bson::BSON) :: BSON
    return BSON(bson_copy(bson.handle))
end

#
# API
#

Base.convert(::Type{String}, oid::BSONObjectId) = bson_oid_to_string(oid)
Base.string(oid::BSONObjectId) = convert(String, oid)
Base.convert(::Type{BSONObjectId}, oid_string::String) = BSONObjectId(oid_string)

Base.convert(::Type{String}, code::BSONCode) = code.code
Base.string(code::BSONCode) = convert(String, code)
Base.convert(::Type{BSONCode}, code_string::String) = BSONCode(code_string)

@static if VERSION < v"0.7-"
    Base.:(==)(c1::BSONCode, c2::BSONCode) = c1.code == c2.code
end

Base.show(io::IO, oid::BSONObjectId) = print(io, "BSONObjectId(\"", string(oid), "\")")
Base.show(io::IO, bson::BSON) = print(io, "BSON(\"", as_json(bson), "\")")
Base.show(io::IO, code::BSONCode) = print(io::IO, "BSONCode(\"$(code.code)\")")

function Base.show(io::IO, err::BSONError)
    for c in err.message
        c_char = Char(c)
        if c_char == '\0'
            break
        else
            print(io, c_char)
        end
    end
end

function get_time(oid::BSONObjectId)
    return Dates.unix2datetime(bson_oid_get_time_t(oid))
end

function BSON()
    handle = bson_new()
    if handle == C_NULL
        error("Failed to create a new BSON document.")
    end
    return BSON(handle)
end

function BSON(json_string::String)
    handle = bson_new_from_json(json_string)
    if handle == C_NULL
        error("Failed parsing JSON to BSON. $json_string.")
    end
    return BSON(handle)
end

function BSON(dict::Dict)
    result = BSON()

    for (k, v) in pairs(dict)
        result[k] = v
    end

    return result
end

function BSON(vector::Vector)
    result = BSON()

    for (i, v) in enumerate(vector)
        result[string(i-1)] = v
    end

    return result
end

"""
    as_json(bson::BSON; canonical::Bool=false) :: String

Converts a `bson` object to a JSON string.

# Example

```julia
julia> document = Mongoc.BSON("{ \"hey\" : 1 }")
BSON("{ "hey" : 1 }")

julia> Mongoc.as_json(document)
"{ \"hey\" : 1 }"

julia> Mongoc.as_json(document, canonical=true)
"{ \"hey\" : { \"\$numberInt\" : \"1\" } }"
```

# C API

* [`bson_as_canonical_extended_json`](http://mongoc.org/libbson/current/bson_as_canonical_extended_json.html)

* [`bson_as_relaxed_extended_json`](http://mongoc.org/libbson/current/bson_as_relaxed_extended_json.html)

"""
function as_json(bson::BSON; canonical::Bool=false) :: String
    cstring = canonical ? bson_as_canonical_extended_json(bson.handle) : bson_as_relaxed_extended_json(bson.handle)
    if cstring == C_NULL
        error("Couldn't convert bson to json.")
    end
    result = unsafe_string(cstring)

    @static if VERSION < v"0.7-"
        bson_free(convert(Ptr{Cvoid}, convert(Ptr{UInt8}, cstring)))
    else
        bson_free(convert(Ptr{Cvoid}, cstring))
    end

    return result
end

#
# Read values from BSON
#

has_field(bson::BSON, key::String) = bson_has_field(bson.handle, key)
Base.haskey(bson::BSON, key::String) = has_field(bson, key)

function bson_iter_init(document::BSON) :: Ref{BSONIter}
    iter_ref = Ref{BSONIter}()
    ok = bson_iter_init(iter_ref, document.handle)
    if !ok
        error("Couldn't create iterator for BSON document.")
    end
    return iter_ref
end

function as_dict(document::BSON)
    iter_ref = bson_iter_init(document)
    return as_dict(iter_ref)
end

function as_dict(iter_ref::Ref{BSONIter})
    result = Dict()
    while bson_iter_next(iter_ref)
        result[unsafe_string(bson_iter_key(iter_ref))] = get_value(iter_ref)
    end
    return result
end

function get_value(iter_ref::Ref{BSONIter})

    local bson_type::BSONType = bson_iter_type(iter_ref)

    if bson_type == BSON_TYPE_UTF8
        return unsafe_string(bson_iter_utf8(iter_ref))
    elseif bson_type == BSON_TYPE_INT64
        return bson_iter_int64(iter_ref)
    elseif bson_type == BSON_TYPE_INT32
        return bson_iter_int32(iter_ref)
    elseif bson_type == BSON_TYPE_DOUBLE
        return bson_iter_double(iter_ref)
    elseif bson_type == BSON_TYPE_OID
        return unsafe_load(bson_iter_oid(iter_ref)) # converts Ptr{BSONObjectId} to BSONObjectId
    elseif bson_type == BSON_TYPE_BOOL
        return bson_iter_bool(iter_ref)
    elseif bson_type == BSON_TYPE_DATE_TIME
        millis = Int64(bson_iter_date_time(iter_ref))
        return isodate2datetime(millis)

    elseif bson_type == BSON_TYPE_ARRAY || bson_type == BSON_TYPE_DOCUMENT

        child_iter_ref = Ref{BSONIter}()
        ok = bson_iter_recurse(iter_ref, child_iter_ref)
        if !ok
            error("Could't iterate array inside BSON.")
        end

        if bson_type == BSON_TYPE_ARRAY
            result_vector = Vector()
            while bson_iter_next(child_iter_ref)
                push!(result_vector, get_value(child_iter_ref))
            end
            return result_vector
        else
            @assert bson_type == BSON_TYPE_DOCUMENT
            return as_dict(child_iter_ref)
        end
    elseif bson_type == BSON_TYPE_CODE
        return BSONCode(unsafe_string(bson_iter_code(iter_ref)))
    else
        error("BSON Type not supported: $bson_type.")
    end
end

function Base.getindex(document::Mongoc.BSON, key::String)
    iter_ref = Ref{BSONIter}()
    ok = bson_iter_init_find(iter_ref, document.handle, key)
    if !ok
        error("Key $key not found.")
    end
    return get_value(iter_ref)
end

#
# Write values to BSON
#

function Base.setindex!(document::BSON, value::BSONObjectId, key::String)
    ok = bson_append_oid(document.handle, key, -1, value)
    if !ok
        error("Couldn't append oid to BSON document.")
    end
    nothing
end

function Base.setindex!(document::BSON, value::Int64, key::String)
    ok = bson_append_int64(document.handle, key, -1, value)
    if !ok
        error("Couldn't append Int64 to BSON document.")
    end
    nothing
end

function Base.setindex!(document::BSON, value::Int32, key::String)
    ok = bson_append_int32(document.handle, key, -1, value)
    if !ok
        error("Couldn't append Int32 to BSON document.")
    end
    nothing
end

function Base.setindex!(document::BSON, value::String, key::String)
    ok = bson_append_utf8(document.handle, key, -1, value, -1)
    if !ok
        error("Couldn't append String to BSON document.")
    end
    nothing
end

function Base.setindex!(document::BSON, value::Bool, key::String)
    ok = bson_append_bool(document.handle, key, -1, value)
    if !ok
        error("Couldn't append Bool to BSON document.")
    end
    nothing
end

function Base.setindex!(document::BSON, value::Float64, key::String)
    ok = bson_append_double(document.handle, key, -1, value)
    if !ok
        error("Couldn't append Float64 to BSON document.")
    end
    nothing
end

function Base.setindex!(document::BSON, value::DateTime, key::String)
    ok = bson_append_date_time(document.handle, key, -1, datetime2isodate(value))
    if !ok
        error("Couldn't append DateTime to BSON document.")
    end
    nothing
end

function Base.setindex!(document::BSON, value::BSON, key::String)
    ok = bson_append_document(document.handle, key, -1, value.handle)
    if !ok
        error("Couldn't append Sub-Document BSON to BSON document.")
    end
    nothing
end

Base.setindex!(document::BSON, value::Dict, key::String) = setindex!(document, BSON(value), key)

function Base.setindex!(document::BSON, value::Vector, key::String)
    sub_document = BSON(value)
    ok = bson_append_array(document.handle, key, -1, sub_document.handle)
    if !ok
        error("Couldn't append array to BSON document.")
    end
    nothing
end

function Base.setindex!(document::BSON, value::BSONCode, key::String)
    ok = bson_append_code(document.handle, key, -1, value.code)
    if !ok
        error("Couldn't append String to BSON document.")
    end
    nothing
end

function Base.setindex!(document::BSON, value::Date, key::String)
    error("BSON format does not support `Date` type. Use `DateTime` instead.")
end


text:
module MultiAssign

export @multiassign

const msg = "@multiassign expression must be in the form x, y, ... = z"

macro multiassign(x)
    #must be doing variable assignment
    @assert isa(x, Expr) msg
    @assert x.head == :(=) msg
    @assert length(x.args) == 2 msg
    #handle single assignment case
    vars = isa(x.args[1], Symbol) ? [x.args[1]] : x.args[1].args
    #get the assignment value
    val = x.args[2]
    #build the multiple assignment expression
    esc(Expr(:block, [:($var = $val) for var ∈ vars]...))
end

end


text:
"""
    drop_axes(x, dims)

Returns all axes of `x` except for those identified by `dims`. Elements of `dims`
must be unique integers or symbols corresponding to the dimensions or names of
dimensions of `x`.

## Examples
```jldoctest
julia> axs = (Index{:a}(1:10), Index{:b}(1:10), Index(1:10));

julia> drop_axes(axs, :a)
(Index{b}(1:10 => Base.OneTo(10)), Index(1:10 => Base.OneTo(10)))

julia> drop_axes(axs, :b)
(Index{a}(1:10 => Base.OneTo(10)), Index(1:10 => Base.OneTo(10)))

julia> to_dims(axs, (:a, :b))
(Index(1:10 => Base.OneTo(10)),)
```
"""
drop_axes(x::AbstractArray, dims) = drop_axes(axes(x), dims)
drop_axes(x::Tuple, dims) = _drop_axes(x, to_dims(x, dims))

_drop_axes(x, dims::Tuple) = __drop_axes(x, dims)
_drop_axes(x, dims) = __drop_axes(x, (dims,))

function __drop_axes(axs::Tuple{Vararg{<:Any,D}}, dims::NTuple{N,Int}) where {D,N}
    for i in 1:N
        1 <= dims[i] <= D || throw(ArgumentError("dropped dims must be in range 1:ndims(A)"))
        for j = 1:i-1
            dims[j] == dims[i] && throw(ArgumentError("dropped dims must be unique"))
        end
    end
    d = ()
    for (i,axisi) in zip(1:D,axs)
        if !in(i, dims)
            d = tuple(d..., axisi)
        end
    end
    return d
end


text:
#
# Mixture Sampler
#
# a mixture of two or more samplers
#

type MixtureSampler <: ModifyingSampler
	selectionsampler::CategoricalSampler
	subsamplers
	function MixtureSampler(subsamplers...)
		length(subsamplers) >= 2 || error("At least two subsamplers are required")
		for sampler in subsamplers
			typeof(sampler) <: Sampler || error("Not all parameters are samplers")
		end
		new(CategoricalSampler(length(subsamplers)), subsamplers)
	end
end

# numparams(s::MixtureSampler) = numparams(s.selectionsampler) + sum(numparams, s.subsamplers)

function paramranges(s::MixtureSampler)
	pr = paramranges(s.selectionsampler)
	for subsampler in s.subsamplers
		pr = [pr; paramranges(subsampler)]
	end
	pr
end

function setparams!(s::MixtureSampler, params)
	nparams = numparams(s)
	length(params) == nparams || error("expected $(nparams) parameters but got $(length(params))")
	paramstart = 1
	paramcount = numparams(s.selectionsampler)
	setparams!(s.selectionsampler, params[paramstart:(paramstart+paramcount-1)])
	for subsampler in s.subsamplers
		paramstart += paramcount
		paramcount = numparams(subsampler)
		setparams!(subsampler, params[paramstart:(paramstart+paramcount-1)])
	end
end

function getparams(s::MixtureSampler)
	ps = getparams(s.selectionsampler)
	for subsampler in s.subsamplers
		ps = [ps; getparams(subsampler)]
	end
	ps
end

function sample(s::MixtureSampler, support, cc::ChoiceContext)
	selectionindex, selectiontrace = sample(s.selectionsampler, (1,length(s.subsamplers)), cc)
	x, trace = sample(s.subsamplers[selectionindex], support, cc)
	x, Dict{Symbol, Any}(:idx=>selectionindex, :sel=>selectiontrace, :sub=>trace)
end

function estimateparams!(s::MixtureSampler, traces)
	estimateparams!(s.selectionsampler, map(trace->trace[:sel], traces))
	subsamplertraces = map(i->[], 1:length(s.subsamplers))
	for trace in traces
		selectionindex = trace[:idx]
		push!(subsamplertraces[selectionindex], trace[:sub])
	end
	for i in 1:length(s.subsamplers)
		estimateparams!(s.subsamplers[i], subsamplertraces[i])
	end
end

amendtrace(s::MixtureSampler, trace, x) = amendtrace(s.subsamplers[trace[:idx]], trace[:sub], x)


text:
#functions
aNb2A = (a,b) -> return(rad2deg(atan(b/a))) #when you know a & b what is A
aNb2B = (a,b) -> return(90-aNb2A(a,b)) #when you know a & b what is B
aNb2c = (a,b) -> return(√(a^2+b^2)) #when you know a & b what is c

aNc2b = (a,c) -> return(√(c^2-a^2)) #when you know a & c what is b
aNc2A = (a,c) -> return(rad2deg(acos(a/c))) #when you know a & c what is A
aNc2B = (a,c) -> return(90-aNc2A(a,c)) #when you know a & c what is B

aNA2b = (a,A) -> return(a*tan(deg2rad(A))) #when you know a & A what is b
aNA2c = (a,A) -> return(a/cos(deg2rad(A))) #when you know a & A what is c
aNA2B = (a,A) -> return(90-A) #when you know a & A what is B

bNA2a = (a,A) -> return(b/tan(deg2rad(A))) #when you know b & A what is a
bNA2c = (a,A) -> return(b/sin(deg2rad(A))) #when you know b & A what is c
bNA2B = (a,A) -> return(90-A) #when you know b & A what is B

bNc2A = (b,c) -> return(rad2deg(asin(b/c))) #when you know b & c what is A
bNc2B = (b,c) -> return(90-bNc2A(b,c)) #when you know b & c what is B
bNc2a = (b,c) -> return(√(c^2-b^2))  #when you know b & c what is a

cNA2a = (c,A) -> return(c*cos(deg2rad(A))) #when you know c & A what is a
cNA2b= (c,A) -> return(c*sin(deg2rad(A))) #when you know c & A what is b
cNA2B = (c,A) -> return(90-A) #when you know c & A what is B
a= 270; b= 175
c =350; A= 30
<<<<<<< Updated upstream
function input(prompt::AbstractString="")
    print(prompt)
    return chomp(readline())
end

response = input("Do you know a & b ? (Y/N)")
println(response)
if(response=="y")
=======

>>>>>>> Stashed changes
    println(aNb2A(a,b))
    println(aNb2B(a,b))
    println(aNb2c(a,b))
    println()
    println(aNc2b(a,c))
    println(aNc2A(a,c))
    println(aNc2B(a,c))
    println()
    println(aNA2b(a,A))
    println(aNA2c(a,A))
    println(aNA2B(a,A))
    println()
    println(bNA2a(a,A))
    println(bNA2c(a,A))
    println(bNA2B(a,A))
    println()
    println(bNc2A(b,c))
    println(bNc2B(b,c))
    println(bNc2a(b,c))
    println()
    println(cNA2a(c,A))
    println(cNA2b(c,A))
    println(cNA2B(c,A))
    println("finished")
<<<<<<< Updated upstream
end
=======
>>>>>>> Stashed changes

max_stars_repo_name: JaewooJoung/julia

max_issues_repo_name: JaewooJoung/julia

max_forks_repo_name: JaewooJoung/julia


text:
# mul! implementations for cases when output/input arrays represent vector elements
function LinearAlgebra.mul!(x_temp::AbstractArray, has_vector::Bool, A::Union{GradientOperator, DivergenceOperator}, M::AbstractArray; overwrite = false)
    for L in A.ops
        mul!(x_temp, has_vector, L, M)
    end
end

function LinearAlgebra.mul!(x_temp::AbstractArray, has_vector::Bool, A::DerivativeOperator{T,N}, M::AbstractArray; overwrite = false) where {T,N}

    # transformation for convenience of computing derivates along particular axis
    if has_vector

        # For Divergence, input would be a physical vector and output a scalar matrix with one less dimension than input
        # Check that axis of differentiation is in the dimensions of M and x_temp
        ndims_M = ndims(M)
        @assert N < ndims_M

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

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

        dims_M = [axes(M)[1:end-1]...]
        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...,N)

        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...), true, A, view(minimally_padded_M, idx...), overwrite = false)
        end
    else
        # For Gradient, input would be a scalar matrix and output a phycial vector with an extra dimension
        # Check that axis of differentiation is in the dimensions of M and x_temp
        ndims_M = ndims(M)
        @assert N <= ndims_M

        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)[1:end-1]...]
        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...,N), false, A, view(minimally_padded_M, idx...), overwrite = false)
        end
    end
end

##################################################################################
# Divergence and Gradient convolutions for general N-dim functions/vectors
##################################################################################

function LinearAlgebra.mul!(x_temp::AbstractVector, has_vector::Bool, A::DerivativeOperator,x::AbstractVector; overwrite = false)
    convolve_BC_left!(x_temp, A, x, overwrite = overwrite)
    convolve_interior!(x_temp, A, x, overwrite = overwrite)
    convolve_BC_right!(x_temp, A, x, overwrite = overwrite)
end

function convolve_interior!(x_temp::AbstractVector{T1},  A::DerivativeOperator{T2}, x::AbstractVector{T1}; overwrite = true) where {T1, T2}

    T = promote_type(T1,T2)
    @assert length(x_temp)+2 == length(x)
    stencil = A.stencil_coefs
    coeff = typeof(A.coefficients)   <: AbstractArray ? A.coefficients : A.coefficients*ones(SVector{length(x_temp)})

    mid = div(A.stencil_length,2)
    if eltype(stencil) <: AbstractVector
        for i in (1+A.boundary_point_count) : (length(x_temp)-A.boundary_point_count)
            xtempi = zero(T)
            cur_stencil = stencil[i-A.boundary_point_count]
            cur_coeff   = coeff[i]
            for idx in 1:A.stencil_length
                xtempi += cur_coeff * cur_stencil[idx] * x[i - mid + idx]
            end
            x_temp[i] = xtempi + !overwrite*x_temp[i]
        end
    else
        @turbo for i in (1+A.boundary_point_count) : (length(x_temp)-A.boundary_point_count)
            xtempi = zero(T)
            cur_coeff   = coeff[i]
            for idx in 1:A.stencil_length
                xtempi += cur_coeff * stencil[idx] * x[i - mid + idx]
            end
            x_temp[i] = xtempi + !overwrite*x_temp[i]
        end
    end
end

function convolve_BC_left!(x_temp::AbstractVector{T1}, A::DerivativeOperator{T2}, x::AbstractVector{T1}; overwrite = true) where {T1,T2}

    T = promote_type(T1,T2)
    stencil = A.low_boundary_coefs
    coeff   = A.coefficients
    for i in 1 : A.boundary_point_count
        cur_stencil = stencil[i]
        cur_coeff   = typeof(coeff)   <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
        xtempi = zero(T)
        for idx in 1:A.boundary_stencil_length
            xtempi += cur_coeff * cur_stencil[idx] * x[idx]
        end
        x_temp[i] = xtempi + !overwrite*x_temp[i]
    end
end

function convolve_BC_right!(x_temp::AbstractVector{T1}, A::DerivativeOperator{T2}, x::AbstractVector{T1}; overwrite = true) where {T1,T2}

    T = promote_type(T1,T2)
    stencil = A.high_boundary_coefs
    coeff   = A.coefficients
    for i in 1 : A.boundary_point_count
        cur_stencil = stencil[i]
        cur_coeff   = typeof(coeff)   <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
        xtempi = zero(T)
        for idx in (A.boundary_stencil_length-1):-1:0
            xtempi += cur_coeff * cur_stencil[end-idx] * x[end-idx]
        end
        x_temp[end-A.boundary_point_count+i] = xtempi + !overwrite*x_temp[end-A.boundary_point_count+i]
    end
end

#########################################################################################
# More efficient Divergence and Gradient convolutions for 2D & 3D functions/vectors
#########################################################################################

for MT in [2,3]
    @eval begin
        function LinearAlgebra.mul!(x_temp::AbstractArray, has_vector::Bool, A::Union{GradientOperator{T,$MT}, DivergenceOperator{T,$MT}}, M::AbstractArray; overwrite = false) where {T}
            convolve_BC_left!(x_temp, A, M, overwrite = overwrite)
            convolve_interior!(x_temp, A, M, overwrite = overwrite)
            convolve_BC_right!(x_temp, A, M, overwrite = overwrite)
        end
    end
end

####### Gradient Convolutions for  2D functions #######
function convolve_interior!(x_temp::AbstractArray{T1}, A::GradientOperator{T2,2}, x::AbstractArray{T3}; overwrite = false) where {T1,T2,T3}

    T = promote_type(T1,T2,T3)
    s = size(x_temp)
    bpc = A.ops[1].boundary_point_count
    mid = div(A.ops[1].stencil_length,2)

    stencil_1 = A.ops[1].stencil_coefs
    stencil_2 = A.ops[2].stencil_coefs

    coeff_1 = typeof(A.ops[1].coefficients)   <: AbstractArray ? A.ops[1].coefficients : A.ops[1].coefficients*ones(SVector{s[1]})
    coeff_2 = typeof(A.ops[2].coefficients)   <: AbstractArray ? A.ops[2].coefficients : A.ops[2].coefficients*ones(SVector{s[2]})

    # Compute derivatives along particular axis and aggregate the outputs

    if eltype(stencil_1) <: AbstractArray
        # Along 1st axis
        for j in 1:s[2], i in 1+bpc:s[1]-bpc
            cur_stencil_1 = stencil_1[i-bpc]
            cur_coeff_1   = coeff_1[i]
            x_temp1 = zero(T)
                for idx in 1:(A.ops[1].stencil_length)
                    x_temp1 += cur_coeff_1 * cur_stencil_1[idx] * x[i+idx-mid,j]
                end
                x_temp[i,j,1] += x_temp1
        end
        # Along 2nd axis
        for j in 1+bpc : s[2]-bpc , i in 1:s[1]
            cur_stencil_2 = stencil_2[j-bpc]
            cur_coeff_2   = coeff_2[j]
            x_temp2 = zero(T)
            for idx in 1:(A.ops[2].stencil_length)
                x_temp2 += cur_coeff_2 * cur_stencil_2[idx] * x[i,j+idx-mid]
            end
            x_temp[i,j,2] += x_temp2
        end
    else
        # Along 1st axis
        @turbo for j in 1:s[2], i in 1+bpc:s[1]-bpc
            cur_coeff_1   = coeff_1[i]
            x_temp1 = zero(T)
                for idx in 1:(A.ops[1].stencil_length)
                    x_temp1 += cur_coeff_1 * stencil_1[idx] * x[i+idx-mid,j]
                end
            x_temp[i,j,1] += x_temp1

        end
        # Along 2nd axis
        @turbo for j in 1+bpc : s[2]-bpc , i in 1:s[1]
            cur_coeff_2   = coeff_2[j]
            x_temp2 = zero(T)
                for idx in 1:(A.ops[2].stencil_length)
                    x_temp2 += cur_coeff_2 * stencil_2[idx] * x[i,j+idx-mid]
                end
                x_temp[i,j,2] += x_temp2
        end
    end
end

function convolve_BC_left!(x_temp::AbstractArray{T1}, A::GradientOperator{T2,2}, x::AbstractArray{T3}; overwrite = false) where {T1,T2,T3}

    T = promote_type(T1,T2,T3)
    s = size(x_temp)
    stencil_1 = A.ops[1].low_boundary_coefs
    stencil_2 = A.ops[2].low_boundary_coefs

    coeff_1 = typeof(A.ops[1].coefficients)   <: AbstractArray ? A.ops[1].coefficients : A.ops[1].coefficients*ones(SVector{s[1]})
    coeff_2 = typeof(A.ops[2].coefficients)   <: AbstractArray ? A.ops[2].coefficients : A.ops[2].coefficients*ones(SVector{s[2]})

    bpc = A.ops[1].boundary_point_count

    # Compute derivatives along particular axis and aggregate the outputs

    # Along 1st axis
    for j in 1:s[2], i in 1:bpc
        cur_stencil_1 = stencil_1[i]
        cur_coeff_1   = coeff_1[i]
        x_temp1 = zero(T)
        for idx in 1:(A.ops[1].stencil_length)
            x_temp1 += cur_coeff_1 * cur_stencil_1[idx] * x[idx,j]
        end
        x_temp[i,j,1] += x_temp1
    end
    # Along 2nd axis
    for j in 1:bpc , i in 1:s[1]
        cur_stencil_2 = stencil_2[j]
        cur_coeff_2   = coeff_2[j]
        x_temp2 = zero(T)
        for idx in 1:(A.ops[2].stencil_length)
            x_temp2 += cur_coeff_2 * cur_stencil_2[idx] * x[i,idx]
        end
        x_temp[i,j,2] += x_temp2
    end
end

function convolve_BC_right!(x_temp::AbstractArray{T1}, A::GradientOperator{T2,2}, x::AbstractArray{T3}; overwrite = false) where {T1,T2,T3}

    T = promote_type(T1,T2,T3)
    s = size(x_temp)
    stencil_1 = A.ops[1].high_boundary_coefs
    stencil_2 = A.ops[2].high_boundary_coefs

    coeff_1 = typeof(A.ops[1].coefficients)   <: AbstractArray ? A.ops[1].coefficients : A.ops[1].coefficients*ones(SVector{s[1]})
    coeff_2 = typeof(A.ops[2].coefficients)   <: AbstractArray ? A.ops[2].coefficients : A.ops[2].coefficients*ones(SVector{s[2]})

    bpc = A.ops[1].boundary_point_count
    bstl = A.ops[1].boundary_stencil_length

    # Compute derivatives along particular axis and aggregate the outputs
    # Along 1st axis
    for j in 1:s[2], i in s[1]-bpc+1 : s[1]
        cur_stencil_1 = stencil_1[i - s[1] + bpc]
        cur_coeff_1   = coeff_1[i]
        x_temp1 = zero(T)
        for idx in 1:bstl
            x_temp1 += cur_coeff_1 * cur_stencil_1[idx] * x[s[1]-bstl+idx+2,j]
        end
        x_temp[i,j,1] += x_temp1
    end
    # Along 2nd axis
    for j in s[2]-bpc+1 : s[2] , i in 1:s[1]
        cur_stencil_2 = stencil_2[j - s[2] + bpc]
        cur_coeff_2   = coeff_2[j]
        x_temp2 = zero(T)
        for idx in 1:bstl
            x_temp2 += cur_coeff_2 * cur_stencil_2[idx] * x[i,s[2]-bstl+idx+2]
        end
        x_temp[i,j,2] += x_temp2
    end
end

####### Divergence Convolutions for  2D Vectors #######
function convolve_interior!(x_temp::AbstractArray{T1}, A::DivergenceOperator{T2,2}, x::AbstractArray{T3}; overwrite = false) where {T1,T2,T3}

    T = promote_type(T1,T2,T3)
    s = size(x_temp)
    bpc = A.ops[1].boundary_point_count
    mid = div(A.ops[1].stencil_length,2)

    stencil_1 = A.ops[1].stencil_coefs
    stencil_2 = A.ops[2].stencil_coefs

    coeff_1 = typeof(A.ops[1].coefficients)   <: AbstractArray ? A.ops[1].coefficients : A.ops[1].coefficients*ones(SVector{s[1]})
    coeff_2 = typeof(A.ops[2].coefficients)   <: AbstractArray ? A.ops[2].coefficients : A.ops[2].coefficients*ones(SVector{s[2]})

    # Compute derivatives along particular axis and aggregate the outputs

    if eltype(stencil_1) <: AbstractArray

        # Along 1st axis
        for j in 1:s[2], i in 1+bpc:s[1]-bpc
            cur_stencil_1 = stencil_1[i-bpc]
            cur_coeff_1   = coeff_1[i]
            x_temp1 = zero(T)
            for idx in 1:(A.ops[1].stencil_length)
                x_temp1 += cur_coeff_1 * cur_stencil_1[idx] * x[i+idx-mid,j,1]
            end
            x_temp[i,j] += x_temp1
        end
        # Along 2nd axis
        for j in 1+bpc : s[2]-bpc , i in 1:s[1]
            cur_stencil_2 = stencil_2[j-bpc]
            cur_coeff_2   = coeff_2[j]
            x_temp2 = zero(T)
            for idx in 1:(A.ops[2].stencil_length)
                x_temp2 += cur_coeff_2 * cur_stencil_2[idx] * x[i,j+idx-mid,2]
            end
            x_temp[i,j] += x_temp2
        end
    else
        # Along 1st axis
        @turbo for j in 1:s[2], i in 1+bpc:s[1]-bpc
            cur_coeff_1   = coeff_1[i]
            x_temp1 = zero(T)
                for idx in 1:(A.ops[1].stencil_length)
                    x_temp1 += cur_coeff_1 * stencil_1[idx] * x[i+idx-mid,j,1]
                end
            x_temp[i,j] += x_temp1

        end
        # Along 2nd axis
        @turbo for j in 1+bpc : s[2]-bpc , i in 1:s[1]
            cur_coeff_2   = coeff_2[j]
            x_temp2 = zero(T)
                for idx in 1:(A.ops[2].stencil_length)
                    x_temp2 += cur_coeff_2 * stencil_2[idx] * x[i,j+idx-mid,2]
                end
                x_temp[i,j] += x_temp2
        end
    end
end

function convolve_BC_left!(x_temp::AbstractArray{T1}, A::DivergenceOperator{T2,2}, x::AbstractArray{T3}; overwrite = false) where {T1,T2,T3}

    T = promote_type(T1,T2,T3)
    s = size(x_temp)
    stencil_1 = A.ops[1].low_boundary_coefs
    stencil_2 = A.ops[2].low_boundary_coefs

    coeff_1 = typeof(A.ops[1].coefficients)   <: AbstractArray ? A.ops[1].coefficients : A.ops[1].coefficients*ones(SVector{s[1]})
    coeff_2 = typeof(A.ops[2].coefficients)   <: AbstractArray ? A.ops[2].coefficients : A.ops[2].coefficients*ones(SVector{s[2]})

    bpc = A.ops[1].boundary_point_count

    # Compute derivatives along particular axis and aggregate the outputs
    # Along 1st axis
    for j in 1:s[2], i in 1:bpc
        cur_stencil_1 = stencil_1[i]
        cur_coeff_1   = coeff_1[i]
        x_temp1 = zero(T)
        for idx in 1:(A.ops[1].stencil_length)
            x_temp1 += cur_coeff_1 * cur_stencil_1[idx] * x[idx,j,1]
        end
        x_temp[i,j] += x_temp1
    end
    # Along 2nd axis
    for j in 1:bpc , i in 1:s[1]
        cur_stencil_2 = stencil_2[j]
        cur_coeff_2   = coeff_2[j]
        x_temp2 = zero(T)
        for idx in 1:(A.ops[2].stencil_length)
            x_temp2 += cur_coeff_2 * cur_stencil_2[idx] * x[i,idx,2]
        end
        x_temp[i,j] += x_temp2
    end
end

function convolve_BC_right!(x_temp::AbstractArray{T1}, A::DivergenceOperator{T2,2}, x::AbstractArray{T3}; overwrite = false) where {T1,T2,T3}

    T = promote_type(T1,T2,T3)
    s = size(x_temp)
    stencil_1 = A.ops[1].high_boundary_coefs
    stencil_2 = A.ops[2].high_boundary_coefs

    coeff_1 = typeof(A.ops[1].coefficients)   <: AbstractArray ? A.ops[1].coefficients : A.ops[1].coefficients*ones(SVector{s[1]})
    coeff_2 = typeof(A.ops[2].coefficients)   <: AbstractArray ? A.ops[2].coefficients : A.ops[2].coefficients*ones(SVector{s[2]})

    bpc = A.ops[1].boundary_point_count
    bstl = A.ops[1].boundary_stencil_length

    # Compute derivatives along particular axis and aggregate the outputs

    # Along 1st axis
    for j in 1:s[2], i in s[1]-bpc+1 : s[1]
        cur_stencil_1 = stencil_1[i - s[1] + bpc]
        cur_coeff_1   = coeff_1[i]
        x_temp1 = zero(T)
        for idx in 1:bstl
            x_temp1 += cur_coeff_1 * cur_stencil_1[idx] * x[s[1]-bstl+idx+2,j,1]
        end
        x_temp[i,j] += x_temp1
    end
    # Along 2nd axis
    for j in s[2]-bpc+1 : s[2] , i in 1:s[1]
        cur_stencil_2 = stencil_2[j - s[2] + bpc]
        cur_coeff_2   = coeff_2[j]
        x_temp2 = zero(T)
        for idx in 1:bstl
            x_temp2 += cur_coeff_2 * cur_stencil_2[idx] * x[i,s[2]-bstl+idx+2,2]
        end
        x_temp[i,j] += x_temp2
    end
end

####### Gradient Convolutions for  3D functions #######
function convolve_interior!(x_temp::AbstractArray{T1}, A::GradientOperator{T2,3}, x::AbstractArray{T3}; overwrite = false) where {T1,T2,T3}

    T = promote_type(T1,T2,T3)
    s = size(x_temp)
    bpc = A.ops[1].boundary_point_count
    mid = div(A.ops[1].stencil_length,2)

    stencil_1 = A.ops[1].stencil_coefs
    stencil_2 = A.ops[2].stencil_coefs
    stencil_3 = A.ops[3].stencil_coefs

    coeff_1 = typeof(A.ops[1].coefficients)   <: AbstractArray ? A.ops[1].coefficients : A.ops[1].coefficients*ones(SVector{s[1]})
    coeff_2 = typeof(A.ops[2].coefficients)   <: AbstractArray ? A.ops[2].coefficients : A.ops[2].coefficients*ones(SVector{s[2]})
    coeff_3 = typeof(A.ops[3].coefficients)   <: AbstractArray ? A.ops[3].coefficients : A.ops[3].coefficients*ones(SVector{s[3]})

    # Compute derivatives along particular axis and aggregate the outputs

    if eltype(stencil_1) <: AbstractArray
        # Along 1st axis
        for k in 1:s[3], j in 1:s[2], i in 1+bpc:s[1]-bpc
            cur_stencil_1 = stencil_1[i-bpc]
            cur_coeff_1   = coeff_1[i]
            x_temp1 = zero(T)
            for idx in 1:(A.ops[1].stencil_length)
                x_temp1 += cur_coeff_1 * cur_stencil_1[idx] * x[i+idx-mid,j,k]
            end
            x_temp[i,j,k,1] += x_temp1
        end
        # Along 2nd axis
        for k in 1:s[3], j in 1+bpc : s[2]-bpc , i in 1:s[1]
            cur_stencil_2 = stencil_2[j-bpc]
            cur_coeff_2   = coeff_2[j]
            x_temp2 = zero(T)
            for idx in 1:(A.ops[2].stencil_length)
                x_temp2 += cur_coeff_2 * cur_stencil_2[idx] * x[i,j+idx-mid,k]
            end
            x_temp[i,j,k,2] += x_temp2
        end
        # Along 3rd axis
        for k in 1+bpc : s[3]-bpc, j in 1:s[2] , i in 1:s[1]
            cur_stencil_3 = stencil_3[k-bpc]
            cur_coeff_3   = coeff_3[k]
            x_temp3 = zero(T)
            for idx in 1:(A.ops[3].stencil_length)
                x_temp3 += cur_coeff_3 * cur_stencil_3[idx] * x[i,j,k+idx-mid]
            end
            x_temp[i,j,k,3] += x_temp3
        end
    else
        # Along 1st axis
        @turbo for k in 1:s[3], j in 1:s[2], i in 1+bpc:s[1]-bpc
            cur_coeff_1   = coeff_1[i]
            x_temp1 = zero(T)
            for idx in 1:(A.ops[1].stencil_length)
                x_temp1 += cur_coeff_1 * stencil_1[idx] * x[i+idx-mid,j,k]
            end
            x_temp[i,j,k,1] += x_temp1
        end
        # Along 2nd axis
        @turbo for k in 1:s[3], j in 1+bpc : s[2]-bpc , i in 1:s[1]
            cur_coeff_2   = coeff_2[j]
            x_temp2 = zero(T)
            for idx in 1:(A.ops[2].stencil_length)
                x_temp2 += cur_coeff_2 * stencil_2[idx] * x[i,j+idx-mid,k]
            end
            x_temp[i,j,k,2] += x_temp2
        end
        # Along 3rd axis
        @turbo for k in 1+bpc : s[3]-bpc, j in 1:s[2] , i in 1:s[1]
            cur_coeff_3   = coeff_3[k]
            x_temp3 = zero(T)
            for idx in 1:(A.ops[3].stencil_length)
                x_temp3 += cur_coeff_3 * stencil_3[idx] * x[i,j,k+idx-mid]
            end
            x_temp[i,j,k,3] += x_temp3
        end
    end
end

function convolve_BC_left!(x_temp::AbstractArray{T1}, A::GradientOperator{T2,3}, x::AbstractArray{T3}; overwrite = false) where {T1,T2,T3}

    T = promote_type(T1,T2,T3)
    s = size(x_temp)
    stencil_1 = A.ops[1].low_boundary_coefs
    stencil_2 = A.ops[2].low_boundary_coefs
    stencil_3 = A.ops[3].low_boundary_coefs

    coeff_1 = typeof(A.ops[1].coefficients)   <: AbstractArray ? A.ops[1].coefficients : A.ops[1].coefficients*ones(SVector{s[1]})
    coeff_2 = typeof(A.ops[2].coefficients)   <: AbstractArray ? A.ops[2].coefficients : A.ops[2].coefficients*ones(SVector{s[2]})
    coeff_3 = typeof(A.ops[3].coefficients)   <: AbstractArray ? A.ops[3].coefficients : A.ops[3].coefficients*ones(SVector{s[3]})

    bpc = A.ops[1].boundary_point_count

    # Compute derivatives along particular axis and aggregate the outputs

    # Along 1st axis
    for k in 1:s[3], j in 1:s[2], i in 1:bpc
        cur_stencil_1 = stencil_1[i]
        cur_coeff_1   = coeff_1[i]
        x_temp1 = zero(T)
        for idx in 1:(A.ops[1].stencil_length)
            x_temp1 += cur_coeff_1 * cur_stencil_1[idx] * x[idx,j,k]
        end
        x_temp[i,j,k,1] += x_temp1
    end
    # Along 2nd axis
    for k in 1:s[3], j in 1:bpc , i in 1:s[1]
        cur_stencil_2 = stencil_2[j]
        cur_coeff_2   = coeff_2[j]
        x_temp2 = zero(T)
        for idx in 1:(A.ops[2].stencil_length)
            x_temp2 += cur_coeff_2 * cur_stencil_2[idx] * x[i,idx,k]
        end
        x_temp[i,j,k,2] += x_temp2
    end
    # Along 3rd axis
    for k in 1:bpc, j in 1:s[2] , i in 1:s[1]
        cur_stencil_3 = stencil_3[k]
        cur_coeff_3   = coeff_3[k]
        x_temp3 = zero(T)
        for idx in 1:(A.ops[3].stencil_length)
            x_temp3 += cur_coeff_3 * cur_stencil_3[idx] * x[i,j,idx]
        end
        x_temp[i,j,k,3] += x_temp3
    end
end

function convolve_BC_right!(x_temp::AbstractArray{T1}, A::GradientOperator{T2,3}, x::AbstractArray{T3}; overwrite = false) where {T1,T2,T3}

    T = promote_type(T1,T2,T3)
    s = size(x_temp)
    stencil_1 = A.ops[1].high_boundary_coefs
    stencil_2 = A.ops[2].high_boundary_coefs
    stencil_3 = A.ops[3].high_boundary_coefs

    coeff_1 = typeof(A.ops[1].coefficients)   <: AbstractArray ? A.ops[1].coefficients : A.ops[1].coefficients*ones(SVector{s[1]})
    coeff_2 = typeof(A.ops[2].coefficients)   <: AbstractArray ? A.ops[2].coefficients : A.ops[2].coefficients*ones(SVector{s[2]})
    coeff_3 = typeof(A.ops[3].coefficients)   <: AbstractArray ? A.ops[3].coefficients : A.ops[3].coefficients*ones(SVector{s[3]})

    bpc = A.ops[1].boundary_point_count
    bstl = A.ops[1].boundary_stencil_length

    # Compute derivatives along particular axis and aggregate the outputs
    # Along 1st axis
    for k in 1:s[3], j in 1:s[2], i in s[1]-bpc+1 : s[1]
        cur_stencil_1 = stencil_1[i - s[1] + bpc]
        cur_coeff_1   = coeff_1[i]
        x_temp1 = zero(T)
        for idx in 1:bstl
            x_temp1 += cur_coeff_1 * cur_stencil_1[idx] * x[s[1]-bstl+idx+2,j,k]
        end
        x_temp[i,j,k,1] += x_temp1
    end
    # Along 2nd axis
    for k in 1:s[3], j in s[2]-bpc+1 : s[2] , i in 1:s[1]
        cur_stencil_2 = stencil_2[j - s[2] + bpc]
        cur_coeff_2   = coeff_2[j]
        x_temp2 = zero(T)
        for idx in 1:bstl
            x_temp2 += cur_coeff_2 * cur_stencil_2[idx] * x[i,s[2]-bstl+idx+2,k]
        end
        x_temp[i,j,k,2] += x_temp2
    end
    # Along 3rd axis
    for k in s[3]-bpc+1 : s[3], j in 1:s[2] , i in 1:s[1]
        cur_stencil_3 = stencil_3[k - s[3] + bpc]
        cur_coeff_3   = coeff_3[k]
        x_temp3 = zero(T)
        for idx in 1:bstl
            x_temp3 += cur_coeff_3 * cur_stencil_3[idx] * x[i,j,s[3]-bstl+idx+2]
        end
        x_temp[i,j,k,3] += x_temp3
    end
end

####### Divergence Convolutions for  3D Vectors #######
function convolve_interior!(x_temp::AbstractArray{T1}, A::DivergenceOperator{T2,3}, x::AbstractArray{T3}; overwrite = false) where {T1,T2,T3}

    T = promote_type(T1,T2,T3)
    s = size(x_temp)
    bpc = A.ops[1].boundary_point_count
    mid = div(A.ops[1].stencil_length,2)

    stencil_1 = A.ops[1].stencil_coefs
    stencil_2 = A.ops[2].stencil_coefs
    stencil_3 = A.ops[3].stencil_coefs

    coeff_1 = typeof(A.ops[1].coefficients)   <: AbstractArray ? A.ops[1].coefficients : A.ops[1].coefficients*ones(SVector{s[1]})
    coeff_2 = typeof(A.ops[2].coefficients)   <: AbstractArray ? A.ops[2].coefficients : A.ops[2].coefficients*ones(SVector{s[2]})
    coeff_3 = typeof(A.ops[3].coefficients)   <: AbstractArray ? A.ops[3].coefficients : A.ops[3].coefficients*ones(SVector{s[3]})

    # Compute derivatives along particular axis and aggregate the outputs

    if eltype(stencil_1) <: AbstractArray
        # Along 1st axis
        for k in 1:s[3], j in 1:s[2], i in 1+bpc:s[1]-bpc
            cur_stencil_1 = stencil_1[i-bpc]
            cur_coeff_1   = coeff_1[i]
            x_temp1 = zero(T)
                for idx in 1:(A.ops[1].stencil_length)
                    x_temp1 += cur_coeff_1 * cur_stencil_1[idx] * x[i+idx-mid,j,k,1]
                end
                x_temp[i,j,k] += x_temp1
        end
        # Along 2nd axis
        for k in 1:s[3], j in 1+bpc : s[2]-bpc , i in 1:s[1]
            cur_stencil_2 = stencil_2[j-bpc]
            cur_coeff_2   = coeff_2[j]
            x_temp2 = zero(T)
            for idx in 1:(A.ops[2].stencil_length)
                x_temp2 += cur_coeff_2 * cur_stencil_2[idx] * x[i,j+idx-mid,k,2]
            end
            x_temp[i,j,k] += x_temp2
        end
        # Along 3rd axis
        for k in 1+bpc : s[3]-bpc, j in 1:s[2] , i in 1:s[1]
            cur_stencil_3 = stencil_3[k-bpc]
            cur_coeff_3   = coeff_3[k]
            x_temp3 = zero(T)
            for idx in 1:(A.ops[3].stencil_length)
                x_temp3 += cur_coeff_3 * cur_stencil_3[idx] * x[i,j,k+idx-mid,3]
            end
            x_temp[i,j,k] += x_temp3
        end
    else
        # Along 1st axis
        @turbo for k in 1:s[3], j in 1:s[2], i in 1+bpc:s[1]-bpc
            cur_coeff_1   = coeff_1[i]
            x_temp1 = zero(T)
                for idx in 1:(A.ops[1].stencil_length)
                    x_temp1 += cur_coeff_1 * stencil_1[idx] * x[i+idx-mid,j,k,1]
                end
                x_temp[i,j,k] += x_temp1
        end
        # Along 2nd axis
        @turbo for k in 1:s[3], j in 1+bpc : s[2]-bpc , i in 1:s[1]
            cur_coeff_2   = coeff_2[j]
            x_temp2 = zero(T)
            for idx in 1:(A.ops[2].stencil_length)
                x_temp2 += cur_coeff_2 * stencil_2[idx] * x[i,j+idx-mid,k,2]
            end
            x_temp[i,j,k] += x_temp2
        end
        # Along 3rd axis
        @turbo for k in 1+bpc : s[3]-bpc, j in 1:s[2] , i in 1:s[1]
            cur_coeff_3   = coeff_3[k]
            x_temp3 = zero(T)
            for idx in 1:(A.ops[3].stencil_length)
                x_temp3 += cur_coeff_3 * stencil_3[idx] * x[i,j,k+idx-mid,3]
            end
            x_temp[i,j,k] += x_temp3
        end
    end
end

function convolve_BC_left!(x_temp::AbstractArray{T1}, A::DivergenceOperator{T2,3}, x::AbstractArray{T3}; overwrite = false) where {T1,T2,T3}

    T = promote_type(T1,T2,T3)
    s = size(x_temp)
    stencil_1 = A.ops[1].low_boundary_coefs
    stencil_2 = A.ops[2].low_boundary_coefs
    stencil_3 = A.ops[3].low_boundary_coefs

    coeff_1 = typeof(A.ops[1].coefficients)   <: AbstractArray ? A.ops[1].coefficients : A.ops[1].coefficients*ones(SVector{s[1]})
    coeff_2 = typeof(A.ops[2].coefficients)   <: AbstractArray ? A.ops[2].coefficients : A.ops[2].coefficients*ones(SVector{s[2]})
    coeff_3 = typeof(A.ops[3].coefficients)   <: AbstractArray ? A.ops[3].coefficients : A.ops[3].coefficients*ones(SVector{s[3]})

    bpc = A.ops[1].boundary_point_count

    # Compute derivatives along particular axis and aggregate the outputs

    # Along 1st axis
    for k in 1:s[3], j in 1:s[2], i in 1:bpc
        cur_stencil_1 = stencil_1[i]
        cur_coeff_1   = coeff_1[i]
        x_temp1 = zero(T)
        for idx in 1:(A.ops[1].stencil_length)
            x_temp1 += cur_coeff_1 * cur_stencil_1[idx] * x[idx,j,k,1]
        end
        x_temp[i,j,k] += x_temp1
    end
    # Along 2nd axis
    for k in 1:s[3], j in 1:bpc , i in 1:s[1]
        cur_stencil_2 = stencil_2[j]
        cur_coeff_2   = coeff_2[j]
        x_temp2 = zero(T)
        for idx in 1:(A.ops[2].stencil_length)
            x_temp2 += cur_coeff_2 * cur_stencil_2[idx] * x[i,idx,k,2]
        end
        x_temp[i,j,k] += x_temp2
    end
    # Along 3rd axis
    for k in 1:bpc, j in 1:s[2] , i in 1:s[1]
        cur_stencil_3 = stencil_3[k]
        cur_coeff_3   = coeff_3[k]
        x_temp3 = zero(T)
        for idx in 1:(A.ops[3].stencil_length)
            x_temp3 += cur_coeff_3 * cur_stencil_3[idx] * x[i,j,idx,3]
        end
        x_temp[i,j,k] += x_temp3
    end
end

function convolve_BC_right!(x_temp::AbstractArray{T1}, A::DivergenceOperator{T2,3}, x::AbstractArray{T3}; overwrite = false) where {T1,T2,T3}

    T = promote_type(T1,T2,T3)
    s = size(x_temp)
    stencil_1 = A.ops[1].high_boundary_coefs
    stencil_2 = A.ops[2].high_boundary_coefs
    stencil_3 = A.ops[3].high_boundary_coefs

    coeff_1 = typeof(A.ops[1].coefficients)   <: AbstractArray ? A.ops[1].coefficients : A.ops[1].coefficients*ones(SVector{s[1]})
    coeff_2 = typeof(A.ops[2].coefficients)   <: AbstractArray ? A.ops[2].coefficients : A.ops[2].coefficients*ones(SVector{s[2]})
    coeff_3 = typeof(A.ops[3].coefficients)   <: AbstractArray ? A.ops[3].coefficients : A.ops[3].coefficients*ones(SVector{s[3]})

    bpc = A.ops[1].boundary_point_count
    bstl = A.ops[1].boundary_stencil_length

    # Compute derivatives along particular axis and aggregate the outputs
    # Along 1st axis
    for k in 1:s[3], j in 1:s[2], i in s[1]-bpc+1 : s[1]
        cur_stencil_1 = stencil_1[i - s[1] + bpc]
        cur_coeff_1   = coeff_1[i]
        x_temp1 = zero(T)
        for idx in 1:bstl
            x_temp1 += cur_coeff_1 * cur_stencil_1[idx] * x[s[1]-bstl+idx+2,j,k,1]
        end
        x_temp[i,j,k] += x_temp1
    end
    # Along 2nd axis
    for k in 1:s[3], j in s[2]-bpc+1 : s[2] , i in 1:s[1]
        cur_stencil_2 = stencil_2[j - s[2] + bpc]
        cur_coeff_2   = coeff_2[j]
        x_temp2 = zero(T)
        for idx in 1:bstl
            x_temp2 += cur_coeff_2 * cur_stencil_2[idx] * x[i,s[2]-bstl+idx+2,k,2]
        end
        x_temp[i,j,k] += x_temp2
    end
    # Along 3rd axis
    for k in s[3]-bpc+1 : s[3], j in 1:s[2] , i in 1:s[1]
        cur_stencil_3 = stencil_3[k - s[3] + bpc]
        cur_coeff_3   = coeff_3[k]
        x_temp3 = zero(T)
        for idx in 1:bstl
            x_temp3 += cur_coeff_3 * cur_stencil_3[idx] * x[i,j,s[3]-bstl+idx+2,3]
        end
        x_temp[i,j,k] += x_temp3
    end
end

##################################################################################
# Curl convolutions
##################################################################################

# Against a standard vector, assume already padded and just apply the stencil

function LinearAlgebra.mul!(x_temp::AbstractArray{T,4}, A::CurlOperator, x::AbstractArray{T,4}; overwrite = false) where {T<:Number}
    convolve_BC_left!(x_temp, x, A, overwrite = overwrite)
    convolve_interior!(x_temp, x, A, overwrite = overwrite)
    convolve_BC_right!(x_temp, x, A, overwrite = overwrite)
end

function convolve_interior!(x_temp::AbstractArray{T,4}, u::AbstractArray{T,4}, A::CurlOperator; overwrite = false) where {T<:Number}

    s = size(x_temp)
    stencil_1 = A.ops[1].stencil_coefs
    stencil_2 = A.ops[2].stencil_coefs
    stencil_3 = A.ops[3].stencil_coefs

    coeff_1 = typeof(A.ops[1].coefficients)   <: AbstractArray ? A.ops[1].coefficients : A.ops[1].coefficients*ones(SVector{s[1]})
    coeff_2 = typeof(A.ops[2].coefficients)   <: AbstractArray ? A.ops[2].coefficients : A.ops[2].coefficients*ones(SVector{s[2]})
    coeff_3 = typeof(A.ops[3].coefficients)   <: AbstractArray ? A.ops[3].coefficients : A.ops[3].coefficients*ones(SVector{s[3]})

    bpc = A.ops[1].boundary_point_count
    mid = div(A.ops[1].stencil_length,2)

    # Compute derivatives in particular dimensions and aggregate the outputs

    if eltype(stencil_1) <: AbstractArray
        # Along 1st axis
        for k in 1:s[3], j in 1:s[2], i in 1+bpc:s[1]-bpc
            cur_stencil_1 = stencil_1[i-bpc]
            cur_coeff_1   = coeff_1[i]
            x_temp2 = zero(T)
            x_temp3 = zero(T)
            for idx in 1:(A.ops[1].stencil_length)
                x_temp2 += -cur_coeff_1*cur_stencil_1[idx]*u[i+idx-mid,j,k,3]
                x_temp3 += cur_coeff_1*cur_stencil_1[idx]*u[i+idx-mid,j,k,2]
            end
            x_temp[i,j,k,3] += x_temp3
            x_temp[i,j,k,2] += x_temp2
        end

        # Along 2nd axis
        for k in 1:s[3], j in 1+bpc : s[2]-bpc , i in 1:s[1]
            cur_stencil_2 = stencil_2[j-bpc]
            cur_coeff_2   = coeff_2[j]
            x_temp1 = zero(T)
            x_temp3 = zero(T)
            for idx in 1:(A.ops[2].stencil_length)
                x_temp1 += cur_coeff_2*cur_stencil_2[idx]*u[i,j+idx-mid,k,3]
                x_temp3 += -cur_coeff_2*cur_stencil_2[idx]*u[i,j+idx-mid,k,1]
            end
            x_temp[i,j,k,3] += x_temp3
            x_temp[i,j,k,1] += x_temp1
        end

        # Along 3rd axis
        for k in 1+bpc : s[3]-bpc, j in 1 : s[2] , i in 1 : s[1]
            cur_stencil_3 = stencil_3[k-bpc]
            cur_coeff_3   = coeff_3[k]
            x_temp2 = zero(T)
            x_temp1 = zero(T)
            for idx in 1:(A.ops[3].stencil_length)
                x_temp1 += -cur_coeff_3*cur_stencil_3[idx]*u[i,j,k+idx-mid,2]
                x_temp2 += cur_coeff_3*cur_stencil_3[idx]*u[i,j,k+idx-mid,1]
            end
            x_temp[i,j,k,1] += x_temp1
            x_temp[i,j,k,2] += x_temp2
        end
    else
        # Along 1st axis
        @turbo for k in 1:s[3], j in 1:s[2], i in 1+bpc:s[1]-bpc
            cur_coeff_1   = coeff_1[i]
            x_temp2 = zero(T)
            x_temp3 = zero(T)
            for idx in 1:(A.ops[1].stencil_length)
                x_temp2 += -cur_coeff_1*stencil_1[idx]*u[i+idx-mid,j,k,3]
                x_temp3 += cur_coeff_1*stencil_1[idx]*u[i+idx-mid,j,k,2]
            end
            x_temp[i,j,k,3] += x_temp3
            x_temp[i,j,k,2] += x_temp2
        end

        # Along 2nd axis
        @turbo for k in 1:s[3], j in 1+bpc : s[2]-bpc , i in 1:s[1]
            cur_coeff_2   = coeff_2[j]
            x_temp1 = zero(T)
            x_temp3 = zero(T)
            for idx in 1:(A.ops[2].stencil_length)
                x_temp1 += cur_coeff_2*stencil_2[idx]*u[i,j+idx-mid,k,3]
                x_temp3 += -cur_coeff_2*stencil_2[idx]*u[i,j+idx-mid,k,1]
            end
            x_temp[i,j,k,3] += x_temp3
            x_temp[i,j,k,1] += x_temp1
        end

        # Along 3rd axis
        @turbo for k in 1+bpc : s[3]-bpc, j in 1 : s[2] , i in 1 : s[1]
            cur_coeff_3   = coeff_3[k]
            x_temp2 = zero(T)
            x_temp1 = zero(T)
            for idx in 1:(A.ops[3].stencil_length)
                x_temp1 += -cur_coeff_3*stencil_3[idx]*u[i,j,k+idx-mid,2]
                x_temp2 += cur_coeff_3*stencil_3[idx]*u[i,j,k+idx-mid,1]
            end
            x_temp[i,j,k,1] += x_temp1
            x_temp[i,j,k,2] += x_temp2
        end
    end
end

function convolve_BC_left!(x_temp::AbstractArray{T,4}, u::AbstractArray{T,4}, A::CurlOperator; overwrite = false) where {T<:Number}

    s = size(x_temp)
    stencil_1 = A.ops[1].low_boundary_coefs
    stencil_2 = A.ops[2].low_boundary_coefs
    stencil_3 = A.ops[3].low_boundary_coefs

    coeff_1 = typeof(A.ops[1].coefficients)   <: AbstractArray ? A.ops[1].coefficients : A.ops[1].coefficients*ones(SVector{s[1]})
    coeff_2 = typeof(A.ops[2].coefficients)   <: AbstractArray ? A.ops[2].coefficients : A.ops[2].coefficients*ones(SVector{s[2]})
    coeff_3 = typeof(A.ops[3].coefficients)   <: AbstractArray ? A.ops[3].coefficients : A.ops[3].coefficients*ones(SVector{s[3]})

    bpc = A.ops[1].boundary_point_count

    # Compute derivatives in particular dimensions and aggregate the outputs
    for  k in 1:s[3], j in 1:s[2] , i in 1 : bpc
        cur_stencil_1 = stencil_1[i]
        cur_coeff_1   = coeff_1[i]
        x_temp2 = zero(T)
        x_temp3 = zero(T)
        for idx in 1:(A.ops[1].stencil_length)
            x_temp2 += -cur_coeff_1*cur_stencil_1[idx]*u[idx,j,k,3]
            x_temp3 += cur_coeff_1*cur_stencil_1[idx]*u[idx,j,k,2]
        end
        x_temp[i,j,k,3] += x_temp3
        x_temp[i,j,k,2] += x_temp2
    end

    for k in 1:s[3], j in 1:bpc , i in 1:s[1]
        cur_stencil_2 = stencil_2[j]
        cur_coeff_2   = coeff_2[j]
        x_temp1 = zero(T)
        x_temp3 = zero(T)
        for idx in 1:(A.ops[2].stencil_length)
            x_temp1 += cur_coeff_2*cur_stencil_2[idx]*u[i,idx,k,3]
            x_temp3 += -cur_coeff_2*cur_stencil_2[idx]*u[i,idx,k,1]
        end
        x_temp[i,j,k,1] += x_temp1
        x_temp[i,j,k,3] += x_temp3
    end

    for k in 1:bpc, j in 1:s[2] , i in 1:s[1]
        cur_stencil_3 = stencil_3[k]
        cur_coeff_3   = coeff_3[k]
        x_temp2 = zero(T)
        x_temp1 = zero(T)
        for idx in 1:(A.ops[1].stencil_length)
            x_temp1 += -cur_coeff_3*cur_stencil_3[idx]*u[i,j,idx,2]
            x_temp2 += cur_coeff_3*cur_stencil_3[idx]*u[i,j,idx,1]
        end
        x_temp[i,j,k,1] += x_temp1
        x_temp[i,j,k,2] += x_temp2
    end
end

function convolve_BC_right!(x_temp::AbstractArray{T,4}, u::AbstractArray{T,4}, A::CurlOperator; overwrite = false) where {T<:Number}

    s = size(x_temp)
    stencil_1 = A.ops[1].high_boundary_coefs
    stencil_2 = A.ops[2].high_boundary_coefs
    stencil_3 = A.ops[3].high_boundary_coefs

    coeff_1 = typeof(A.ops[1].coefficients)   <: AbstractArray ? A.ops[1].coefficients : A.ops[1].coefficients*ones(SVector{s[1]})
    coeff_2 = typeof(A.ops[2].coefficients)   <: AbstractArray ? A.ops[2].coefficients : A.ops[2].coefficients*ones(SVector{s[2]})
    coeff_3 = typeof(A.ops[3].coefficients)   <: AbstractArray ? A.ops[3].coefficients : A.ops[3].coefficients*ones(SVector{s[3]})

    bpc = A.ops[1].boundary_point_count
    bstl = A.ops[1].boundary_stencil_length

    # Compute derivatives in particular dimensions and aggregate the outputs
    for k in 1 : s[3], j in 1 : s[2] ,i in s[1]-bpc+1 : s[1]
        cur_stencil_1 = stencil_1[i - s[1] + bpc]
        cur_coeff_1   = coeff_1[i]
        x_temp2 = zero(T)
        x_temp3 = zero(T)
        for idx in 1:bstl
            x_temp2 += -cur_coeff_1*cur_stencil_1[idx]*u[s[1]-bstl+idx+2,j,k,3]
            x_temp3 += cur_coeff_1*cur_stencil_1[idx]*u[s[1]-bstl+idx+2,j,k,2]
        end
        x_temp[i,j,k,3] += x_temp3
        x_temp[i,j,k,2] += x_temp2
    end

    for k in 1 : s[3], j in s[2]-bpc+1 : s[2] , i in 1 : s[1]
        cur_stencil_2 = stencil_2[j - s[2] + bpc]
        cur_coeff_2   = coeff_2[j]
        x_temp1 = zero(T)
        x_temp3 = zero(T)
        for idx in 1:bstl
            x_temp1 += cur_coeff_2*cur_stencil_2[idx]*u[i,s[2]-bstl+idx+2,k,3]
            x_temp3 += -cur_coeff_2*cur_stencil_2[idx]*u[i,s[2]-bstl+idx+2,k,1]
        end
        x_temp[i,j,k,1] += x_temp1
        x_temp[i,j,k,3] += x_temp3
    end

    for k in s[3]-bpc+1 : s[3], j in 1 : s[2] , i in 1 : s[1]
        cur_stencil_3 = stencil_3[k - s[3] + bpc]
        cur_coeff_3   = coeff_3[k]
        x_temp1 = zero(T)
        x_temp2 = zero(T)
        for idx in 1:bstl
            x_temp1 += -cur_coeff_3*cur_stencil_3[idx]*u[i,j,s[3]-bstl+idx+2,2]
            x_temp2 += cur_coeff_3*cur_stencil_3[idx]*u[i,j,s[3]-bstl+idx+2,1]
        end
        x_temp[i,j,k,1] += x_temp1
        x_temp[i,j,k,2] += x_temp2
    end
end


text:
using NeuralSimulation

# Hodgkin-Huxley channel gating functions
αn(V) = V == -55.0 ? 0.1 : (0.01*(V + 55.0))/(1.0 - exp(-(V + 55.0)/10.0))
βn(V) = 0.125 * exp(-(V + 65.0)/80.0)
αm(V) = V == -40.0 ? 1.0 : (0.1*(V + 40.0))/(1.0 - exp(-(V + 40.0)/10.0))
βm(V) = 4.0*exp(-(V + 65.0)/18.0)
αh(V) = 0.07*exp(-(V+65.0)/20.0)
βh(V) = 1.0/(1.0 + exp(-(V + 35.0)/10.0))

# Steady-state channel gating
m∞(V) = αm(V)/(αm(V) + βm(V))
h∞(V) = αh(V)/(αh(V) + βh(V))
n∞(V) = αn(V)/(αn(V) + βn(V))

# Updating gates (rewritten as explicit solution to backwards euler)
@inline function BEgates!(u, p, dt)
    V = view(u, :, 1)
    m = view(u, :, 2)
    h = view(u, :, 3)
    n = view(u, :, 4)
    @inbounds for i in 1:length(V)
        m[i] = (m[i] + αm(V[i])*dt)/(1 + (αm(V[i]) + βm(V[i]))*dt)
        h[i] = (h[i] + αh(V[i])*dt)/(1 + (αh(V[i]) + βh(V[i]))*dt)
        n[i] = (n[i] + αn(V[i])*dt)/(1 + (αn(V[i]) + βn(V[i]))*dt)
    end
end

# In-place update of the Hines coefficients + RHS
@inline function f_hines!(Hines::HinesMatrix{S, BackwardsEuler}, u, p, dt) where {S}
    V = view(u, :, 1)
    m = view(u, :, 2)
    h = view(u, :, 3)
    n = view(u, :, 4)
    a, rl, cm, gna, gk, gl, Ena, Ek, El, Iapp, dx = p
    @inbounds for k in 1:length(V)
        Hines.D[k] = 1 + (a*dt)/(rl*cm*dx^2) + (dt/cm)*(gna*m[k]^3*h[k] + gk*n[k]^4 + gl)
        V[k] +=  (dt/cm)*(gna*m[k]^3*h[k]*Ena + gk*n[k]^4*Ek + gl*El)
    end
    # current injected to first compartment; for testing only!
    V[1] += (Iapp*dt)/(pi*a*cm*dx)
    return nothing
end

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

dt = 0.025 # ms
Iapp = 1.0e-6 # mA (1 nA)
cable_length = 2000.0e-3 # 2000 µm length
N = 2000
dx = cable_length/N
radius = 2.0e-3 # µm in mm
area = 2*pi*radius*dx

rl = 1.0e6 # 1 kΩ in milliohms*mm
gl  = 0.003 #mS/mm^2
gna = 1.2 # mS/mm^2 (modified from 1.2 mS)
gk  = .36 # mS/mm^2
El  = -54.4 # mV
Ena = 50.0 # mV
Ek  = -77.0 # mV
cm  = 1.0e-5 # mF/mm^2

Vinit = -55.0 # mV
tspan = (0.0, 200.0)

p = [radius, rl, cm, gna, gk, gl, Ena, Ek, El, Iapp, dx]

cable = CableSection(radius, cable_length, rl, cm, N)

u0 = zeros(2000,4)
u0[:,1] .= Vinit
u0[:,2] .= m∞(Vinit)
u0[:,3] .= h∞(Vinit)
u0[:,4] .= n∞(Vinit)

fn = HHCableFunction{true}(f_hines!, BEgates!)
prob = HHCableProblem(fn, cable, u0, tspan, p)
@time sol = solve(prob, BackwardsEuler(), dt = 0.025)

max_stars_repo_name: wsphillips/FutureNeuralSimulation

max_issues_repo_name: wsphillips/FutureNeuralSimulation

max_forks_repo_name: wsphillips/FutureNeuralSimulation