text:
module TestNAArrayOperators

using FactCheck
using DataCubes
using DataCubes: AbstractArrayWrapper,FloatNAArray,wrap_array

facts("NAArrayOperators tests") do
  context("AbstractArrayWrapper tests") do
    @fact broadcast(+, 1, AbstractArrayWrapper(nalift([1,2,3]))) --> AbstractArrayWrapper(nalift([2,3,4]))
    @fact broadcast(-, AbstractArrayWrapper(nalift([1,2,3])), 1) --> AbstractArrayWrapper(nalift([0,1,2]))
    @fact broadcast(*, AbstractArrayWrapper(nalift([1,2,3])),
            AbstractArrayWrapper(nalift([5.0,6.0,7.0]))) --> AbstractArrayWrapper(nalift([5.0,12.0,21.0]))
    @fact map(x->Nullable(x.value*2), AbstractArrayWrapper(nalift([1,2,3]))) --> AbstractArrayWrapper(nalift([2,4,6]))
    @fact (x=AbstractArrayWrapper([5,1,3,2,4]);sort!(x);x) --> AbstractArrayWrapper([1,2,3,4,5])
    @fact sort(AbstractArrayWrapper([5 1 3 2 4]),2) --> AbstractArrayWrapper([1 2 3 4 5])
    @fact sort(AbstractArrayWrapper([5,1,3,2,4])) --> AbstractArrayWrapper([1,2,3,4,5])
    @fact map((x,y)->DataCubes.naop_plus(x,y), nalift([1,2,3,4,5]), @nalift([1,2,NA,4,5])) --> @nalift([2,4,NA,8,10])
    @fact broadcast(>, @nalift([1,2,NA]), 1) --> @nalift([false, true, NA])
    @fact broadcast(<, 1, @nalift([1,2,NA])) --> @nalift([false, true, NA])
    @fact broadcast(<, @nalift([1,2,3,NA]), @nalift([NA,3,1,2])) --> @nalift([NA,true,false,NA])
    @fact broadcast(-, @nalift([1.0,2.0,NA]), 1) --> @nalift([0.0,1.0,NA])
    @fact broadcast(-, 1, @nalift([1.0,2.0,NA])) --> @nalift([0.0,-1.0,NA])
    @fact broadcast(*, @nalift([1.0,2.0,NA]), @nalift([NA,2.0,3.0])) --> @nalift([NA,4.0,NA])
    @fact -@nalift([1.0,2.0,NA]) --> @nalift([-1.0,-2.0,NA])
    @fact ~@nalift([1,2,NA]) --> @nalift([~1,~2,NA])
    @fact broadcast(>, AbstractArrayWrapper([1,2,3]), 1) --> AbstractArrayWrapper([false,true,true])
    @fact broadcast(>, AbstractArrayWrapper([1,2,3]), 1.5) --> AbstractArrayWrapper([false,true,true])
    @fact broadcast(<, 1, AbstractArrayWrapper([1,2,3])) --> AbstractArrayWrapper([false,true,true])
    @fact broadcast(>, AbstractArrayWrapper([1,2,3]), Nullable(1)) --> nalift([false,true,true])
    @fact AbstractArrayWrapper([1,2,3]) + 1.0 --> AbstractArrayWrapper([2.0,3.0,4.0])
    @fact AbstractArrayWrapper([1,2,3]) + Nullable(1.0) --> nalift([2.0,3.0,4.0])
    @fact broadcast(<, Nullable(1), AbstractArrayWrapper([1,2,3])) --> nalift([false,true,true])
    @fact broadcast(<, AbstractArrayWrapper([3.0,2.0,1.0]), AbstractArrayWrapper([1.0,2.0,3.0])) --> AbstractArrayWrapper([false,false,true])
    @fact broadcast(<, AbstractArrayWrapper([3,2,1]), nalift([1,2,3])) --> nalift([false,false,true])
    @fact broadcast(<, @nalift([3,NA,1]), AbstractArrayWrapper([1.0,2.0,3.0])) --> @nalift([false,NA,true])
    @fact broadcast(<, nalift([3,2,1]), nalift([1.0,2.0,3.0])) --> nalift([false,false,true])
    @fact @nalift([NA,2,1]) + @nalift([1.0,NA,3.0]) --> @nalift([NA,NA,4.0])
    @fact nalift([3,2,1]) + AbstractArrayWrapper([1.0,2.0,3.0]) --> nalift([4.0,4.0,4.0])
    @fact AbstractArrayWrapper([3,2,1]) + AbstractArrayWrapper([1.0,2.0,3.0]) --> AbstractArrayWrapper([4.0,4.0,4.0])
    @fact AbstractArrayWrapper([3,2,1]) + 5.0 --> AbstractArrayWrapper([8.0,7.0,6.0])
    @fact AbstractArrayWrapper([3,2,1]) + Nullable(5.0) --> nalift([8.0,7.0,6.0])
    @fact -AbstractArrayWrapper([3,2,1]) --> AbstractArrayWrapper([-3,-2,-1])
    @fact -AbstractArrayWrapper([3.0,2.0,1.0]) --> AbstractArrayWrapper([-3.0,-2.0,-1.0])
    @fact -nalift([3.0,2.0,1.0]) --> nalift([-3.0,-2.0,-1.0])
    @fact -darr(a=[3.0,2.0,1.0]) --> darr(a=[-3.0,-2.0,-1.0])
    @fact -larr(a=[3.0,2.0,1.0], axis1=[:x,:y,:z]) --> larr(a=[-3.0,-2.0,-1.0], axis1=[:x,:y,:z])
    @fact nalift([1.0,2.0,3.0]) + nalift([1,2,3]) --> nalift([2.0,4.0,6.0])
    @fact nalift([1.0,2.0,3.0]) + nalift([1.0,2.0,3.0]) --> nalift([2.0,4.0,6.0])
    @fact AbstractArrayWrapper(map(Nullable, [1.0,2.0,3.0])) + nalift([1.0,2.0,3.0]) --> nalift([2.0,4.0,6.0])
    @fact AbstractArrayWrapper(map(Nullable, [1.0,2.0,3.0])) + @nalift([1.0,2.0,NA]) --> @nalift([2.0,4.0,NA])
    @fact @nalift([1.0,2.0,NA]) + AbstractArrayWrapper([Nullable(1.0),Nullable{Float64}(), Nullable(3.0)]) --> @nalift([2.0,NA,NA])
    @fact AbstractArrayWrapper(map(Nullable, [1.0,2.0,3.0])) + nalift([1,2,3]) --> nalift([2.0,4.0,6.0])
    @fact AbstractArrayWrapper([Nullable(1.0),Nullable(2.0),Nullable{Float64}()]) + AbstractArrayWrapper([1,2,3]) --> @nalift([2.0,4.0,NA])
    @fact AbstractArrayWrapper([Nullable(1), Nullable(2), Nullable{Int}()]) + AbstractArrayWrapper([1.0,2.0,3.0]) --> @nalift([2.0,4.0,NA])
    @fact AbstractArrayWrapper([1.0,2.0,3.0]) + AbstractArrayWrapper([Nullable(1), Nullable(2), Nullable{Int}()]) --> @nalift([2.0,4.0,NA])
    @fact AbstractArrayWrapper([1.0,2.0,3.0]) + AbstractArrayWrapper([1,2,3]) --> [2.0,4.0,6.0]
    @fact broadcast(/, AbstractArrayWrapper([1,2,3]), AbstractArrayWrapper([1,2,3])) --> [1.0,1.0,1.0]
    @fact wrap_array([true,true,false]) & wrap_array([false,true,true]) --> wrap_array([false,true,false])
    @fact wrap_array([true,true,false]) & true --> wrap_array([true,true,false])
    @fact wrap_array([true,true,false]) & false --> wrap_array([false,false,false])
    @fact (2+3im) * @nalift([1,2,NA]) --> @nalift([2+3im,4+6im,NA])
    @fact @nalift([1,2,NA]) * (2+3im)--> @nalift([2+3im,4+6im,NA])
    @fact Nullable(2+3im) * @nalift([1,2,NA]) --> @nalift([2+3im,4+6im,NA])
    @fact @nalift([1,2,NA]) * Nullable(2+3im)--> @nalift([2+3im,4+6im,NA])
    @fact Nullable(im) / @nalift([1,2,NA]) --> @nalift([im,0.5im,NA])
    @fact @nalift([1,2,NA]) / Nullable(im)--> @nalift([-im,-2im,NA])
    @fact Nullable(e) / @nalift([1,2,NA]) --> @nalift([convert(Float64,e),convert(Float64,e/2),NA])
    @fact @nalift([1,2,NA]) / Nullable(e)--> @nalift([convert(Float64,1/e),convert(Float64,2/e),NA])
    @fact e / @nalift([1,2,NA]) --> @nalift([convert(Float64,e),convert(Float64,e/2),NA])
    @fact @nalift([1,2,NA]) / e--> @nalift([convert(Float64,1/e),convert(Float64,2/e),NA])
    @fact nalift([1.0,2.0,3.0]) + 1.0 --> nalift([2.0,3.0,4.0])
    @fact 1.0 + nalift([1.0,2.0,3.0]) --> nalift([2.0,3.0,4.0])
    @fact nalift([1.0,2.0,3.0]) + 1 --> nalift([2.0,3.0,4.0])
    @fact 1 + nalift([1.0,2.0,3.0]) --> nalift([2.0,3.0,4.0])
    @fact 1.0 + nalift([1,2,3]) --> nalift([2.0,3.0,4.0])
    @fact nalift([1,2,3]) + 1.0 --> nalift([2.0,3.0,4.0])
    @fact nalift([1.0,2.0,3.0]) + Nullable(1.0) --> nalift([2.0,3.0,4.0])
    @fact Nullable(1.0) + nalift([1.0,2.0,3.0]) --> nalift([2.0,3.0,4.0])
    @fact nalift([1.0,2.0,3.0]) + Nullable(1) --> nalift([2.0,3.0,4.0])
    @fact Nullable(1) + nalift([1.0,2.0,3.0]) --> nalift([2.0,3.0,4.0])
    @fact Nullable(1.0) + nalift([1,2,3]) --> nalift([2.0,3.0,4.0])
    @fact nalift([1,2,3]) + Nullable(1.0) --> nalift([2.0,3.0,4.0])
    @fact nalift([1.0,2.0,3.0]) * 1.0 --> nalift([2.0,3.0,4.0]) - 1
    @fact 1.0 * nalift([1.0,2.0,3.0]) --> nalift([2.0,3.0,4.0]) - 1
    @fact nalift([1.0,2.0,3.0]) * 1 --> nalift([2.0,3.0,4.0]) - 1
    @fact 1 * nalift([1.0,2.0,3.0]) --> nalift([2.0,3.0,4.0]) - 1
    @fact 1.0 * nalift([1,2,3]) --> nalift([2.0,3.0,4.0]) - 1
    @fact nalift([1,2,3]) * 1.0 --> nalift([2.0,3.0,4.0]) - 1
    @fact nalift([1.0,2.0,3.0]) * Nullable(1.0) --> nalift([2.0,3.0,4.0]) - 1
    @fact Nullable(1.0) * nalift([1.0,2.0,3.0]) --> nalift([2.0,3.0,4.0]) - 1
    @fact @nalift([1.0,2.0,NA]) * Nullable(1) --> @nalift([2.0,3.0,NA]) - 1
    @fact Nullable(1) * nalift([1.0,2.0,3.0]) --> nalift([2.0,3.0,4.0]) - 1
    @fact Nullable(1.0) * @nalift([NA,2,3]) --> @nalift([NA,3.0,4.0]) - 1
    @fact nalift([1,2,3]) * Nullable(1.0) --> nalift([2.0,3.0,4.0]) - 1
    @fact nalift([1,2,3]) * Nullable{Int}() --> nalift([Nullable{Int}(), Nullable{Int}(), Nullable{Int}()])
    @fact broadcast(/, @nalift([1,2,3,NA,NA,5]), @nalift([2,1,2,3,NA,NA])) --> @nalift([0.5,2.0,1.5,NA,NA,NA])
    @fact broadcast(/, @nalift([1,2,3,NA,NA,5]), @nalift([2,1,2,3,NA,NA])) --> @nalift([0.5,2.0,1.5,NA,NA,NA])
    @fact broadcast(/, @nalift([1,2,3,NA,NA,5]), DataCubes.wrap_array([2,1,2,3,2,5])) --> @nalift([0.5,2.0,1.5,NA,NA,1.0])
    @fact broadcast(/, DataCubes.wrap_array([1,2,3,4,5,6]), @nalift([2,1,2,2,NA,NA])) --> @nalift([0.5,2.0,1.5,2.0,NA,NA])
    @fact broadcast(/, @nalift([1,2,3,NA,NA,5]), 2) --> @nalift([0.5,1.0,1.5,NA,NA,2.5])
    @fact broadcast(/, @nalift([1,2,3,NA,NA,5]), 2.0) --> @nalift([0.5,1.0,1.5,NA,NA,2.5])
    @fact @nalift([1,2,3,NA,NA,5]) / 2 --> @nalift([0.5,1.0,1.5,NA,NA,2.5])
    @fact @nalift([1,2,3,NA,NA,5]) / 2.0 --> @nalift([0.5,1.0,1.5,NA,NA,2.5])
    @fact 1/@nalift([1,2,NA]) --> @nalift([1.0,0.5,NA])
    @fact 1.0/@nalift([1,2,NA]) --> @nalift([1.0,0.5,NA])
    @fact broadcast(/, 1, @nalift([1,2,NA])) --> @nalift([1.0,0.5,NA])
    @fact broadcast(/, 1.0, @nalift([1,2,NA])) --> @nalift([1.0,0.5,NA])
    @fact broadcast(^, e, @nalift([1,2,NA])) --> map(x -> isnull(x) ? Nullable{Float64}() : Nullable(e^x.value), @nalift([1,2,NA]))
    @fact broadcast(^, e, DataCubes.wrap_array([1.0,2.0])) --> DataCubes.wrap_array(map(x->e^x, [1.0,2.0]))
    @fact broadcast(^, e, DataCubes.wrap_array(Real[1.0,2])) --> DataCubes.wrap_array(map(x->e^x, Real[1.0,2]))
    @fact @nalift([1,3,NA]) == @nalift([1,NA,3]) --> false
    @fact @nalift([1,NA,3]) == @nalift([1,3,NA]) --> false
    @fact @nalift([1,NA,3]) == @nalift([1,NA,3]) --> true
    @fact @nalift([1,NA,3]) == @nalift([1,NA,2]) --> false
    @fact @nalift([1.0,3.0,NA]) == @nalift([1.0,NA,3.0]) --> false
    @fact @nalift([1.0,NA,3.0]) == @nalift([1.0,3.0,NA]) --> false
    @fact @nalift([1.0,NA,3.0]) == @nalift([1.0,NA,3.0]) --> true
    @fact @nalift([1.0,NA,3.0]) == @nalift([1.0,NA,2.0]) --> false
    @fact larr(a=[1,2,3]) + Nullable(1) --> larr(a=[2,3,4])
    @fact larr(a=1.0*[1,2,3]) + Nullable(1) --> larr(a=1.0*[2,3,4])
    @fact larr(a=1.0*[1,2,3]) + Nullable(1.0) --> larr(a=1.0*[2,3,4])
    @fact broadcast(+, larr(a=1.0*[1,2,3]), Nullable(1.0)) --> larr(a=1.0*[2,3,4])
    @fact broadcast(+, larr(a=1.0*[1,2,3]), Nullable{Int}()) --> larr(a=fill(Nullable{Float64}(),3))
    @fact broadcast(*, larr(a=1.0*[1,2,3]), Nullable(2.0)) --> larr(a=2.0*[1,2,3])
    @fact larr(a=1.0*[1,2,3]) * Nullable(2.0) --> larr(a=2.0*[1,2,3])
    @fact Nullable(2.0) * larr(a=1.0*[1,2,3]) --> larr(a=2.0*[1,2,3])
    @fact broadcast(*, Nullable(2), larr(a=1.0*[1,2,3])) --> larr(a=2.0*[1,2,3])
    @fact broadcast(*, Nullable{Int}(), larr(a=1.0*[1,2,3])) --> larr(a=fill(Nullable{Float64}(),3))
    @fact larr(a=[1,2,3]) + 1 --> larr(a=[2,3,4])
    @fact larr(a=1.0*[1,2,3]) + 1 --> larr(a=1.0*[2,3,4])
    @fact larr(a=1.0*[1,2,3]) * im --> larr(a=im*1.0*[1,2,3])
    @fact larr(a=1.0*[1,2,3]) / im --> larr(a=-im*1.0*[1,2,3])
    @fact larr(a=1.0*[1,2,3]) * Nullable(im) --> larr(a=1.0im*[1,2,3])
    @fact larr(a=1.0*[1,2,3]) / Nullable(im) --> larr(a=-1.0im*[1,2,3])
    @fact im * larr(a=1.0*[1,2,3]) --> larr(a=im*[1,2,3]*1.0)
    @fact im / larr(a=1.0*[1,2,3]) --> larr(a=im./(1.0*[1,2,3]))
    @fact Nullable(im) * larr(a=1.0*[1,2,3]) --> larr(a=im.*[1,2,3])
    @fact Nullable(im) / larr(a=1.0*[1,2,3]) --> larr(a=im./(1.0*[1,2,3]))
    @fact im * nalift(1.0*[1,2,3]) --> nalift(im*[1,2,3]*1.0)
    @fact im / nalift(1.0*[1,2,3]) --> nalift(im./(1.0*[1,2,3]))
    @fact Nullable(im) * nalift(1.0*[1,2,3]) --> nalift(im.*[1,2,3])
    @fact Nullable(im) / nalift(1.0*[1,2,3]) --> nalift(im./(1.0*[1,2,3]))
    @fact broadcast(*, nalift(1.0*[1,2,3]), im) --> nalift(im*1.0*[1,2,3])
    @fact broadcast(/, nalift(1.0*[1,2,3]), im) --> nalift(-im*1.0*[1,2,3])
    @fact nalift(1.0*[1,2,3]) * im --> nalift(im*1.0*[1,2,3])
    @fact nalift(1.0*[1,2,3]) / im --> nalift(-im*1.0*[1,2,3])
    @fact broadcast(*, wrap_array(1.0*[1,2,3]), im) --> wrap_array(im*1.0*[1,2,3])
    @fact broadcast(/, wrap_array(1.0*[1,2,3]), im) --> wrap_array(-im*1.0*[1,2,3])
    @fact larr(a=1.0*[1,2,3]) + 1.0 --> larr(a=1.0*[2,3,4])
    @fact broadcast(+, larr(a=1.0*[1,2,3]), 1.0) --> larr(a=1.0*[2,3,4])
    @fact broadcast(*, larr(a=1.0*[1,2,3]), 2.0) --> larr(a=2.0*[1,2,3])
    @fact larr(a=1.0*[1,2,3]) * 2.0 --> larr(a=2.0*[1,2,3])
    @fact 2.0 * larr(a=1.0*[1,2,3]) --> larr(a=2.0*[1,2,3])
    @fact broadcast(*, 2, larr(a=1.0*[1,2,3])) --> larr(a=2.0*[1,2,3])
    @fact darr(a=[1,2,3]) + 1 --> darr(a=[2,3,4])
    @fact darr(a=1.0*[1,2,3]) + 1 --> darr(a=1.0*[2,3,4])
    @fact darr(a=1.0*[1,2,3]) + 1.0 --> darr(a=1.0*[2,3,4])
    @fact broadcast(+, darr(a=1.0*[1,2,3]), 1.0) --> darr(a=1.0*[2,3,4])
    @fact broadcast(*, darr(a=1.0*[1,2,3]), 2.0) --> darr(a=2.0*[1,2,3])
    @fact darr(a=1.0*[1,2,3]) * 2.0 --> darr(a=2.0*[1,2,3])
    @fact 2.0 * darr(a=1.0*[1,2,3]) --> darr(a=2.0*[1,2,3])
    @fact broadcast(*, 2, darr(a=1.0*[1,2,3])) --> darr(a=2.0*[1,2,3])
    @fact darr(a=[1,2,3]) / 1 --> darr(a=[1,2,3])
    @fact darr(a=1.0*[1,2,3]) / 2 --> darr(a=[1,2,3]/2.0)
    @fact darr(a=1.0*[1,2,3]) / Nullable(2) --> darr(a=[1,2,3]/2.0)
    @fact broadcast(/, darr(a=1.0*[1,2,3]), 2) --> darr(a=[1,2,3]/2.0)
    @fact broadcast(/, darr(a=1.0*[1,2,3]), Nullable(2)) --> darr(a=[1,2,3]/2.0)
    @fact darr(a=1.0*[1,2,3]) / 3.0 --> darr(a=[1,2,3]/3.0)
    @fact broadcast(/, darr(a=1.0*[1,2,3]), 3.0) --> darr(a=[1,2,3]/3.0)
    @fact broadcast(/, darr(a=1.0*[1,2,3]), 2.0) --> darr(a=[1,2,3]/2.0)
    @fact darr(a=1.0*[1,2,3]) / 2.0 --> darr(a=[1,2,3]/2.0)
    @fact broadcast(/, darr(a=1.0*[1,2,3]), Nullable(2.0)) --> darr(a=[1,2,3]/2.0)
    @fact darr(a=1.0*[1,2,3]) / Nullable(2.0) --> darr(a=[1,2,3]/2.0)
    @fact broadcast(/, darr(a=[1,2,3]), darr(a=[3,2,1])) --> darr(a=[1,2,3]./[3,2,1])
    @fact 2.0 / darr(a=1.0*[1,2,3]) --> darr(a=2.0./[1,2,3])
    @fact broadcast(/, 2.0, darr(a=1.0*[1,2,3])) --> darr(a=2.0./[1,2,3])
    @fact Nullable(2.0) / darr(a=1.0*[1,2,3]) --> darr(a=2.0./[1,2,3])
    @fact broadcast(/, Nullable(2.0), darr(a=1.0*[1,2,3])) --> darr(a=2.0./[1,2,3])
    @fact darr(a=[1 2 3;4 5 6],b=1.0*[1 2 3;4 5 6])+darr(b=[1 2 3;4 5 6],a=[11 12 13;14 15 16]) --> darr(a=[12 14 16;18 20 22],b=2.0*[1 2 3;4 5 6])
    @fact darr(a=[1 2 3;4 5 6],b=1.0*[1 2 3;4 5 6])+larr(b=[1 2 3;4 5 6],a=[11 12 13;14 15 16],axis1=[:X,:Y]) --> larr(a=[12 14 16;18 20 22],b=2.0*[1 2 3;4 5 6],axis1=[:X,:Y])
    @fact darr(a=[1 2 3;4 5 6],b=1.0*[1 2 3;4 5 6])+larr(b=[1 2 3;4 5 6],a=[11 12 13;14 15 16],axis1=[:X,:Y]) --> larr(a=[12 14 16;18 20 22],b=2.0*[1 2 3;4 5 6],axis1=[:X,:Y])
    @fact larr(b=[1 2 3;4 5 6],a=[11 12 13;14 15 16],axis1=[:X,:Y]) + darr(a=[1 2 3;4 5 6],b=1.0*[1 2 3;4 5 6]) --> larr(b=2.0*[1 2 3;4 5 6],a=[12 14 16;18 20 22],axis1=[:X,:Y])
    @fact_throws larr(b=[1 2 3;4 5 6],a=[11 12 13;14 15 16],axis1=[:X,:Y]) + larr(a=[1 2 3;4 5 6],b=1.0*[1 2 3;4 5 6])
    @fact darr(a=[1.0 2.0]) * darr(a=[2.0,3.0]) --> darr(a=[8.0])
    @fact darr(a=[1.0 2.0;3.0 4.0]) / darr(a=[2.0 3.0]) --> darr(a=[1.0 2.0;3.0 4.0]/[2.0 3.0])
    @fact nalift([1.0 2.0]) * nalift([2.0,3.0]) --> nalift([8.0])
    @fact nalift([1.0 2.0;3.0 4.0]) / nalift([2.0 3.0]) --> nalift([1.0 2.0;3.0 4.0]/[2.0 3.0])
    @fact wrap_array(FloatNAArray([1.0 2.0]) * FloatNAArray([2.0,3.0])) --> wrap_array(FloatNAArray([8.0]))
    @fact wrap_array(FloatNAArray([1.0 2.0;3.0 4.0]) / FloatNAArray([2.0 3.0])) --> wrap_array(FloatNAArray([1.0 2.0;3.0 4.0]/[2.0 3.0]))
    @fact vcat(nalift([1 2 3]), nalift([4 5 6])) --> nalift([1 2 3;4 5 6])
    @fact hcat(nalift([1 2 3]), nalift([4 5 6])) --> nalift([1 2 3 4 5 6])
    @fact broadcast(/, nalift([1.0,2.0,3.0]), wrap_array([1,2,3])) --> nalift([1.0,1.0,1.0])
    @fact broadcast(/, wrap_array([1,2,3]), nalift([1.0,2.0,3.0]))  --> nalift([1.0,1.0,1.0])
    @fact broadcast(/, wrap_array([1,2,3]), wrap_array([1.0,2.0,3.0]))  --> wrap_array([1.0,1.0,1.0])
    @fact 3 * wrap_array([1,2,3]) --> wrap_array([3,6,9])
    @fact wrap_array([1,2,3]) * 3 --> wrap_array([3,6,9])
    # currently a matrix multiplication between non floating arrays are not supported.
    @fact larr(a=1.0*[1 2;3 4]) * (1.0 * larr(a=[5 6;7 8])) --> larr(a=[19.0 22.0;43.0 50.0])
    @fact_throws larr(a=[1 2;3 4]) * 1.0 * larr(a=[5 6;7 8])
    @fact_throws larr(a=[1 2;3 4]) * larr(a=[5 6;7 8])
    @fact_throws larr(rand(2,3),axis=[:a,:b],axis=['x','y','z']) * larr(rand(3,4),axis=[:m,:n,:p],axis=[10,11,12,13])
    @fact size(larr(rand(2,3),axis=[:a,:b],axis=['x','y','z']) * larr(rand(3,4),axis=['x','y','z'],axis=[10,11,12,13])) --> (2,4)
    @fact size(larr(rand(2,3),axis=[:a,:b],axis=['x','y','z']) / larr(rand(4,3),axis=[10,11,12,13],axis=['x','y','z'])) --> (2,4)
    @fact_throws size(larr(rand(2,3),axis=[:a,:b],axis=['x','y','z']) / larr(rand(4,3),axis=[10,11,12,13]))
  end
  context("ifelse tests") do
    @fact ifelse(@nalift([true,false,NA]), 3, 5) --> @nalift([3,5,NA])
    @fact ifelse(@nalift([true,false,NA]), Nullable(3), Nullable(5)) --> @nalift([3,5,NA])
    @fact ifelse(@nalift([true,false,NA]), Nullable(3), 5) --> @nalift([3,5,NA])
    @fact ifelse(@nalift([true,false,NA]), 3, Nullable(5)) --> @nalift([3,5,NA])
    @fact ifelse(@nalift([true,true,false,NA]), @nalift([1,2,3,4]), 5) --> @nalift([1,2,5,NA])
    @fact ifelse(@nalift([true,true,false,NA]), 5, @nalift([1,2,3,4])) --> @nalift([5,5,3,NA])
    @fact ifelse(@nalift([true,true,false,NA]), 1.0*nalift([1,2,3,4]), 5) --> 1.0*@nalift([1,2,5,NA])
    @fact ifelse(@nalift([true,true,false,NA]), 5, 1.0*nalift([1,2,3,4])) --> 1.0*@nalift([5,5,3,NA])
    @fact ifelse(@nalift([true,true,false,NA]), @nalift([1,2,3,4]), 5) --> @nalift([1,2,5,NA])
    @fact ifelse(@nalift([true,true,false,NA]), Nullable(5), @nalift([1,2,3,4])) --> @nalift([5,5,3,NA])
    @fact ifelse(@nalift([true,true,false,NA]), 1.0*nalift([1,2,3,4]), Nullable(5)) --> 1.0*@nalift([1,2,5,NA])
    @fact ifelse(@nalift([true,true,false,NA]), Nullable(5), 1.0*nalift([1,2,3,4])) --> 1.0*@nalift([5,5,3,NA])
    @fact ifelse(@nalift([true,true,false,NA]), Nullable(5), 1.0*[1,2,3,4]) --> 1.0*@nalift([5,5,3,NA])
    @fact ifelse(@nalift([true,true,false,NA]), [1,2,3,4], Nullable(5)) --> @nalift([1,2,5,NA])
    @fact ifelse(@nalift([true,true,false,NA]), 5, [1,2,3,4]) --> @nalift([5,5,3,NA])
    @fact ifelse(@nalift([true,true,false,NA]), 1.0*[1,2,3,4], 5) --> 1.0*@nalift([1,2,5,NA])
    @fact ifelse(@nalift([true,true,false,NA]), 1.0*AbstractArrayWrapper([1,2,3,4]), 5) --> 1.0*@nalift([1,2,5,NA])
    @fact ifelse(@nalift([true,true,false,NA]), 5, 1.0*[1,2,3,4]) --> 1.0*@nalift([5,5,3,NA])
    @fact ifelse(@nalift([true,true,false,NA]), 5, 1.0*AbstractArrayWrapper([1,2,3,4])) --> 1.0*@nalift([5,5,3,NA])
    @fact ifelse(@nalift([true,true,false,NA]), [1,2,3,4], Nullable(5)) --> @nalift([1,2,5,NA])
    @fact ifelse(@nalift([true,true,false,NA]), Nullable(5), [1,2,3,4]) --> @nalift([5,5,3,NA])
    @fact ifelse(@nalift([true,true,false,NA]), 1.0*[1,2,3,4], Nullable(5)) --> 1.0*@nalift([1,2,5,NA])
    @fact ifelse(@nalift([true,true,false,NA]), 1.0*AbstractArrayWrapper([1,2,3,4]), Nullable(5)) --> 1.0*@nalift([1,2,5,NA])
    @fact ifelse(@nalift([true,true,false,NA]), Nullable(5), 1.0*[1,2,3,4]) --> 1.0*@nalift([5,5,3,NA])
    @fact ifelse(@nalift([true,true,false,NA]), Nullable(5), 1.0*AbstractArrayWrapper([1,2,3,4])) --> 1.0*@nalift([5,5,3,NA])
    @fact ifelse(@nalift([true true false NA]), AbstractArrayWrapper([1 2 3 4]), AbstractArrayWrapper([11 12 13 14])*1.0) --> 1.0*@nalift([1 2 13 NA])
    @fact ifelse(@nalift([true true false NA]), nalift([1 2 3 4]), [11 12 13 14]*1.0) --> 1.0*@nalift([1 2 13 NA])
    @fact ifelse(@nalift([true true false NA]), nalift([1 2 3 4]), AbstractArrayWrapper([11 12 13 14])*1.0) --> 1.0*@nalift([1 2 13 NA])
    @fact ifelse(@nalift([true true false NA]), [1 2 3 4], nalift([11 12 13 14])*1.0) --> 1.0*@nalift([1 2 13 NA])
    @fact ifelse(@nalift([true true false NA]), AbstractArrayWrapper([1 2 3 4]), nalift([11 12 13 14])*1.0) --> 1.0*@nalift([1 2 13 NA])
    @fact ifelse(@nalift([true true false NA]), nalift([1 2 3 4]), nalift([11 12 13 14])*1.0) --> 1.0*@nalift([1 2 13 NA])

    @fact ifelse(DataCubes.AbstractArrayWrapper([true,false,false]), 3, 5) --> @nalift([3,5,5])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,false,false]), Nullable(3), 5) --> @nalift([3,5,5])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,false,false]), 3, Nullable(5)) --> @nalift([3,5,5])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,false,false]), Nullable(3), Nullable(5)) --> @nalift([3,5,5])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,true,false,false]), @nalift([1,2,3,4]), 5) --> @nalift([1,2,5,5])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,true,false,false]), [1,2,3,4], 5) --> @nalift([1,2,5,5])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,true,false,false]), AbstractArrayWrapper([1,2,3,4]), 5) --> @nalift([1,2,5,5])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,true,false,false]), 5, @nalift([1,2,3,4])) --> @nalift([5,5,3,4])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,true,false,false]), 5, [1,2,3,4]) --> @nalift([5,5,3,4])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,true,false,false]), 5, AbstractArrayWrapper([1,2,3,4])) --> @nalift([5,5,3,4])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,true,false,false]), 1.0*nalift([1,2,3,4]), 5) --> 1.0*nalift([1,2,5,5])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,true,false,false]), 5, 1.0*nalift([1,2,3,4])) --> 1.0*nalift([5,5,3,4])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,true,false,false]), @nalift([1,2,3,4]), 5) --> @nalift([1,2,5,5])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,true,false,false]), Nullable(5), @nalift([1,2,3,4])) --> @nalift([5,5,3,4])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,true,false,false]), Nullable(5), [1,2,3,4]) --> @nalift([5,5,3,4])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,true,false,false]), Nullable(5), AbstractArrayWrapper([1,2,3,4])) --> @nalift([5,5,3,4])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,true,false,false]), 1.0*nalift([1,2,3,4]), Nullable(5)) --> 1.0*nalift([1,2,5,5])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,true,false,false]), 1.0*[1,2,3,4], Nullable(5)) --> 1.0*nalift([1,2,5,5])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,true,false,false]), 1.0*AbstractArrayWrapper([1,2,3,4]), Nullable(5)) --> 1.0*nalift([1,2,5,5])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,true,false,false]), Nullable(5), 1.0*nalift([1,2,3,4])) --> 1.0*nalift([5,5,3,4])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,true,false,false]), [1,2,3,4], Nullable(5)) --> @nalift([1,2,5,5])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,true,false,false]), 5, [1,2,3,4]) --> @nalift([5,5,3,4])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,true,false,false]), 1.0*[1,2,3,4], 5) --> 1.0*nalift([1,2,5,5])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,true,false,false]), 5, 1.0*[1,2,3,4]) --> 1.0*nalift([5,5,3,4])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,true,false,false]), [1,2,3,4], Nullable(5)) --> @nalift([1,2,5,5])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,true,false,false]), Nullable(5), [1,2,3,4]) --> @nalift([5,5,3,4])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,true,false,false]), 1.0*[1,2,3,4], Nullable(5)) --> 1.0*nalift([1,2,5,5])
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,true,false,false]), Nullable(5), 1.0*[1,2,3,4]) --> 1.0*nalift([5,5,3,4])

    @fact ifelse(true, @nalift([1,2,3,4]), 5) --> @nalift([1,2,3,4])
    @fact ifelse(true, AbstractArrayWrapper([1,2,3,4]), 5) --> @nalift([1,2,3,4])
    @fact ifelse(Nullable(true), @nalift([1,2,3,4]), 5) --> @nalift([1,2,3,4])
    @fact ifelse(false, @nalift([1,2,3,4]), 5) --> @nalift([5,5,5,5])
    @fact ifelse(Nullable(false), @nalift([1,2,3,4]), 5) --> @nalift([5,5,5,5])
    @fact ifelse(Nullable(false), AbstractArrayWrapper([1,2,3,4]), 5) --> @nalift([5,5,5,5])
    @fact ifelse(Nullable{Bool}(), @nalift([1,2,3,4]), 5) --> @nalift(fill(Nullable{Int}(),4))
    @fact ifelse(true, 5, @nalift([1,2,3,4])) --> @nalift([5,5,5,5])
    @fact ifelse(Nullable(true), 5, @nalift([1,2,3,4])) --> @nalift([5,5,5,5])
    @fact ifelse(false, 5, @nalift([1,2,3,4])) --> @nalift([1,2,3,4])
    @fact ifelse(false, 5, wrap_array([1,2,3,4])) --> @nalift([1,2,3,4])
    @fact ifelse(Nullable(false), 5, @nalift([1,2,3,4])) --> @nalift([1,2,3,4])
    @fact ifelse(Nullable(false), 5, wrap_array([1,2,3,4])) --> @nalift([1,2,3,4])
    @fact ifelse(Nullable{Bool}(), 5, @nalift([1,2,3,4])) --> @nalift(fill(Nullable{Int}(),4))
    @fact ifelse(true, 1.0*nalift([1,2,3,4]), 5) --> 1.0*@nalift([1,2,3,4])
    @fact ifelse(Nullable(true), 1.0*nalift([1,2,3,4]), 5) --> 1.0*@nalift([1,2,3,4])
    @fact ifelse(false, 1.0*nalift([1,2,3,4]), 5) --> @nalift([5,5,5,5])
    @fact ifelse(Nullable(false), 1.0*nalift([1,2,3,4]), 5) --> @nalift([5,5,5,5])
    @fact ifelse(Nullable{Bool}(), 1.0*nalift([1,2,3,4]), 5) --> @nalift(fill(Nullable{Float64}(),4))
    @fact ifelse(true, 5, 1.0*nalift([1,2,3,4])) --> @nalift([5,5,5,5])
    @fact ifelse(Nullable(true), 5, 1.0*nalift([1,2,3,4])) --> @nalift([5,5,5,5])
    @fact ifelse(false, 5, 1.0*nalift([1,2,3,4])) --> 1.0*nalift([1,2,3,4])
    @fact ifelse(Nullable(false), 5, 1.0*nalift([1,2,3,4])) --> 1.0*nalift([1,2,3,4])
    @fact ifelse(Nullable{Bool}(), 5, 1.0*nalift([1,2,3,4])) --> 1.0*nalift(fill(Nullable{Int}(), 4))
    @fact ifelse(true, Nullable(5), @nalift([1,2,3,4])) --> @nalift([5,5,5,5])
    @fact ifelse(Nullable(true), Nullable(5), @nalift([1,2,3,4])) --> @nalift([5,5,5,5])
    @fact ifelse(false, Nullable(5), @nalift([1,2,3,4])) --> @nalift([1,2,3,4])
    @fact ifelse(Nullable(false), Nullable(5), @nalift([1,2,3,4])) --> @nalift([1,2,3,4])
    @fact ifelse(Nullable{Bool}(), Nullable(5), @nalift([1,2,3,4])) --> @nalift(fill(Nullable{Int}(),4))
    @fact ifelse(true, 1.0*nalift([1,2,3,4]), Nullable(5)) --> 1.0*nalift([1,2,3,4])
    @fact ifelse(Nullable(true), 1.0*nalift([1,2,3,4]), Nullable(5)) --> 1.0*nalift([1,2,3,4])
    @fact ifelse(false, 1.0*nalift([1,2,3,4]), Nullable(5)) --> @nalift([5,5,5,5])
    @fact ifelse(Nullable(false), 1.0*nalift([1,2,3,4]), Nullable(5)) --> @nalift([5,5,5,5])
    @fact ifelse(Nullable{Bool}(), 1.0*nalift([1,2,3,4]), Nullable(5)) --> 1.0*@nalift(fill(Nullable{Int}(),4))
    @fact ifelse(true, Nullable(5), 1.0*nalift([1,2,3,4])) --> @nalift([5,5,5,5])
    @fact ifelse(Nullable(true), Nullable(5), 1.0*nalift([1,2,3,4])) --> @nalift([5,5,5,5])
    @fact ifelse(false, Nullable(5), 1.0*nalift([1,2,3,4])) --> 1.0*nalift([1,2,3,4])
    @fact ifelse(Nullable(false), Nullable(5), 1.0*nalift([1,2,3,4])) --> 1.0*nalift([1,2,3,4])
    @fact ifelse(Nullable{Bool}(), Nullable(5), 1.0*nalift([1,2,3,4])) --> 1.0*nalift(fill(Nullable{Int}(),4))
    @fact ifelse(true, 5, DataCubes.AbstractArrayWrapper([1,2,3,4])) --> @nalift([5,5,5,5])
    @fact ifelse(Nullable(true), 5, DataCubes.AbstractArrayWrapper([1,2,3,4])) --> @nalift([5,5,5,5])
    @fact ifelse(false, 5, DataCubes.AbstractArrayWrapper([1,2,3,4])) --> @nalift([1,2,3,4])
    @fact ifelse(Nullable(false), 5, DataCubes.AbstractArrayWrapper([1,2,3,4])) --> @nalift([1,2,3,4])
    @fact ifelse(Nullable{Bool}(), 5, DataCubes.AbstractArrayWrapper([1,2,3,4])) --> @nalift(fill(Nullable{Int}(),4))
    @fact ifelse(true, DataCubes.AbstractArrayWrapper([1,2,3,4]), Nullable(5)) --> @nalift([1,2,3,4])
    @fact ifelse(Nullable(true), DataCubes.AbstractArrayWrapper([1,2,3,4]), Nullable(5)) --> @nalift([1,2,3,4])
    @fact ifelse(false, Nullable(5), DataCubes.AbstractArrayWrapper([1,2,3,4])) --> @nalift([1,2,3,4])
    @fact ifelse(Nullable(false), Nullable(5), DataCubes.AbstractArrayWrapper([1,2,3,4])) --> @nalift([1,2,3,4])
    @fact ifelse(Nullable{Bool}(), Nullable(5), DataCubes.AbstractArrayWrapper([1,2,3,4])) --> @nalift(fill(Nullable{Int}(),4))
    @fact ifelse(true, DataCubes.AbstractArrayWrapper(1.0*[1,2,3,4]), Nullable(5)) --> 1.0*nalift([1,2,3,4])
    @fact ifelse(Nullable(true), DataCubes.AbstractArrayWrapper(1.0*[1,2,3,4]), Nullable(5)) --> 1.0*nalift([1,2,3,4])
    @fact ifelse(false, DataCubes.AbstractArrayWrapper(1.0*[1,2,3,4]), Nullable(5)) --> @nalift([5,5,5,5])
    @fact ifelse(Nullable(false), DataCubes.AbstractArrayWrapper(1.0*[1,2,3,4]), Nullable(5)) --> @nalift([5,5,5,5])
    @fact ifelse(Nullable{Bool}(), DataCubes.AbstractArrayWrapper(1.0*[1,2,3,4]), Nullable(5)) --> @nalift(fill(Nullable{Int}(),4))

    @fact ifelse(true, @nalift([1,2,NA,NA,5]), @nalift([10,NA,8,NA,6])) --> @nalift([1,2,NA,NA,5])
    @fact ifelse(Nullable(true), @nalift([1,2,NA,NA,5]), @nalift([10,NA,8,NA,6])) --> @nalift([1,2,NA,NA,5])
    @fact ifelse(false, @nalift([1,2,NA,NA,5]), @nalift([10,NA,8,NA,6])) --> @nalift([10,NA,8,NA,6])
    @fact ifelse(Nullable(false), @nalift([1,2,NA,NA,5]), @nalift([10,NA,8,NA,6])) --> @nalift([10,NA,8,NA,6])
    @fact ifelse(Nullable{Bool}(), @nalift([1,2,NA,NA,5]), @nalift([10,NA,8,NA,6])) --> @nalift(fill(Nullable{Int}(),5))
    @fact ifelse(DataCubes.AbstractArrayWrapper([true,false,true,false,true]), @nalift([1,2,NA,NA,5]), @nalift([10,NA,8,NA,6])) --> @nalift([1,NA,NA,NA,5])
    @fact ifelse(@nalift([true,NA,true,false,NA]), @nalift([1,2,NA,NA,5]), @nalift([10,NA,8,NA,6])) --> @nalift([1,NA,NA,NA,NA])

    @fact ifelse(true, @nalift([1,2,NA,NA,5]), AbstractArrayWrapper([10,9,8,7,6])) --> @nalift([1,2,NA,NA,5])
    @fact ifelse(true, @nalift([1,2,NA,NA,5]), [10,9,8,7,6]) --> @nalift([1,2,NA,NA,5])
    @fact ifelse(Nullable(true), @nalift([1,2,NA,NA,5]), AbstractArrayWrapper([10,9,8,7,6])) --> @nalift([1,2,NA,NA,5])
    @fact ifelse(Nullable(true), @nalift([1,2,NA,NA,5]), [10,9,8,7,6]) --> @nalift([1,2,NA,NA,5])
    @fact ifelse(false, @nalift([1,2,NA,NA,5]), AbstractArrayWrapper([10,9,8,7,6])) --> @nalift([10,9,8,7,6])
    @fact ifelse(Nullable(false), @nalift([1,2,NA,NA,5]), AbstractArrayWrapper([10,9,8,7,6])) --> @nalift([10,9,8,7,6])
    @fact ifelse(Nullable{Bool}(), @nalift([1,2,NA,NA,5]), AbstractArrayWrapper([10,9,8,7,6])) --> @nalift(fill(Nullable{Int}(),4))
    @fact ifelse(AbstractArrayWrapper([true,false,true,false,true]), @nalift([1,2,NA,NA,5]), AbstractArrayWrapper([10,9,8,7,6])) --> @nalift([1,9,NA,7,5])
    @fact ifelse(@nalift([true,NA,true,false,true]), @nalift([1,2,NA,NA,5]), AbstractArrayWrapper([10,9,8,7,6])) --> @nalift([1,NA,NA,7,5])

    @fact ifelse(true, AbstractArrayWrapper([1,2,3,4,5]), @nalift([10,NA,8,NA,6])) --> @nalift([1,2,3,4,5])
    @fact ifelse(true, [1,2,3,4,5], @nalift([10,NA,8,NA,6])) --> @nalift([1,2,3,4,5])
    @fact ifelse(Nullable(true), AbstractArrayWrapper([1,2,3,4,5]), @nalift([10,NA,8,NA,6])) --> @nalift([1,2,3,4,5])
    @fact ifelse(Nullable(true), [1,2,3,4,5], @nalift([10,NA,8,NA,6])) --> @nalift([1,2,3,4,5])
    @fact ifelse(false, AbstractArrayWrapper([1,2,3,4,5]), @nalift([10,NA,8,NA,6])) --> @nalift([10,NA,8,NA,6])
    @fact ifelse(Nullable(false), AbstractArrayWrapper([1,2,3,4,5]), @nalift([10,NA,8,NA,6])) --> @nalift([10,NA,8,NA,6])
    @fact ifelse(Nullable{Bool}(), AbstractArrayWrapper([1,2,3,4,5]), @nalift([10,NA,8,NA,6])) --> @nalift(fill(Nullable{Int}(),4))
    @fact ifelse(AbstractArrayWrapper([true,false,true,false,true]), AbstractArrayWrapper([1,2,3,4,5]), @nalift([10,NA,8,NA,6])) --> @nalift([1,NA,3,NA,5])
    @fact ifelse(AbstractArrayWrapper([true,false,true,false,true]), AbstractArrayWrapper([1,2,3,4,5]), [10,9,8,7,6]) --> @nalift([1,9,3,7,5])
    @fact ifelse(AbstractArrayWrapper([true,false,true,false,true]), [10,9,8,7,6], AbstractArrayWrapper([1,2,3,4,5])) --> @nalift([10,2,8,4,6])
    @fact ifelse(AbstractArrayWrapper([true,false,true,false,true]), [10,9,8,7,6], nalift([1,2,3,4,5])) --> @nalift([10,2,8,4,6])
    @fact ifelse(AbstractArrayWrapper([true,false,true,false,true]), [10,9,8,7,6], [1,2,3,4,5]) --> @nalift([10,2,8,4,6])
    @fact ifelse(AbstractArrayWrapper([true,false,true,false,true]), nalift([10,9,8,7,6]), [1,2,3,4,5]) --> @nalift([10,2,8,4,6])
    @fact ifelse(@nalift([NA,NA,true,false,true]), AbstractArrayWrapper([1,2,3,4,5]), @nalift([10,NA,8,NA,6])) --> @nalift([NA,NA,3,NA,5])

    @fact ifelse(true, AbstractArrayWrapper([1,2,3,4,5]), AbstractArrayWrapper([10,9,8,7,6])) --> @nalift([1,2,3,4,5])
    @fact ifelse(true, AbstractArrayWrapper([1,2,3,4,5]), [10,9,8,7,6]) --> @nalift([1,2,3,4,5])
    @fact ifelse(true, [1,2,3,4,5], AbstractArrayWrapper([10,9,8,7,6])) --> @nalift([1,2,3,4,5])
    @fact ifelse(Nullable(true), AbstractArrayWrapper([1,2,3,4,5]), AbstractArrayWrapper([10,9,8,7,6])) --> @nalift([1,2,3,4,5])
    @fact ifelse(false, AbstractArrayWrapper([1,2,3,4,5]), AbstractArrayWrapper([10,9,8,7,6])) --> @nalift([10,9,8,7,6])
    @fact ifelse(Nullable(false), AbstractArrayWrapper([1,2,3,4,5]), AbstractArrayWrapper([10,9,8,7,6])) --> @nalift([10,9,8,7,6])
    @fact ifelse(Nullable{Bool}(), AbstractArrayWrapper([1,2,3,4,5]), AbstractArrayWrapper([10,9,8,7,6])) --> @nalift(fill(Nullable{Int}(),4))
    @fact ifelse(AbstractArrayWrapper([true,false,true,false,true]), AbstractArrayWrapper([1,2,3,4,5]), AbstractArrayWrapper([10,9,8,7,6])) --> @nalift([1,9,3,7,5])
    @fact ifelse(@nalift([true,NA,true,false,true]), AbstractArrayWrapper([1,2,3,4,5]), AbstractArrayWrapper([10,9,8,7,6])) --> @nalift([1,NA,3,7,5])
  end
end

end


text:
using NeuroCore
using Test
using Unitful
using CoordinateTransformations
using Documenter
using ColorTypes

#using NeuroCore: SPQuat, RotMatrix, quat2mat, mat2qua

using NeuroCore: s, Hz, °, T

using NeuroCore.FieldProperties

include("./ColorChannels/ColorChannels.jl")
include("./SpatialAPI/SpatialAPI.jl")
include("./NeuroMetadata/NeuroMetadata.jl")

#=

@testset "dimensions" begin
    x = NeuroArray(rand(2,3,4,10);
                      left = 1:2,
                      anterior = 1:3,
                      superior=1:4,
                      time = range(1, stop=10, length=10));

    @test NeuroCore.indices(x) == (1:2, 1:3, 1:4, 1.0:1.0:10.0)

    @test NeuroCore.spatial_offset(x) == (1, 1, 1)
    @test NeuroCore.spatial_eltype(x) == (Int64, Int64, Int64)

    # TODO get rid of once we don't need AxisArrays

    @test dimnames(typeof(NeuroCore.Axis{:time}(1:2))) == (:time,)
end

@testset "Magnetization Transfer" begin
    mt = NeuroCore.MagnetizationTransferMetadata(true, 1,2,3, :FERMI, 4)
    @test @inferred(NeuroCore.mt_offset_frequency(mt)) == 1.0Hz
    @test @inferred(NeuroCore.mt_pulse_bandwidth(mt)) == 2.0Hz
    @test @inferred(NeuroCore.mt_npulses(mt)) == 3
    @test @inferred(NeuroCore.mt_pulse_duration(mt)) == 4.0s

    #=
    @test @inferred(mt_offset_frequency(mt)) == 1.0Hz
    @test @inferred(mt_pulse_bandwidth(mt)) == 2.0Hz
    @test @inferred(mt_npulses(mt)) == 3
    @test @inferred(mt_pulse_duration(mt)) == 4.0s
    =#
end

@testset "Spoiling Gradiant" begin
    m = Metadata()

    @testset "spoiling_rf_phase_increment" begin
        m.spoiling_rf_phase_increment = 1
        @test @inferred(NeuroCore.spoiling_rf_phase_increment(m)) == (1*°)
        NeuroCore.spoiling_rf_phase_increment!(m, 1)
        @test m.spoiling_rf_phase_increment == (1*°)
    end

    @testset "duration" begin
        m.duration = 2.0
        @test @inferred(duration(m)) == 2.0s
        duration!(m, 2)
        @test m.duration == 2.0s
    end

    @testset "stop_time" begin
        m.stop_time = 2.0
        @test @inferred(stop_time(m)) == 2.0s
        stop_time!(m, 2)
        m.stop_time == 2.0s
    end

    @testset "sampling_rate" begin
        m.sampling_rate = .5
        @test @inferred(sampling_rate(m)) == 0.5s^-1
        sampling_rate!(m, .5)
        @test m.sampling_rate == 0.5s^-1
    end

    @testset "spoiling_gradient_duration" begin
        m.spoiling_gradient_duration = 2
        @test @inferred(NeuroCore.spoiling_gradient_duration(m)) == 2s
        NeuroCore.spoiling_gradient_duration!(m, 2)
        @test m.spoiling_gradient_duration == 2s
    end

    @testset "flip_angle" begin
        m.flip_angle = 1
        @test @inferred(NeuroCore.flip_angle(m)) == (1* °)
        NeuroCore.flip_angle!(m, 1)
        @test m.flip_angle == (1*°)
    end

    @testset "magnetic_field_strength" begin
        m.magnetic_field_strength = 2
        @test @inferred(NeuroCore.magnetic_field_strength(m)) == (2*NeuroCore.Unitful.T)
        NeuroCore.magnetic_field_strength!(m, 2)
        @test m.magnetic_field_strength == (2*NeuroCore.Unitful.T)
    end

    @testset "slice_encoding_direction" begin
        NeuroCore.slice_encoding_direction!(m, 1)
        @test @inferred(slicedim(m)) == 1
        slicedim!(m, 1)
        @test m.slicedim == 1
    end

    @testset "phase_encoding_direction" begin
        NeuroCore.phase_encoding_direction!(m, 3)
        @test @inferred(phasedim(m)) == 3
        phasedim!(m, 3)
        @test m.phasedim == 3
    end

    @testset "echo_time" begin
        m.echo_time = 1
        @test NeuroCore.echo_time(m) == 1.0s
        NeuroCore.echo_time!(m, 1)
        @test m.echo_time == 1.0s
    end

    @testset "inversion_time" begin
        m.inversion_time = 2
        @test NeuroCore.inversion_time(m) == 2.0s
        NeuroCore.inversion_time!(m, 2)
        @test m.inversion_time == 2.0s
    end

    # TODO needs to be an array
    #m.slice_timing = 3
    #@test NeuroCore.slice_timing(m) == 3.0s

    @testset "dwell_time" begin
        m.dwell_time = 4
        @test @inferred(NeuroCore.dwell_time(m)) == 4.0s
        NeuroCore.dwell_time!(m, 4)
        @test m.dwell_time == 4.0s
    end

    @testset "acquisition_duration" begin
        m.acquisition_duration = 5
        @test @inferred(NeuroCore.acquisition_duration(m)) == 5.0s
        NeuroCore.acquisition_duration!(m, 5)
        @test m.acquisition_duration == 5.0s
    end

    @testset "delay_after_trigger" begin
        m.delay_after_trigger = 6
        @test @inferred(NeuroCore.delay_after_trigger(m)) == 6.0s
        NeuroCore.delay_after_trigger!(m, 6)
        @test m.delay_after_trigger == 6.0s
    end

    # TODO needs to be an array
    #m.volume_timing = 7
    #@test NeuroCore.volume_timing(m) == 7.0s

    @testset "repetition_time" begin
        m.repetition_time = 8
        @test @inferred(NeuroCore.repetition_time(m)) == 8.0s
        NeuroCore.repetition_time!(m, 8)
        @test m.repetition_time == 8.0s
    end

    @testset "effective_echo_spacing" begin
        m.effective_echo_spacing = 9
        @test @inferred(NeuroCore.effective_echo_spacing(m)) == 9.0s
        NeuroCore.effective_echo_spacing!(m, 9)
        @test m.effective_echo_spacing == 9.0s
    end

    @testset "total_readout_time" begin
        m.total_readout_time = 10
        @test @inferred(NeuroCore.total_readout_time(m)) == 10.0s
        NeuroCore.total_readout_time!(m, 10)
        @test m.total_readout_time == 10.0s
    end

    @testset "high_cutoff" begin
        m.high_cutoff = 10
        @test @inferred(NeuroCore.high_cutoff(m)) == 10.0Hz
        NeuroCore.high_cutoff!(m, 10)
        @test m.high_cutoff == 10.0Hz
    end

    @testset "low_cutoff" begin
        m.low_cutoff = 10
        @test @inferred(NeuroCore.low_cutoff(m)) == 10.0Hz
        NeuroCore.low_cutoff!(m, 10)
        @test m.low_cutoff == 10.0Hz
    end

    @testset "notch_filter" begin
        m.notch_filter = 10
        @test @inferred(NeuroCore.notch_filter(m)) == 10.0Hz
        NeuroCore.notch_filter!(m, 10)
        @test m.notch_filter == 10.0Hz
    end

    @testset "power_line_frequency" begin
        m.power_line_frequency = 5
        @test @inferred(NeuroCore.power_line_frequency(m)) == 5.0Hz
        NeuroCore.power_line_frequency!(m, 10)
        @test m.power_line_frequency == 10.0Hz
    end

    @testset "dewar_position" begin
        m.dewar_position = 1
        @test @inferred(NeuroCore.dewar_position(m)) == (1* °)
        NeuroCore.dewar_position!(m, 1)
        @test m.dewar_position == (1* °)
    end

    @testset "epoch_length" begin
        m.epoch_length = 1
        @test @inferred(NeuroCore.epoch_length(m)) == 1.0s
        NeuroCore.epoch_length!(m, 1)
        @test m.epoch_length == 1.0s
    end
end

x = NAPArray(ones(2,2,2,2), x = 1:2, y = 1:2, time=1:2, z = 1:2)
@test @inferred(spatial_order(x)) == (:x, :y, :z)

include("semantic_positions.jl")
include("orientation.jl")
include("encoding_directions.jl")
include("time_tests.jl")
include("contrast_ingredient_tests.jl")

=#

@testset "docs" begin
    doctest(NeuroCore)
end

ext: jl


text:
include("Hilbert.jl")
include("OffHilbert.jl")
include("Convolution.jl")
include("LowRankIntegralOperator.jl")
include("NonlocalOperators.jl")


text:
using RobotDynamics
using ForwardDiff
using FiniteDiff
using Test
const RD = RobotDynamics

"""
My Docs
"""
RD.@autodiff Base.@kwdef struct MyKwModel <: RD.AbstractFunction
    a::Float64 = 1.0
    b::Float64 = 1.0
    int::Int = 10
end
RD.state_dim(::MyKwModel) = 4
RD.control_dim(::MyKwModel) = 2
RD.output_dim(::MyKwModel) = 4

function RD.evaluate(model::MyKwModel, x, u)
    return model.a * x + [u;u * model.int] * model.b
end

model = MyKwModel()
@test model.a == 1
@test model.int == 10
@test hasfield(MyKwModel, :cfg)
model = MyKwModel(a = 10, int=4)
@test model.a == 10
@test model.int == 4

x,u = rand(model)
y = zeros(length(x))
J = zeros(4,6)
z = RD.KnotPoint(x, u, 0., 0.1)
RD.evaluate(model, x, u)
RD.jacobian!(RD.StaticReturn(), RD.ForwardAD(), model, J, y, z)
A = ForwardDiff.jacobian(x->RD.evaluate(model, x, u), x)
B = ForwardDiff.jacobian(u->RD.evaluate(model, x, u), u)
@test J ≈ [A B]


text:
@testset "PertBeta" begin
    @testset "basic" begin
        dist = PertBeta(4, 10, 20)
        @test minimum(dist) == 4
        @test maximum(dist) == 20
        @test mode(dist) == 10
        @test mean(dist) == (4 + 4*10 + 20)/6
    end

    @testset "$dist" for dist in (PertBeta(4, 10, 20), PertBeta(0.0,0.5,1.0), PertBeta(1.0, 1.01, 10.0), )
        samples = rand(dist, 1_000_000)
        @test minimum(dist) < minimum(samples)
        @test maximum(samples) < maximum(dist)
        @test mean(samples) ≈ mean(dist) atol=0.01
        @test std(samples) ≈ std(dist) atol=0.01
    end

    @testset "Constructor" begin
        @test_throws Exception PertBeta(0, 10, 5)  # mode must be less than max
        @test_throws Exception PertBeta(5, 1, 7)  # mode must be more than min
    end
end


text:
using Lens

function bubblesort{T}(a::AbstractArray{T,1})
    b = copy(a)
    isordered = false
    span = length(b)
    iter = 0
    while !isordered && span > 1
        lens(:start_of_loop, b, iter) # <--- lens here!!
        isordered = true
        for i in 2:span
            if b[i] < b[i-1]
                t = b[i]
                b[i] = b[i-1]
                b[i-1] = t
                isordered = false
            end
        end
        span -= 1
        iter += 1
    end
    lens(:after_loop, sorteddata=b, niters=iter) # <--- and here!!
    return b
end

register!(println, :print_data, :start_of_loop, false)
bubblesort([1,2,0,3])
register!(:start_of_loop, Listener(:print_data, print, true, true))
bubblesort([1,2,0,3])

clear_all_listeners!()
register!(:print_data, :start_of_loop, false) do data, index
  println("on $index-th iteration the first element is $(index[1])")
end

bubblesort([50,2,0,20])

clear_all_listeners!()
function many_bubbles()
  for i = 1:100
    bubblesort(rand(100))
  end
end

iters = capture(many_bubbles, [(:after_loop, [:niters])])
get(iters)


text:
#when optimizing over the state psi1 we sometimes have terms like <psi1 | psi2>
#can reuse a lot of environments
#(we asume psi2 remains invariant)

struct OvlEnv{S,C<:MPSBondTensor,D <: GenericMPSTensor} <: Cache
    above :: S

    ldependencies::Vector{D} #the data we used to calculate leftenvs/rightenvs
    rdependencies::Vector{D}

    leftenvs::Vector{C}
    rightenvs::Vector{C}
end

function params(above::S,below::S,leftstart::C,rightstart::C) where S <: Union{<:FiniteMPS,<:MPSComoving} where C <: GenericMPSTensor
    leftenvs = [leftstart]
    rightenvs = [rightstart]

    for i in 1:length(above)
        push!(leftenvs,similar(leftstart))
        push!(rightenvs,similar(rightstart))
    end

    t = below.AC[1];
    return OvlEnv{S,C,eltype(below)}(above,fill(t,length(below)),fill(t,length(below)),leftenvs,reverse(rightenvs))
end

function params(above::S,below::S) where S <:FiniteMPS
    params(above,below,l_LL(above),r_RR(above))
end
function params(above::S,below::S) where S <:MPSComoving
    above.left_gs == below.left_gs || throw(ArgumentError("left gs differs"))
    above.right_gs == below.right_gs || throw(ArgumentError("right gs differs"))

    params(above,below,l_LL(above),r_RR(above))
end
#notify the cache that we updated in-place, so it should invalidate the dependencies
function poison!(ca::OvlEnv,ind)
    ca.ldependencies[ind] = similar(ca.ldependencies[ind])
    ca.rdependencies[ind] = similar(ca.rdependencies[ind])
end

#rightenv[ind] will be contracteable with the tensor on site [ind]
function rightenv(ca::OvlEnv,ind,state)
    a = findfirst(i -> !(state.AR[i] === ca.rdependencies[i]), length(state):-1:(ind+1))
    a = a == nothing ? nothing : length(state)-a+1

    if a != nothing
        #we need to recalculate
        for j = a:-1:ind+1
            ca.rightenvs[j] = transfer_right(ca.rightenvs[j+1],ca.above.AR[j],state.AR[j])
            ca.rdependencies[j] = state.AR[j]
        end
    end

    return ca.rightenvs[ind+1]
end

function leftenv(ca::OvlEnv,ind,state)
    a = findfirst(i -> !(state.AL[i] === ca.ldependencies[i]), 1:(ind-1))

    if a != nothing
        #we need to recalculate
        for j = a:ind-1
            ca.leftenvs[j+1] = transfer_left(ca.leftenvs[j],ca.above.AL[j],state.AL[j])
            ca.ldependencies[j] = state.AL[j]
        end
    end

    return ca.leftenvs[ind]
end


text:
using Cassette

Cassette.@context CounterCtx;

const ternops = (
    (:fma, Core.Intrinsics.fma_float, 2), # 2 flops per FMA instruction
    (:muladd, Core.Intrinsics.muladd_float, 2), # 2 flops per muladd instruction
)

const binops = (
    (:add, Core.Intrinsics.add_float, 1),
    (:sub, Core.Intrinsics.sub_float, 1),
    (:mul, Core.Intrinsics.mul_float, 1),
    (:div, Core.Intrinsics.div_float, 1),
    (:rem, Core.Intrinsics.rem_float, 1),
)

const unops = (
    (:abs, Core.Intrinsics.abs_float, 1),
    (:neg, Core.Intrinsics.neg_float, 1),
    (:sqrt, Core.Intrinsics.sqrt_llvm, 1),
)

const ops = Iterators.flatten((ternops, binops, unops)) |> collect

const typs = (
    (Float32, :32),
    (Float64, :64),
)

function gen_count(ops, suffix)
    body = Expr(:block)
    for (name, op) in ops
        fieldname = Symbol(name, suffix)
        e = quote
            if op == $op
                ctx.metadata.$fieldname += 1
                return
            end
        end
        push!(body.args, e)
    end
    body
end

for (typ, suffix) in typs
    @eval function Cassette.prehook(ctx::CounterCtx,
                                    op::Core.IntrinsicFunction,
                                    ::$typ,
                                    ::$typ,
                                    ::$typ)
        $(gen_count(ternops, suffix))
    end

    @eval function Cassette.prehook(ctx::CounterCtx,
                                    op::Core.IntrinsicFunction,
                                    ::$typ,
                                    ::$typ)
        $(gen_count(binops, suffix))
    end

    @eval function Cassette.prehook(ctx::CounterCtx,
                                    op::Core.IntrinsicFunction,
                                    ::$typ)
        $(gen_count(unops, suffix))
    end
end


text:
using SIMD

function mygemm!(C, A, B,
                 ::Val{MC}, ::Val{NC}, ::Val{KC},
                 ::Val{MR}, ::Val{NR},
                ) where {MC, NC, KC, MR, NR}
  (m, n) = size(C)
  (~, k) = size(A)
  @assert mod(m, MR) == 0
  @assert mod(n, NR) == 0
  @inbounds begin
    for J = 1:NC:n
      J2 = min(n, J-1+NC)
      for P = 1:KC:k
        P2 = min(k, P-1+KC)
        for I = 1:MC:m
          I2 = min(n, I-1+MC)
          mygemm_NRxMR!(view(C, I:I2, J:J2), m,
                        view(A, I:I2, P:P2), m,
                        view(B, P:P2, J:J2), k,
                        Val(MR), Val(NR))
        end
      end
    end
  end
end

function mygemm_NRxMR!(C, ldc, A, lda, B, ldb, ::Val{MR}, ::Val{NR}) where {MR, NR}
  (m, n) = size(C)
  (~, k) = size(A)
  @inbounds for J = 1:NR:n
    for I = 1:MR:m
      @views k_mygemm(k,
                      pointer(C, I + (J-1) * m), ldc,
                      pointer(A, I), lda,
                      pointer(B, 1 + (J-1) * k), ldb,
                      Val(MR), Val(NR))
    end
  end
end

function k_mygemm(k,
                  C, ldc,
                  A, lda,
                  B, ldb,
                  ::Val{MR}, ::Val{NR}) where {MR, NR}
  T = eltype(A)
  NR > 4 && error("NR too big")

  NR>0 && (c1 = vload(Vec{MR, T}, C + 0ldc*sizeof(T)))
  NR>1 && (c2 = vload(Vec{MR, T}, C + 1ldc*sizeof(T)))
  NR>2 && (c3 = vload(Vec{MR, T}, C + 2ldc*sizeof(T)))
  NR>3 && (c4 = vload(Vec{MR, T}, C + 3ldc*sizeof(T)))

  @inbounds for p = 1:k
    αv = vload(Vec{MR, T}, A + (p-1)*lda*sizeof(T))

    β = unsafe_load(B + (p-1 + 0ldb)*sizeof(T))
    c1 = muladd(β, αv, c1)

    if NR > 1
      β = unsafe_load(B + (p-1 + 1ldb)*sizeof(T))
      c2 = muladd(β, αv, c2)
    end

    if NR > 2
      β = unsafe_load(B + (p-1 + 2ldb)*sizeof(T))
      c3 = muladd(β, αv, c3)
    end

    if NR > 3
      β = unsafe_load(B + (p-1 + 3ldb)*sizeof(T))
      c4 = muladd(β, αv, c4)
    end
  end
  NR>0 && vstore(c1, C + 0ldc*sizeof(T))
  NR>1 && vstore(c2, C + 1ldc*sizeof(T))
  NR>2 && vstore(c3, C + 2ldc*sizeof(T))
  NR>3 && vstore(c4, C + 3ldc*sizeof(T))
end

max_stars_repo_name: jkozdon/Julia-LAFF-On-PfHP

max_issues_repo_name: jkozdon/Julia-LAFF-On-PfHP

max_forks_repo_name: jkozdon/Julia-LAFF-On-PfHP


text:
@testset "Alphabets" begin
    import KnuthBendix.set_inversion!

    @test Alphabet{Char}() isa Alphabet{Char}
    @test_throws AssertionError Alphabet{Int}()
    @test_throws AssertionError Alphabet([1, 2, 3])
    @test_throws AssertionError Alphabet(['a', 'b', 'a'])

    let B = Alphabet(['a', 'b', 'c'])
        @test B isa Alphabet{Char}
        @test sprint(show, B) isa String
        @test sprint(show, MIME"text/plain"(), B) isa String
        @test length(B.letters) == 3 && length(B.inversions) == 3
        @test findfirst(i -> i != 0, B.inversions) === nothing

        @test B[1] == 'a' && B[2] == 'b' && B[3] == 'c'
        @test B['a'] == 1 && B['b'] == 2 && B['c'] == 3
        @test_throws DomainError B['d']
    end

    let A = Alphabet{Char}()
        @test sprint(show, A) isa String
        @test sprint(show, MIME"text/plain"(), A) isa String

        @test length(A.letters) == 0 && length(A.inversions) == 0

        push!(A, 'a', 'b', 'c')
        @test length(A.letters) == 3 && length(A.inversions) == 3
        @test findfirst(i -> i != 0, A.inversions) === nothing
        @test_throws ErrorException push!(A, 'a')

        @test A[1] == 'a' && A[2] == 'b' && A[3] == 'c'
        @test A['a'] == 1 && A['b'] == 2 && A['c'] == 3

        @test_throws ErrorException set_inversion!(A, 'd', 'e')
        @test_throws ErrorException set_inversion!(A, 'a', 'e')

        set_inversion!(A, 'a', 'b')
        @test A[-2] == 'a' && A[-1] == 'b'
        set_inversion!(A, 'b', 'c')
        @test A[-2] == 'c' && A[-3] == 'b'
        set_inversion!(A, 'a', 'c')
        @test A[-1] == 'c' && A[-3] == 'a'
        @test_throws DomainError A[-2]

        @test sprint(show, A) isa String
        @test sprint(show, MIME"text/plain"(), A) isa String

        @test A[A[1]] == 1

        A = Alphabet(["a₁", "a₁^-1"], [2, 1])
        w = Word([1,2,2])
        @test sprint(KnuthBendix.print_repr, w, A) == "a₁*a₁^-2"
    end

    @testset "Inverting using Alphabet" begin
        A = Alphabet(['a', 'b', 'A'])
        KnuthBendix.set_inversion!(A, 'a', 'A')

        w = KnuthBendix.Word([1, 2])
        @test_throws DomainError inv(A, w) # b is not invertible
        w = KnuthBendix.Word([1, 1, 3])
        @test inv(A, w) == Word([1, 3, 3]) # inv(a*a*c)
        @test inv(A, inv(A, w)) == w

        @test sprint(KnuthBendix.print_repr, w, A) == "a^2*A"
        @test sprint(KnuthBendix.print_repr, inv(A, w), A) == "a*A^2"

        B = Alphabet(["a", "a^-1", "c"], [2,1,0])
        w = KnuthBendix.Word([1])
        @test sprint(KnuthBendix.print_repr, w, B) == "a"
        w = KnuthBendix.Word([1, 1])
        @test sprint(KnuthBendix.print_repr, w, B) == "a^2"
        w = KnuthBendix.Word([2])
        @test sprint(KnuthBendix.print_repr, w, B) == "a^-1"
        w = KnuthBendix.Word([2,2])
        @test sprint(KnuthBendix.print_repr, w, B) == "a^-2"

        w = KnuthBendix.Word([3, 1])
        @test sprint(KnuthBendix.print_repr, w, B) == "c*a"
        w = KnuthBendix.Word([3, 1, 1])
        @test sprint(KnuthBendix.print_repr, w, B) == "c*a^2"
        w = KnuthBendix.Word([3, 2])
        @test sprint(KnuthBendix.print_repr, w, B) == "c*a^-1"
        w = KnuthBendix.Word([3, 2,2])
        @test sprint(KnuthBendix.print_repr, w, B) == "c*a^-2"
    end
end


text:
module Solution
export spinwords

spinwords(s::String) = join( map( w -> length(w)≥5 ? reverse(w) : w, split(s) ), " " )

end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
using ResourceManagement
using DataFrames
A = DisciplineLabor()
B = DisciplineLabor()

dflabor = ReadLaborTracker("C:\\Users\\junqueg\\Documents\\My Documents\\15. Programming\\Projects\\ResourceManagement.jl\\src\\TEAM_PLANNED_FWD24_NOV.csv");
vh, pv = _getEmployeePlannedHours(dflabor, "HIGA ANTHONY", 24)

vh

sum(sum.(vh))

typeof(pv)
for p in pv
    @show p
end

p=unique(pv)
for project in p
    @show project
end

p[1]
p[2]
p[3]
p[4]

Tony = DisciplineLabor("430300", "Higa Anthony",  24);

V1, p1 = fetchAndWritePlannedHours!(dflabor, "HIGA ANTHONY", 24, Tony);

sum(getFwdPlannedHours(Tony, ""))

sum(getFwdPlannedHours(Tony, "153804"))

dfAvail = ReadAvailHours("src\\UTILREPORT_FWD_NOV.csv")
A = _getAvailMonthHours(dfAvail, 24)

for i in A
    @show i[1]
end

dfRev = ReadAvailHours("C:\\Users\\junqueg\\Documents\\My Documents\\15. Programming\\Projects\\ResourceManagement.jl\\src\\TEAM_ACTUALPLAN_REV18.csv");

Tony = DisciplineLabor("430300", "Higa Anthony",  24)
Julie = DisciplineLabor("430300", "BRIONES JULITA",  24)

fetchAndWritePlannedHours!(dfRev, "HIGA ANTHONY", 18, Tony, :plan);

fetchAndWritePlannedHours!(dfRev, "HIGA ANTHONY", 18, Tony, :rev);
fetchAndWritePlannedHours!(dflabor, "HIGA ANTHONY", 24, Tony, :fwd);
fetchAndWriteRevActualHours!(dfRev, "HIGA ANTHONY", 18, Tony);

Julie = fetchAndWritePlannedHours!(dfRev, "BRIONES JULITA", 18, Julie, :rev);
Julie = fetchAndWritePlannedHours!(dflabor, "BRIONES JULITA", 24, Julie, :fwd);
fetchAndWriteRevActualHours!(dfRev, "BRIONES JULITA", 18, Julie);

writeAvailableFwdHours!(Tony, dfAvail, 24);
writeAvailableFwdHours!(Julie, dfAvail, 24);

Team1 = Tony + Julie

take(n::Int) = xs -> xs[1:min(n,end)]

t = TeamDump(Team1, :TotalFwdPlannedHours)
TeamDump(Team1, :TotalRevPlannedHours)
TeamDump(Team1, :TotalRevActualHours)
TeamDump(Team1, :FwdUtilization)

#pipe examples
TeamDump(Team1, :FwdUtilization) |> values |> collect .|> mean
TeamDump(Team1, :TotalFwdPlannedHours) |> values |> collect .|> sum

A = TeamDump(Team1, :TotalFwdPlannedHours) |> values |>collect |>minimum |> take(6)  # collect hours and take the minimum monthly values from the resources and take the first 6

B = TeamDump(Team1, :TotalFwdPlannedHours) |> values |>collect

TeamDump(Team1, :TotalRevActualHours) |> values |> take(3)

"TT = combine(Tony.FwdHoursForecast, months .=>sum)

size(Tony.FwdHoursForecast)[2]

TT = Tony.FwdHoursForecast[:, 2: end]
months= names(Tony.FwdHoursForecast)[2:end]

Vector(TT[1,:])

Vector(TT)

Vector(TT)

combine(Tony.FwdHoursForecast, months .=>sum)

Tony.FwdHoursForecast.Project

filter = Tony.FwdHoursForecast.Project .== "153804"

df = Tony.FwdHoursForecast

v = combine(df[filter, :] , months .=>sum)
v

v = Vector(v[1, :])

sum(getFwdPlannedHours(Tony, "153804"))

Tony.FwdHoursForecast.Project
Julie.FwdHoursForecast.Project

Tr = DisciplineLabor("430300", "Higa Anthony",  24);

Tr.FwdHoursForecast.Project

Tony.FwdHoursForecast ∪ Julie.FwdHoursForecast
S1 = Set(Julie.FwdHoursForecast.Project)

S2 = Set(Tony.FwdHoursForecast.Project)

Tony.FwdHoursForecast.Project

S1 ∪ S2

Tr.FwdHoursForecast.Project = String.(S1 ∪ S2)

Tr.FwdHoursForecast[!,:Project] = String.(S1 ∪ S2)

String.(S1 ∪ S2)

Tr.FwdHoursForecast

Tony.FwdHoursForecast

names(Tony.FwdHoursForecast)
Tony.FwdHoursForecast[:,2] .+ Julie.FwdHoursForecast[:,2]

Tr.FwdHoursForecast = vcat(Tony.FwdHoursForecast, Julie.FwdHoursForecast)

Tr.FwdHoursForecast.Project

sum(getFwdPlannedHours(Julie))
getFwdPlannedHours(Tony)
dflabor

sum(getFwdPlannedHours(Tr))


text:
"""
    $(FUNCTIONNAME)(model::HitranModel, grid::AbstractRange{<:Real}, pressure::Real, temperature::Real; wavelength_flag::Bool=false)

Calculate absorption cross-section at the given pressure, temperature, and grid of wavelengths
(or wavenumbers), and have the option to perform auto-differentiation

"""
function absorption_cross_section(model::AbstractCrossSectionModel,          # Model to use
                                  grid::Union{AbstractRange{<:Real}, AbstractArray}, # Wavelength [nm] or wavenumber [cm-1] grid
                                  pressure::Real,              # actual pressure [hPa]
                                  temperature::Real;
                                  autodiff::Bool=false,           # actual temperature [K]
                                  wavelength_flag::Bool=false)

    # This function takes in the "x-vector" along with the input model so that ForwardDiff will work
    function absorption_cross_section_autodiff(x ; model::AbstractCrossSectionModel = model, grid=grid, pressure=pressure,
                                               temperature=temperature, wavelength_flag=wavelength_flag)

        if length(x) !== 2
            @error "Must receive two cross-section parameters for auto-differentiation (p, T)" x
        end

        # Make sure that 𝐱 and parameter match
        @assert (pressure== x[1])
        @assert (temperature == x[2])

        cross_section = compute_absorption_cross_section(model, grid, x[1], x[2], wavelength_flag=wavelength_flag);

        return cross_section
    end

    if (autodiff)

        x = [pressure, temperature]

        result = DiffResults.JacobianResult(zeros(length(collect(grid))), x);
        ForwardDiff.jacobian!(result, absorption_cross_section_autodiff, x);

        return (result.value, result.derivs[1]);

    else
        return compute_absorption_cross_section(model, grid, pressure, temperature, wavelength_flag=wavelength_flag)
    end

end


text:
using Test

# Thin category
###############

x,y,z = Ob(FreeThinCategory, :x, :y, :z)
f, g = Hom(:f, x, y), Hom(:g, y, z)
@test dom(compose(f,g)) == x
@test codom(compose(f,g)) == z

x,y,z,zz = Ob(FreeThinSymmetricMonoidalCategory, :x, :y, :z, :zz)
f, g, h= Hom(:f, x, y), Hom(:g, y, z), Hom(:h, otimes(y,y), z)
j = Hom(:j, y⊗z, zz)
@test dom(compose(f,g)) == x
@test codom(compose(f,g)) == z
@test dom(otimes(f,g)) == otimes(x,y)
@test codom(otimes(f,g)) == otimes(y,z)
@test dom((f⊗id(y))⋅h) == x⊗y
@test codom((f⊗id(y))⋅h) == z
@test dom((f⊗g)⋅j) == x⊗y
@test codom((f⊗g)⋅j) == zz

# Preorder
##########

x,y,z = (El(FreePreorder.El, sym) for sym in (:x, :y, :z))
f, g = Leq(:f, x, y), Leq(:g, y, z)
@test lhs(reflexive(x)) == x
@test rhs(reflexive(x)) == x
@test lhs(transitive(f,g)) == x
@test rhs(transitive(f,g)) == z


text:
"""
    Paragraph() <: Block

[`Block`](#) for a text paragraph containing one or more
sentences (basically, a single observation in the textual dataset).
`data` is valid for `Paragraph` if it is of type string.

Example valid Paragraphs:

```julia
@test checkblock(Paragraph(), "Hello world!")
@test checkblock(Paragraph(), "Hello world!, How are you?")
```

You can create a random observation using [`mockblock`](#):

{cell=main}
```julia
using FastAI
FastAI.mockblock(Paragraph())
```

"""

struct Paragraph <: Block end

FastAI.checkblock(::Paragraph, ::String) = true
FastAI.mockblock(::Paragraph) = randstring(" ABCEEFGHIJKLMNOPQESRUVWXYZ 1234567890 abcdefghijklmnopqrstynwxyz\n\t.,", rand(10:40))


text:
include("crowding-distance.jl")

abstract type AbstractNSGA <: AbstractParameters end

# Abstracts for algorithm parameter
mutable struct NSGA2 <: AbstractNSGA
    N::Int
    η_cr::Float64
    p_cr::Float64
    η_m::Float64
    p_m::Float64
end

"""
    NSGA2(;
        N = 100,
        η_cr = 20,
        p_cr = 0.9,
        η_m = 20,
        p_m = 1.0 / D,
        information = Information(),
        options = Options(),
    )

Parameters for the metaheuristic NSGA-II.

Parameters:

- `N` Population size.
- `η_cr`  η for the crossover.
- `p_cr` Crossover probability.
- `η_m`  η for the mutation operator.
- `p_m` Mutation probability (1/D for D-dimensional problem by default).

To use NSGA2, the output from the objective function should be a 3-touple
`(f::Vector, g::Vector, h::Vector)`, where `f` contains the objective functions,
`g` and `h` are inequality, equality constraints respectively.

A feasible solution is such that `g_i(x) ≤ 0 and h_j(x) = 0`.

```julia
using Metaheuristics

# Dimension
D = 2

# Objective function
f(x) = ( x, [sum(x.^2) - 1], [0.0] )

# bounds
bounds = [-1 -1;
           1  1.0
        ]

# define the parameters (use `NSGA2()` for using default parameters)
nsga2 = NSGA2(N = 100, p_cr = 0.85)

# optimize
status = optimize(f, bounds, nsga2)

# show results
display(status)
```

"""
function NSGA2(;
    N = 100,
    η_cr = 20,
    p_cr = 0.9,
    η_m = 20,
    p_m = -1,
    information = Information(),
    options = Options(),
)

    parameters = NSGA2(N, promote( Float64(η_cr), p_cr, η_m, p_m )...)
    Algorithm(
        parameters,
        information = information,
        options = options,
    )

end

function update_state!(
    status::State,
    parameters::AbstractNSGA,
    problem::AbstractProblem,
    information::Information,
    options::Options,
    args...;
    kargs...
    )

    I = randperm(parameters.N)
    J = randperm(parameters.N)
    Q = empty(status.population)
    for i = 1:2:parameters.N

        pa = tournament_selection(status.population, I[i])
        pb = tournament_selection(status.population, J[i])

        offspring1, offspring2 = reproduction(pa, pb, parameters, problem)

        # save offsprings
        push!(Q, offspring1, offspring2)
    end

    status.population = vcat(status.population, Q)

    # non-dominated sort, crowding distance, elitist removing
    environmental_selection!(status.population, parameters)

end

function reproduction(pa, pb, parameters::AbstractNSGA, problem)
    # crossover and mutation
    c1, c2 = GA_reproduction(get_position(pa),
                             get_position(pb),
                             problem.bounds;
                             η_cr = parameters.η_cr,
                             p_cr = parameters.p_cr,
                             η_m = parameters.η_m,
                             p_m = parameters.p_m)

    # evaluate offspring
    create_solution(c1, problem), create_solution(c2, problem)
end

function environmental_selection!(population, parameters::AbstractNSGA)
    truncate_population!(population, parameters.N)
end

function initialize!(
    status,
    parameters::AbstractNSGA,
    problem::AbstractProblem,
    information::Information,
    options::Options,
    args...;
    kargs...
)
    D = size(problem.bounds, 2)

    if parameters.p_m < 0.0
        parameters.p_m = 1.0 / D
    end

    if options.iterations == 0
        options.iterations = 500
    end

    if options.f_calls_limit == 0
        options.f_calls_limit = options.iterations * parameters.N + 1
    end

    status = gen_initial_state(problem,parameters,information,options,status)

    status

end

function final_stage!(
    status::State,
    parameters::AbstractNSGA,
    problem::AbstractProblem,
    information::Information,
    options::Options,
    args...;
    kargs...
    )
    status.final_time = time()

end

function tournament_selection(P, a = rand(1:length(P)))
    # chose two different solutions at random
    b = rand(1:length(P))
    while a == b
        b = rand(1:length(P))
    end

    # perform selection
    P[a].rank < P[b].rank || (P[a].rank == P[b].rank && P[a].crowding > P[b].crowding ) ? P[a] : P[b]
end


text:
"""
    DataBunch(train::Flux.Data.DataLoader, valid::Flux.Data.DataLoader)

A `DataBunch` is a bunched `train` and `valid` (validation) dataloader.
"""
struct DataBunch
    train::Flux.Data.DataLoader
    valid::Flux.Data.DataLoader
end

"""
Get the `db.train` dataloader.
"""
train(db::DataBunch) = db.train
"""
Get the `db.valid` (validation) dataloader.
"""
valid(db::DataBunch) = db.valid


text:
using CuthillMcKee
using Test

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


text:
@enum(ProxALTime,
    Origin,
    Final,
    Normal,
)

abstract type AbstractTimeStep end
struct Current <: AbstractTimeStep end
struct Previous <: AbstractTimeStep end
struct Next <: AbstractTimeStep end

"""
    ProxALEvaluator{T, VI, VT, MT} <: AbstractNLPEvaluator

Evaluator wrapping a `ReducedSpaceEvaluator` for use inside the
decomposition algorithm implemented in [ProxAL.jl](https://github.com/exanauts/ProxAL.jl).

"""
mutable struct ProxALEvaluator{T, VI, VT, MT} <: AbstractNLPEvaluator
    inner::ReducedSpaceEvaluator{T, VI, VT, MT}
    s_min::VT
    s_max::VT
    nu::Int
    ng::Int
    # Augmented penalties parameters
    time::ProxALTime
    scale_objective::T
    τ::T
    λf::VT
    λt::VT
    ρf::T
    ρt::T
    pg_f::VT
    pg_ref::VT
    pg_t::VT
    # Buffer
    ramp_link_prev::VT
    ramp_link_next::VT
end
function ProxALEvaluator(
    nlp::ReducedSpaceEvaluator{T, VI, VT, MT},
    time::ProxALTime;
    τ=0.1, ρf=0.1, ρt=0.1, scale_obj=1.0,
) where {T, VI, VT, MT}
    nu = n_variables(nlp)
    ng = get(nlp, PS.NumberOfGenerators())

    s_min = xzeros(VT, ng)
    s_max = xones(VT, ng)
    λf = xzeros(VT, ng)
    λt = xzeros(VT, ng)

    pgf = xzeros(VT, ng)
    pgc = xzeros(VT, ng)
    pgt = xzeros(VT, ng)

    ramp_link_prev = xzeros(VT, ng)
    ramp_link_next = xzeros(VT, ng)

    return ProxALEvaluator(
        nlp, s_min, s_max, nu, ng, time, scale_obj, τ, λf, λt, ρf, ρt,
        pgf, pgc, pgt, ramp_link_prev, ramp_link_next,
    )
end
function ProxALEvaluator(
    pf::PS.PowerNetwork,
    time::ProxALTime;
    device=KA.CPU(),
    options...
)
    # Build network polar formulation
    model = PolarForm(pf, device)
    # Build reduced space evaluator
    x = initial(model, State())
    u = initial(model, Control())
    nlp = ReducedSpaceEvaluator(model; options...)
    return ProxALEvaluator(nlp, time)
end
function ProxALEvaluator(
    datafile::String;
    time::ProxALTime=Normal,
    device=KA.CPU(),
    options...
)
    nlp = ReducedSpaceEvaluator(datafile; device=device, options...)
    return ProxALEvaluator(nlp, time; options...)
end

n_variables(nlp::ProxALEvaluator) = nlp.nu + nlp.ng
n_constraints(nlp::ProxALEvaluator) = n_constraints(nlp.inner)

constraints_type(::ProxALEvaluator) = :inequality

# Getters
get(nlp::ProxALEvaluator, attr::AbstractNLPAttribute) = get(nlp.inner, attr)
get(nlp::ProxALEvaluator, attr::AbstractVariable) = get(nlp.inner, attr)
get(nlp::ProxALEvaluator, attr::PS.AbstractNetworkValues) = get(nlp.inner, attr)
get(nlp::ProxALEvaluator, attr::PS.AbstractNetworkAttribute) = get(nlp.inner, attr)

# Setters
function setvalues!(nlp::ProxALEvaluator, attr::PS.AbstractNetworkValues, values)
    setvalues!(nlp.inner, attr, values)
end
function transfer!(nlp::ProxALEvaluator, vm, va, pg, qg)
    transfer!(nlp.inner, vm, va, pg, qg)
end

# Initial position
function initial(nlp::ProxALEvaluator)
    u0 = initial(nlp.inner)
    s0 = copy(nlp.s_min)
    return [u0; s0]
end

# Bounds
function bounds(nlp::ProxALEvaluator, ::Variables)
    u♭, u♯ = bounds(nlp.inner, Variables())
    return [u♭; nlp.s_min], [u♯; nlp.s_max]
end
bounds(nlp::ProxALEvaluator, ::Constraints) = bounds(nlp.inner, Constraints())

function update!(nlp::ProxALEvaluator, w)
    u = @view w[1:nlp.nu]
    s = @view w[nlp.nu+1:end]
    conv = update!(nlp.inner, u)
    pg = get(nlp.inner, PS.ActivePower())
    # Update terms for augmented penalties
    if nlp.time != Origin
        nlp.ramp_link_prev .= nlp.pg_f .- pg .+ s
    end
    if nlp.time != Final
        nlp.ramp_link_next .= pg .- nlp.pg_t
    end
    return conv
end

## Update penalty terms
function update_multipliers!(nlp::ProxALEvaluator, ::Current, λt)
    copyto!(nlp.λf, λt)
end
function update_multipliers!(nlp::ProxALEvaluator, ::Next, λt)
    copyto!(nlp.λt, λt)
end

function update_primal!(nlp::ProxALEvaluator, ::Previous, pgk)
    copyto!(nlp.pg_f, pgk)
end
function update_primal!(nlp::ProxALEvaluator, ::Current, pgk)
    copyto!(nlp.pg_ref, pgk)
end
function update_primal!(nlp::ProxALEvaluator, ::Next, pgk)
    copyto!(nlp.pg_t, pgk)
end

## Objective
function objective(nlp::ProxALEvaluator, w)
    u = @view w[1:nlp.nu]
    s = @view w[nlp.nu+1:end]
    buffer = get(nlp.inner, PhysicalState())
    pg = get(nlp.inner, PS.ActivePower())
    # Operational costs
    c = nlp.scale_objective * cost_production(nlp.inner.model, buffer)
    # Augmented Lagrangian penalty
    c += 0.5 * nlp.τ * xnorm(pg .- nlp.pg_ref)^2
    if nlp.time != Origin
        c += dot(nlp.λf, nlp.ramp_link_prev)
        c += 0.5 * nlp.ρf * xnorm(nlp.ramp_link_prev)^2
    end
    if nlp.time != Final
        c += dot(nlp.λt, nlp.ramp_link_next)
        c += 0.5 * nlp.ρt * xnorm(nlp.ramp_link_next)^2
    end
    return c
end

function _adjoint_objective!(nlp::ProxALEvaluator, ∂f, pg)
    # Import model
    model = nlp.inner.model
    ## Seed left-hand side vector
    adjoint_cost!(model, ∂f, pg)
    ∂f .*= nlp.scale_objective
    if nlp.time != Origin
        ∂f .-= nlp.λf
        ∂f .-= nlp.ρf .* nlp.ramp_link_prev
    end
    if nlp.time != Final
        ∂f .+= nlp.λt
        ∂f .+= nlp.ρt .* nlp.ramp_link_next
    end
    ∂f .+= nlp.τ .* (pg .- nlp.pg_ref)
    return
end

## Gradient
function full_gradient!(nlp::ProxALEvaluator, jx, ju, w)
    # Import model
    model = nlp.inner.model
    # Import buffer (has been updated previously in update!)
    buffer = get(nlp.inner, PhysicalState())
    # Import AutoDiff objects
    ∂obj = nlp.inner.obj_stack
    # Scaling
    scale_obj = nlp.scale_objective
    # Current active power
    pg = get(nlp.inner, PS.ActivePower())

    u = @view w[1:nlp.nu]

    ## Compute adjoint of objective
    _adjoint_objective!(nlp, ∂obj.stack.∂pg, pg)

    ## Evaluate conjointly
    # ∇fₓ = v' * J,  with J = ∂pg / ∂x
    # ∇fᵤ = v' * J,  with J = ∂pg / ∂u
    put(model, PS.Generators(), PS.ActivePower(), ∂obj, buffer)
    copyto!(jx, ∂obj.stack.∇fₓ)
    copyto!(ju, ∂obj.stack.∇fᵤ)
end

function gradient_slack!(nlp::ProxALEvaluator, grad, w)
    # Gradient wrt s
    g_s = @view grad[nlp.nu+1:end]
    if nlp.time != Origin
        g_s .= nlp.λf .+ nlp.ρf .* nlp.ramp_link_prev
    else
        g_s .= 0.0
    end
end

function reduced_gradient!(nlp::ProxALEvaluator, grad, jvx, jvu, w)
    g = @view grad[1:nlp.nu]
    reduced_gradient!(nlp.inner, g, jvx, jvu, w)
    gradient_slack!(nlp, grad, w)
end

## Gradient
function gradient!(nlp::ProxALEvaluator, g, w)
    # Import AutoDiff objects
    ∂obj = nlp.inner.obj_stack

    jvu = ∂obj.stack.∇fᵤ ; jvx = ∂obj.stack.∇fₓ
    fill!(jvx, 0)  ; fill!(jvu, 0)
    full_gradient!(nlp, jvx, jvu, w)

    ## Evaluation of reduced gradient
    reduced_gradient!(nlp, g, jvx, jvu, w)
    return nothing
end

function constraint!(nlp::ProxALEvaluator, cons, w)
    u = @view w[1:nlp.nu]
    constraint!(nlp.inner, cons, u)
end

function jacobian_structure(nlp::ProxALEvaluator)
    m, nu = n_constraints(nlp), nlp.nu
    nnzj = m * nu
    rows = zeros(Int, nnzj)
    cols = zeros(Int, nnzj)
    jacobian_structure!(nlp.inner, rows, cols)
    return rows, cols
end

## Jacobian
function jacobian_structure!(nlp::ProxALEvaluator, rows, cols)
    jacobian_structure!(nlp.inner, rows, cols)
end

function jacobian!(nlp::ProxALEvaluator, jac, w)
    m = n_constraints(nlp)
    nnj = length(jac)
    u = @view w[1:nlp.nu]
    J = reshape(jac, m, div(nnj, m))
    Jᵤ = @view J[:, 1:nlp.nu]
    jacobian!(nlp.inner, Jᵤ, u)
end

function jprod!(nlp::ProxALEvaluator, jv, w, v)
    u = @view w[1:nlp.nu]
    vu = @view v[1:nlp.nu]
    jprod!(nlp.inner, jv, u, vu)
end

## Transpose Jacobian-vector product
## ProxAL does not add any constraint to the reduced model
function jtprod!(nlp::ProxALEvaluator, jv, w, v)
    u = @view w[1:nlp.nu]
    jvu = @view jv[1:nlp.nu]
    jtprod!(nlp.inner, jvu, u, v)
end
function full_jtprod!(nlp::ProxALEvaluator, jvx, jvu, w, v)
    u = @view w[1:nlp.nu]
    full_jtprod!(nlp.inner, jvx, jvu, u, v)
end

function ojtprod!(nlp::ProxALEvaluator, jv, u, σ, v)
    ∂obj = nlp.inner.obj_stack
    jvx = ∂obj.stack.jvₓ ; fill!(jvx, 0)
    jvu = ∂obj.stack.jvᵤ ; fill!(jvu, 0)
    # compute gradient of objective
    full_gradient!(nlp, jvx, jvu, u)
    jvx .*= σ
    jvu .*= σ
    # compute transpose Jacobian vector product of constraints
    full_jtprod!(nlp, jvx, jvu, u, v)
    # Evaluate gradient in reduced space
    reduced_gradient!(nlp, jv, jvx, jvu, u)
end

function hessprod!(nlp::ProxALEvaluator, hessvec, w, v)
    model = nlp.inner.model
    nx = get(model, NumberOfState())
    nu = get(model, NumberOfControl())
    buffer = get(nlp.inner, PhysicalState())
    H = nlp.inner.hessians

    u = @view w[1:nlp.nu]
    vᵤ = @view v[1:nlp.nu]
    vₛ = @view v[1+nlp.nu:end]

    fill!(hessvec, 0.0)

    hvu = @view hessvec[1:nlp.nu]
    hvs = @view hessvec[1+nlp.nu:end]

    z = H.z
    tgt_xu = H.tmp_tgt  # tangent wrt and u
    hv_xu = H.tmp_hv    # hv wrt x and u

    _second_order_adjoint_z!(nlp.inner, z, vᵤ)

    # Init tangent
    tgt_xu[1:nx] .= z
    tgt_xu[1+nx:nx+nu] .= vᵤ

    # Init adjoint
    ∂f = similar(buffer.pg)
    ∂²f = similar(buffer.pg)
    # Adjoint w.r.t. ProxAL's objective
    _adjoint_objective!(nlp, ∂f, buffer.pg)

    # Hessian of ProxAL's objective
    hessian_cost!(model, ∂²f)
    # Contribution of penalties
    ∂²f .+= nlp.τ
    if nlp.time != Origin
        ∂²f .+= nlp.ρf
    end
    if nlp.time != Final
        ∂²f .+= nlp.ρt
    end

    ## OBJECTIVE HESSIAN
    has_slack = (nlp.time != Origin)
    hessian_prod_objective_proxal!(model, H.obj, nlp.inner.obj_stack, hv_xu, hvs, ∂²f, ∂f, buffer, tgt_xu, vₛ, nlp.ρf, has_slack)
    ∇²fx = hv_xu[1:nx]       # / x
    ∇²fu = hv_xu[nx+1:nx+nu] # / u

    _reduced_hessian_prod!(nlp.inner, hvu, ∇²fx, ∇²fu, tgt_xu)
    return
end

## Utils function
function reset!(nlp::ProxALEvaluator)
    reset!(nlp.inner)
    # Reset multipliers
    fill!(nlp.λf, 0)
    fill!(nlp.λt, 0)
    # Reset proximal centers
    fill!(nlp.pg_f, 0)
    fill!(nlp.pg_ref, 0)
    fill!(nlp.pg_t, 0)
    # Reset buffers
    fill!(nlp.ramp_link_prev, 0)
    fill!(nlp.ramp_link_next, 0)
end

function primal_infeasibility!(nlp::ProxALEvaluator, cons, w)
    @assert length(w) == nlp.nu + nlp.ng
    u = @view w[1:nlp.nu]
    return primal_infeasibility!(nlp.inner, cons, u)
end
function primal_infeasibility(nlp::ProxALEvaluator, w)
    @assert length(w) == nlp.nu + nlp.ng
    u = @view w[1:nlp.nu]
    return primal_infeasibility(nlp.inner, u)
end


text:
# T , TSSM_HANDLE

println("including tssm_fourier.jl for type ",T)
for (METHOD, SUF, COMPLEX_METHOD, DIM) in (
                 (:Fourier1D, :_fourier_1d, true, 1 ),
                 (:Fourier2D, :_fourier_2d, true, 2 ),
                 (:Fourier3D, :_fourier_3d, true, 3 ),
                 (:FourierReal1D, :_fourier_real_1d, false, 1 ),
                 (:FourierReal2D, :_fourier_real_2d, false, 2 ),
                 (:FourierReal3D, :_fourier_real_3d, false, 3 ),
                )
println("    ", METHOD)

if T == :Float128
    SUF = string(SUF, "_wf128")
end
if DIM==1
@eval begin
    function ($METHOD)(T::Type{$T}, nx::Integer, xmin::Real, xmax::Real;
                       boundary_conditions::Integer=periodic)
        c = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"new",SUF))),
            Ptr{Nothing}, (Int32, ($T), ($T), Int32), nx, xmin, xmax, boundary_conditions)
        m = ($METHOD){$T}(c)
        finalizer(x -> ccall( Libdl.dlsym(($TSSM_HANDLE),
                       $(string(PRE,"finalize",SUF))), Nothing, (Ptr{Nothing},), x.m), m)
       m
    end
end # eval
if T == :Float64
@eval begin
    function ($METHOD)(nx::Integer, xmin::Real, xmax::Real;
                       boundary_conditions::Integer=periodic)
        ($METHOD)(($T), nx, xmin, xmax,  boundary_conditions=boundary_conditions)
    end
end # eval
end #if
elseif DIM==2
@eval begin
    function ($METHOD)(T::Type{$T}, nx::Integer, xmin::Real, xmax::Real,
                                       ny::Integer, ymin::Real, ymax::Real;
                       boundary_conditions::Integer=periodic)
        m = ($METHOD){$T}( ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"new",SUF))),
                        Ptr{Nothing}, (Int32, ($T), ($T), Int32, ($T), ($T), Int32),
                        nx, xmin, xmax, ny, ymin, ymax, boundary_conditions))
        finalizer(x -> ccall( Libdl.dlsym(($TSSM_HANDLE),
                       $(string(PRE,"finalize",SUF))), Nothing, (Ptr{Nothing},), x.m), m)
        m
    end
end # eval
if T == :Float64
@eval begin
    function ($METHOD)(nx::Integer, xmin::Real, xmax::Real,
                       ny::Integer, ymin::Real, ymax::Real;
                       boundary_conditions::Integer=periodic)
        ($METHOD)(($T), nx, xmin, xmax,  ny, ymin, ymax,
                  boundary_conditions=boundary_conditions)
    end
end # eval
end #if
elseif DIM==3
@eval begin
    function ($METHOD)(T::Type{$T}, nx::Integer, xmin::Real, xmax::Real,
                                       ny::Integer, ymin::Real, ymax::Real,
                                       nz::Integer, zmin::Real, zmax::Real;
                       boundary_conditions::Integer=periodic)
        m = ($METHOD){$T}( ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"new",SUF))),
                        Ptr{Nothing}, (Int32, ($T), ($T), Int32, ($T), ($T),
                        Int32, ($T), ($T), Int32),
                        nx, xmin, xmax, ny, ymin, ymax, nz, zmin, zmax,
                        boundary_conditions))
        finalizer(x -> ccall( Libdl.dlsym(($TSSM_HANDLE),
                       $(string(PRE,"finalize",SUF))), Nothing, (Ptr{Nothing},), x.m), m)
        m
    end
end # eval
if T == :Float64
@eval begin
    function ($METHOD)(nx::Integer, xmin::Real, xmax::Real,
                       ny::Integer, ymin::Real, ymax::Real,
                       nz::Integer, zmin::Real, zmax::Real;
                       boundary_conditions::Integer=periodic)
        ($METHOD)(($T), nx, xmin, xmax, ny, ymin, ymax, nz, zmin, zmax,
                  boundary_conditions=boundary_conditions)
    end
end # eval
end #if

end # if

@eval begin
    function get_nx(m::($METHOD){$T})
       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_nx",SUF))), Int32,
             (Ptr{Nothing}, ), m.m )
    end

    function get_xmin(m::($METHOD){$T})
       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_xmin",SUF))), ($T),
             (Ptr{Nothing}, ), m.m )
    end

    function get_xmax(m::($METHOD){$T})
       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_xmax",SUF))), $T,
             (Ptr{Nothing}, ), m.m )
    end
end # eval

if DIM>=2
    @eval begin

        function get_ny(m::($METHOD){$T})
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_ny",SUF))), Int32,
                 (Ptr{Nothing}, ), m.m )
        end

        function get_ymin(m::($METHOD){$T})
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_ymin",SUF))), ($T),
                 (Ptr{Nothing}, ), m.m )
        end

        function get_ymax(m::($METHOD){$T})
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_ymax",SUF))), ($T),
                 (Ptr{Nothing}, ), m.m )
        end
    end # eval
end

if DIM>=3
    @eval begin
        function get_nz(m::($METHOD){$T})
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_nz",SUF))), Int32,
                 (Ptr{Nothing}, ), m.m )
        end

        function get_zmin(m::($METHOD){$T})
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_zmin",SUF))), ($T),
                 (Ptr{Nothing}, ), m.m )
        end

        function get_zmax(m::($METHOD){$T})
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_zmax",SUF))), ($T),
                 (Ptr{Nothing}, ), m.m )
        end
    end # eval
end

end # for


text:
macro grab_output(ex)
    quote
        mktemp() do fname, fout
            ret = nothing
            open(fname, "w") do fout
                redirect_stdout(fout) do
                    ret = $(esc(ex))
                end
            end
            ret, read(fname, String)
        end
    end
end

# some tests are mysteriously broken on CI, but nowhere else...
# use a horrible macro to mark those tests as "potentially broken"
@eval Test begin
    export @test_maybe_broken

    const broken = haskey(ENV, "CI")
    macro test_maybe_broken(ex, kws...)
        test_expr!("@test_maybe_broken", ex, kws...)
        orig_ex = Expr(:inert, ex)
        result = get_test_result(ex, __source__)
        if !broken
            :(do_test($result, $orig_ex))
        else
            quote
                x = $result
                if x.value
                    do_test(x, $orig_ex)
                else
                    do_broken_test(x, $orig_ex)
                end
            end
        end
    end
end


text:
using JuliaInterpreter
using Test

@testset "core" begin
    @test JuliaInterpreter.is_quoted_type(QuoteNode(Int32), :Int32)
    @test !JuliaInterpreter.is_quoted_type(QuoteNode(Int32), :Int64)
    @test !JuliaInterpreter.is_quoted_type(QuoteNode(Int32(0)), :Int32)
    @test !JuliaInterpreter.is_quoted_type(Int32, :Int32)
end


text:
module CommunalHelperTrackSwitchTrigger

using ..Ahorn, Maple

modes = ["Alternate", "On", "Off"]

@mapdef Trigger "CommunalHelper/TrackSwitchTrigger" TrackSwitchTrigger(
    x::Integer,
    y::Integer,
    width::Integer=Maple.defaultTriggerWidth,
    height::Integer=Maple.defaultTriggerHeight,
    oneUse::Bool=true,
    flash::Bool=false,
    globalSwitch::Bool=false,
    mode::String="Alternate",
)

const placements = Ahorn.PlacementDict(
    "Track Switch Trigger (Communal Helper)" => Ahorn.EntityPlacement(
        TrackSwitchTrigger,
        "rectangle"
    ),
)

Ahorn.editingOptions(trigger::TrackSwitchTrigger) = Dict{String,Any}(
    "mode" => modes,
)

end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(Tuple{IonBase.SearchCmd.var"#fetch_repo##kw", NamedTuple{(:auth,), Tuple{GitHub.OAuth2}}, typeof(IonBase.SearchCmd.fetch_repo), Pkg.Types.RegistrySpec, Base.Dict{String, Any}})
    precompile(Tuple{IonBase.SearchCmd.var"#fetch_repo_from_url##kw", NamedTuple{(:auth,), Tuple{GitHub.OAuth2}}, typeof(IonBase.SearchCmd.fetch_repo_from_url), String})
    precompile(Tuple{typeof(IonBase.CreateCmd.copy_templates), String})
    precompile(Tuple{typeof(IonBase.CreateCmd.read_template), String})
    precompile(Tuple{typeof(IonBase.Doc.build), String})
    precompile(Tuple{typeof(IonBase.InstallCmd.find_julia_installer_info), String})
    precompile(Tuple{typeof(IonBase.ReleaseCmd.is_version_number), String})
    precompile(Tuple{typeof(IonBase.ReleaseCmd.update_version!), IonBase.ReleaseCmd.Project, String})
    precompile(Tuple{typeof(IonBase.SearchCmd.print_stars), Base.TTY, GitHub.Repo})
    precompile(Tuple{typeof(IonBase.SearchCmd.search_exact_package), String})
    precompile(Tuple{typeof(IonBase.SearchCmd.search_fuzzy_package), String})
    precompile(Tuple{typeof(IonBase.dot_ion), String, String})
    precompile(Tuple{typeof(IonBase.search), String})
    precompile(Tuple{typeof(IonBase.templates), String, Int})
    precompile(Tuple{typeof(IonBase.templates), String, String})
end


text:
module PyFrom
    using Base: Symbol, Tuple
    using PyCall
    export @pyfrom

    struct ImportPhraseException <: Exception
        var::String
    end

    struct DotExpansionException <: Exception
        var::String
    end

    macro pyfrom(modname, importphrase)
        return pyfrom(modname, importphrase)
    end

    pyfrom(modname, importphrase) = begin
        errmessage() = begin
            ("@pyfrom must follow a pattern like `@pyfrom a.b import c as d, e, f as g`, " *
            "got `@pyfrom $modname $importphrase`")
        end

        nodes, paths = try
            dotpath_expansion(modname)
        catch e
            e isa DotExpansionException  && throw(DotExpansionException(errmessage()))
            rethrow(e)
        end

        importmap = try
            importphrase_to_mapping(importphrase)
        catch e
            e isa ImportPhraseException  && throw(ImportPhraseException(errmessage()))
            rethrow(e)
        end

        # The secret souce
        names = [name for (_, name) in importmap]
        esc(:(($(Expr(:tuple, names...)) = $get_module_imports($nodes, $paths, $importmap))[end]))
    end

    get_module_imports(nodes::Vector{Symbol},
                       paths::Vector{Union{Symbol, Expr}},
                       importmap::Vector{Pair{Symbol, Symbol}}) = begin
        # import from deeper nest, i.e. a -> a.b -> a.b.c -> a.b.c.d
        nextleaf(moduleleaf, node, path) = begin
            if moduleleaf === nothing || !hasproperty(moduleleaf, node)
                PyCall._pywrap_pyimport(pyimport(PyCall.modulename(path)))
            else
                getproperty(moduleleaf, node)
            end
        end

        # find module leaf
        moduleleaf = nothing
        for (node, path) in zip(nodes, paths)
            moduleleaf = nextleaf(moduleleaf, node, path)
        end

        # get relevant children from importmap in a list of modules/pyobjects
        imports = []
        for (node, _) in importmap
            path = Expr(:., paths[end], QuoteNode(node))
            push!(imports, nextleaf(moduleleaf, node, path))
        end

        return imports
    end

    dotpath_expansion(expr) = begin
        errmessage = "Expected dotted expression like `a.b.c.d.e`, got $(expr)"

        nodes = Vector{Symbol}()
        paths = Vector{Union{Symbol, Expr}}()
        ex = expr
        while ex isa Expr && ex.head == (:.)
            push!(paths, ex)

            ex, tail = ex.args[1], ex.args[end]
            if tail isa QuoteNode
                tail = tail.value
            end

            if tail isa Symbol
                push!(nodes, tail)
            else
                throw(DotExpansionException(errmessage))
            end
        end
        if ex isa Symbol
            push!(paths, ex)
            push!(nodes, ex)
        else
            throw(DotExpansionException(errmessage))
        end

        nodes = reverse(nodes)
        paths = reverse(paths)

        return nodes, paths
    end

    importphrase_to_mapping(expr::Expr)::Vector{Pair{Symbol, Symbol}} = begin
        err_message = "Expected import statement pattern like `import a as b, c, d as e, f, g`, got `$(expr)`"

        if expr.head != :import
            throw(ImportPhraseException(err_message))
        end

        # group `a as b, c as d, ...` as [[:a, :b], [:c, :d], ...]
        pairs = []
        for i in expr.args
            push!(
                pairs,
                if i.head == :as
                    [i.args[1], i.args[2]]
                else
                    length(i.args) != 1 && throw(ImportPhraseException(err_message))
                    [i.args[1], i.args[1]]
                end
            )
        end

        # sanitize inputs and convert to [:a=>:b, :c=>d, ...] type of list
        pairs_final = Vector{Pair{Symbol, Symbol}}()
        for pair in pairs
            for (j, item) in enumerate(pair)
                if item isa Symbol
                elseif item isa Expr
                    if item.head == :. && length(item.args) == 1 && item.args[1] isa Symbol
                        pair[j] = item.args[1]
                    else
                        throw(ImportPhraseException(err_message))
                    end
                else
                    throw(ImportPhraseException(err_message))
                end
            end
            push!(pairs_final, pair[1] => pair[2])
        end

        return pairs_final
    end

end


text:
include("utils.jl")

build_so("007")

type Foo
   i::Cint
end

x = Foo(13)

y = ccall(
    (:bar, "binary/lib007.so"),
    Void,
    (Ref{Foo}, ),
    Ref(x)
)

@assert x.i == 42

@assert y === nothing


text:
# Conservative 7-color palette
# Wong, Bang. "Points of view: Color blindness." (2011): 441.
# https://www.nature.com/articles/nmeth.1618?WT.ec_id=NMETH-201106

function wongcolors()
    return [
        RGB(0/255, 114/255, 178/255), # blue
        RGB(230/255, 159/255, 0/255), # orange
        RGB(0/255, 158/255, 115/255), # green
        RGB(204/255, 121/255, 167/255), # reddish purple
        RGB(86/255, 180/255, 233/255), # sky blue
        RGB(213/255, 94/255, 0/255), # vermillion
        RGB(240/255, 228/255, 66/255), # yellow
    ]
end

const font_folder = RelocatableFolders.@path joinpath(dirname(@__DIR__), "assets", "fonts")

firasans(weight) = joinpath(font_folder, "FiraSans-$(weight).ttf")
opensans(weight) = joinpath(font_folder, "OpenSans-$(weight).ttf")

# Batlow colormap
# Crameri, Fabio, Grace E. Shephard, and Philip J. Heron. "The misuse of colour in science communication." Nature communications 11.1 (2020): 1-10.
# https://www.nature.com/articles/s41467-020-19160-7

function aog_theme(; fonts=[firasans("Medium"), firasans("Light")])
    mediumfont = first(fonts)
    lightfont = last(fonts)

    marker = :circle

    colormap = :batlow
    linecolor = :gray25
    markercolor = :gray25
    patchcolor = :gray25

    palette = (
        color=wongcolors(),
        patchcolor=wongcolors(),
        marker=[:circle, :utriangle, :cross, :rect, :diamond, :dtriangle, :pentagon, :xcross],
        linestyle=[:solid, :dash, :dot, :dashdot, :dashdotdot],
        side=[:left, :right],
    )

    # setting marker here is a temporary hack
    # it should either respect `marker = :circle` globally
    # or `:circle` and `Circle` should have the same size
    BoxPlot = (mediancolor=:white, marker=:circle)
    Scatter = (marker=:circle,)
    Violin = (mediancolor=:white,)

    Axis = (
        xgridvisible=false,
        ygridvisible=false,
        topspinevisible=false,
        rightspinevisible=false,
        bottomspinecolor=:darkgray,
        leftspinecolor=:darkgray,
        xtickcolor=:darkgray,
        ytickcolor=:darkgray,
        xticklabelfont=lightfont,
        yticklabelfont=lightfont,
        xlabelfont=mediumfont,
        ylabelfont=mediumfont,
        titlefont=mediumfont,
    )
    Axis3 = (
        protrusions=55, # to include label on z axis, should be fixed in Makie
        xgridvisible=false,
        ygridvisible=false,
        zgridvisible=false,
        xspinecolor_1=:darkgray,
        yspinecolor_1=:darkgray,
        zspinecolor_1=:darkgray,
        xspinecolor_2=:transparent,
        yspinecolor_2=:transparent,
        zspinecolor_2=:transparent,
        xspinecolor_3=:transparent,
        yspinecolor_3=:transparent,
        zspinecolor_3=:transparent,
        xtickcolor=:darkgray,
        ytickcolor=:darkgray,
        ztickcolor=:darkgray,
        xticklabelfont=lightfont,
        yticklabelfont=lightfont,
        zticklabelfont=lightfont,
        xlabelfont=mediumfont,
        ylabelfont=mediumfont,
        zlabelfont=mediumfont,
        titlefont=mediumfont,
    )
    Legend = (
        framevisible=false,
        gridshalign=:left,
        padding=(0f0, 0f0, 0f0, 0f0),
        labelfont=lightfont,
        titlefont=mediumfont,
    )
    Colorbar = (
        flip_vertical_label=true,
        spinewidth=0,
        ticklabelfont=lightfont,
        labelfont=mediumfont,
    )

    return (;
        marker,
        colormap,
        linecolor,
        markercolor,
        patchcolor,
        palette,
        BoxPlot,
        Scatter,
        Violin,
        Axis,
        Axis3,
        Legend,
        Colorbar,
    )
end

function set_aog_theme!(; kwargs...)
    theme = aog_theme(; kwargs...)
    return set_theme!(; theme...)
end


text:
module NeuroCore

using AxisIndices
using NamedDims
using TimeAxes
using FieldProperties
using StaticArrays
using GeometryBasics
using Mmap
using CoordinateTransformations
using LinearAlgebra
using Unitful
using Unitful: s, Hz, T, °, mm

using CoordinateTransformations: Rotations
using CoordinateTransformations.Rotations: SPQuat, Rotation

using Base: @propagate_inbounds

using Reexport

include("./ColorChannels/ColorChannels.jl")

include("./SpatialAPI/SpatialAPI.jl")
@reexport using .SpatialAPI

include("./AnatomicalAPI/AnatomicalAPI.jl")
using .AnatomicalAPI
@reexport using .AnatomicalAPI

#include("./Nodes/Nodes.jl")
#using .Nodes
#@reexport using .Nodes

include("./swapstream.jl")
using .SwapStreams
@reexport using .SwapStreams

include("./NeuroMetadata/NeuroMetadata.jl")
using .NeuroMetadata

const F64Second = typeof(one(Float64) * s)
const F64Tesla = typeof(one(Float64) * T)
const F64kOhms = typeof(1.0u"kΩ")
const F64Hertz = typeof(one(Float64) * Hz)
const IntDegree = typeof(1 * °)

include("io.jl")

"""
    CoordinateSpace

Returns an instance of `CoordinateSystem` describing the coordinate system for `x`.
"""
struct CoordinateSpace{S}
    space::S
end

const UnkownCoordinatesSpace = CoordinateSpace(nothing)

coordinate_space(x) = UnkownCoordinatesSystem

#= TODO formalize this interaction
CoordinateSpace(sc::MetadataArray{T,N,<:CoordinateSystem}) = metadata(sc)

CoordinateSpace(sc::AbstractAxisArray) = (parent(sc))

CoordinateSpace(sc::NamedDimsArray) = (parent(sc))
=#

# TODO document SpatialCoordinates
"""
    SpatialCoordinates

"""
const SpatialCoordinates{L,CS,N,Axs} = NamedMetaCartesianAxes{L,N,CoordinateSpace{CS},Axs}

function SpatialCoordinates(x)
    return NamedMetaCartesianAxes{spatial_order(x)}(spatial_axes(x), metadata=CoordinateSpace(x))
end

#= TODO Is this a good name for this
"The anatomical coordinate system."
@defprop AnatomicalSystem{:anatsystem}

"The acquisition coordinate system."
@defprop AcquisitionSystem{:acqsystem}
=#

end


text:
module Stream

export Map, TakeWhile, SplitBy, Enumerate

struct Map
    closure
    kwparams
    buffer_size
    Map(closure; buffer_size = 0,  kwparams...) = new(closure,kwparams, buffer_size)
end

function (m::Map)(iterable  )
    Channel(m.buffer_size;m.kwparams...) do c
        for var in iterable
            push!(c, m.closure(var))
        end
    end
end

struct TakeWhile
    predicate
    kwparams
    buffer_size
    TakeWhile(predicate; buffer_size = 0,  kwparams...) = new(predicate,kwparams, buffer_size )
end

function (tw::TakeWhile)(channel_like )
    Channel(tw.buffer_size; tw.kwparams...) do c
        try
            while tw.predicate(fetch(channel_like))
                push!(c,take!(channel_like))
            end
        catch e
            if isa(e, InvalidStateException) && e.state === :closed
                return
            else
                rethrow()
            end
        end
    end
end

struct SplitBy
    predicate
    kwparams
    buffer_size
    keepend
    keepempty
    SplitBy(predicate; keepend=false, keepempty=false, buffer_size = 0,  kwparams...) = new(predicate,kwparams, buffer_size, keepend, keepempty)
end

function (sb::SplitBy)(iterable)
    Channel(sb.buffer_size; sb.kwparams...) do c
        buffer = []
        for msg in iterable
            if sb.predicate(msg)
                if sb.keepend push!(buffer,msg) end
                if !isempty(buffer) || sb.keepempty
                    push!(c,buffer)
                    buffer = []
                end
            else
                push!(buffer,msg)
            end
        end
        if !isempty(buffer) push!(c,buffer) end
    end
end

struct Enumerate
    from::Int
    buffer_size
    kwparams
    Enumerate(from::Int=1 ;buffer_size=0, kwparams...) = new(from,buffer_size,kwparams)
end

function (en::Enumerate)(iterable)
    Channel() do c
        counter = Int(en.from)
        for item in iterable
            push!(c,(counter,item))
            counter += 1
        end
    end
end

end


text:
using Documenter, AdvancedTopics

makedocs(;
    modules=[AdvancedTopics],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/rjww/AdvancedTopics.jl/blob/{commit}{path}#L{line}",
    sitename="AdvancedTopics.jl",
    authors="Robert Woods",
    assets=String[],
)

deploydocs(;
    repo="github.com/rjww/AdvancedTopics.jl",
)


text:
#
# Date created: 2022-03-10
# Author: aradclif
#
#
############################################################################################

function sumsqbody(N::Int, D)
    body = Expr(:block)
    params = D.parameters
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    b = Expr(:ref, :B, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    δ = Expr(:(=), :Δ, a)
    e = Expr(:(=), b, Expr(:call, :+, b, Expr(:call, :*, :Δ, :Δ)))
    push!(body.args, δ)
    push!(body.args, e)
    body
end

function sumsq_quote(N::Int, D)
    ls = loopgen(N)
    body = sumsqbody(N, D)
    push!(ls.args, body)
    return quote
        @turbo $ls
        return B
    end
end

@generated function sumsq!(B::AbstractArray{T, N}, A::AbstractArray{T, N}, dims::D) where {T, N, D}
    sumsq_quote(N, D)
end

function sqrt!(A::AbstractArray{T, N}) where {T, N}
    @turbo for i ∈ eachindex(A)
        A[i] = sqrt(A[i])
    end
    A
end

function lveuclidean(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    Dᴬ = size(A)
    Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : Dᴬ[d], N)
    B = zeros(Base.promote_op(/, T, Int), Dᴮ′)
    sumsq!(B, A, Dᴮ′)
    sqrt!(B)
    B
end