text:
losses = Float64[10, 8, 9, 10, 11, 12, 12, 13, 14, 15, 16, 17, 16]

# codecov:
@test EarlyStopping._min(nothing, 5) == 5

@testset "Never" begin
    @test stopping_time(Never(), losses) == 0
end

@testset "InvalidValue" begin
    @test stopping_time(InvalidValue(), losses) == 0
    N = 5
    losses2 = fill(123.4, N)
    @test all(reverse(eachindex(losses2))) do j
        losses2[j] = NaN
        stopping_time(InvalidValue(), losses2) == j
    end

    is_training = map(x -> x%3 > 0, 1:length(losses))
    @test stopping_time(InvalidValue(), losses, is_training) == 0
    for n = 1:2:length(losses)
        n_stop = sum(!, is_training[1:n])
        losses2 = copy(losses)
        losses2[n] = Inf
        @test stopping_time(InvalidValue(), losses2, is_training) == n_stop
        losses2[n] = NaN
        @test stopping_time(InvalidValue(), losses2, is_training) == n_stop
    end
end

struct SleepyIterator{T}
    iter::T
    t::Float64
end
SleepyIterator(iter; t=0.1) = SleepyIterator(iter, t)

Base.iterate(iter::SleepyIterator) = (sleep(iter.t); iterate(iter.iter))
Base.iterate(iter::SleepyIterator, state) =
    (sleep(iter.t); iterate(iter.iter, state))

@testset "TimeLimit" begin
    @test_throws ArgumentError TimeLimit(t=0)
    @test TimeLimit(1).t == Millisecond(3_600_000)
    @test TimeLimit(t=Day(2)).t == Millisecond(48*3_600_000)
    sleepy_losses = SleepyIterator(losses; t=0.1)
    @test stopping_time(TimeLimit(t=Millisecond(600)), sleepy_losses) == 7
    # codecov:
    @test EarlyStopping.update_training(TimeLimit(), 42.0) <= now()
end

@testset "GL" begin
    # constructor:
    @test_throws ArgumentError GL(alpha=0)
    @test GL(alpha=1).alpha === 1.0

    # stopping times:
    n = @test_logs((:info, r"loss updates: 1"),
                   (:info, r"state: \(loss = 10.0, min_loss = 10.0\)"),
                   (:info, r"loss updates: 2"),
                   (:info, r"state: \(loss = 8.0, min_loss = 8.0\)"),
                   (:info, r"loss updates: 3"),
                   (:info, r"state: \(loss = 9.0, min_loss = 8.0\)"),
                   stopping_time(GL(alpha=12), losses, verbosity=1))
    @test n == 3
    @test stopping_time(GL(alpha=20), losses) == 4
    @test stopping_time(GL(alpha=40), losses) == 6
    @test stopping_time(GL(alpha=90), losses) == 11
    @test stopping_time(GL(alpha=110), losses) == 12
    @test stopping_time(GL(alpha=1000), losses) == 0
end

@testset "PQ" begin
    v = [:c, :b, :a]
    @test EarlyStopping.prepend(v, :d, 4) == [:d, :c, :b, :a]
    @test EarlyStopping.prepend(v, :d, 3) == [:d, :c, :b]
    @test EarlyStopping.progress(Float64[2, 3, 6, 7]) ≈ 1000*(4.5 - 2)/2
    @test_throws ArgumentError PQ(alpha=0)
    @test_throws ArgumentError PQ(k=-Inf)
    @test_throws ArgumentError PQ(k=1.2)
    @test_throws ArgumentError PQ(k=1)

    c = PQ(alpha=10, k=2)

    state = EarlyStopping.update_training(c, 10.0)
    state = EarlyStopping.update_training(c, 10.0, state)
    state = EarlyStopping.update(c, 10.0, state)
    @test EarlyStopping.done(c, state) # progress = 0

    state = EarlyStopping.update_training(c, 10.0, state)
    # can't be done if last update was a training update:
    @test !EarlyStopping.done(c, state)

    #                 k=2                progress GL    PQ    t
    losses2 = [9.5, 9.3, 10,            # 10.8     0     0     1
              9.3, 9.1, 8.9, 8,        # 11.2     0     0     2
              8.3, 8.4, 9,             # 6.02     12.5  2.08  3
              9.9, 9.5, 10,            # 21.2     25.0  1.18  4
              10.6, 10.4, 11,          # 9.61     37.5  3.90  5
              8888, 11.8, 11.7, 12,    # 4.27     50.0  11.7  6
              11.6, 11.4, 12,          # 8.77     50.0  5.70  7
              12.2, 12.1, 13,          # 4.1      62.5  15.2  8
              14.5, 14.1, 14,          # 14.2     75    5.28  9
              13.9, 13.7, 15,          # 7.30     87.5  11.9  10
              12.5, 12.3, 16,          # 8.13     100   12.3  11
              11.2, 11.0, 17,          # 9.09     112.5 12.4  12
              10.5, 10.3, 16]          # 9.71     100   10.3  13
    is_training = Bool[1, 1, 0,
                       1, 1, 1, 0,
                       1, 1, 0,
                       1, 1, 0,
                       1, 1, 0,
                       1, 1, 1, 0,
                       1, 1, 0,
                       1, 1, 0,
                       1, 1, 0,
                       1, 1, 0,
                       1, 1, 0,
                       1, 1, 0,
                       1, 1, 0]

    @test stopping_time(PQ(alpha=2.0, k=2), losses2, is_training) == 3
    @test stopping_time(PQ(alpha=3.8, k=2), losses2, is_training) == 5
    @test stopping_time(PQ(alpha=11.6, k=2), losses2, is_training) == 6
    @test stopping_time(PQ(alpha=15.1, k=2), losses2, is_training) == 8
    @test stopping_time(PQ(alpha=15.3, k=2), losses2, is_training) == 0
end

@testset "Patience" begin
    @test_throws ArgumentError Patience(n=0)
    @test stopping_time(Patience(n=6), losses) == 0
    @test stopping_time(Patience(n=5), losses) == 12
    @test stopping_time(Patience(n=4), losses) == 6
    @test stopping_time(Patience(n=3), losses) == 5
    @test stopping_time(Patience(n=2), losses) == 4
    @test stopping_time(Patience(n=1), losses) == 3
end

@testset "NumberSinceBest" begin
    @test_throws ArgumentError NumberSinceBest(n=0)
    @test stopping_time(NumberSinceBest(n=6), losses) == 8
    @test stopping_time(NumberSinceBest(n=5), losses) == 7
    @test stopping_time(NumberSinceBest(n=4), losses) == 6
    @test stopping_time(NumberSinceBest(n=3), losses) == 5
    @test stopping_time(NumberSinceBest(n=2), losses) == 4
    @test stopping_time(NumberSinceBest(n=1), losses) == 3

    losses2 = Float64[10, 9, 8, 9, 10, 7, 10, 10, 10, 10]
    @test stopping_time(NumberSinceBest(n=2), losses2) == 5
    @test stopping_time(NumberSinceBest(n=3), losses2) == 9
end

@testset "NumberLimit" begin
    @test_throws ArgumentError NumberLimit(n=0)

    for i in 1:length(losses)
        @test stopping_time(NumberLimit(i), losses) == i
    end
end

@testset "Threshold" begin
    @test Threshold().value == 0.0
    stopping_time(Threshold(2.5), Float64[12, 32, 3, 2, 5, 7]) == 4
end

@testset "robustness to first loss being a training loss" begin
    criteria = filter(subtypes(StoppingCriterion)) do C
        C != NotANumber # as deprecated
    end
    for C in criteria
        losses = float.(4:-1:1)
        is_training = [true, true, false, false]
        stopping_time(C(), losses, is_training)
    end
end

@testset "Warmup" begin
    @test_throws ArgumentError Warmup(Patience(), 0)
    for n in 1:(length(losses)-1)
        @test stopping_time(Warmup(NumberLimit(1), n), losses) == n+1
    end

    # Test message
    @testset "message" begin
        stopper = Warmup(Patience(2); n = 2)
        stopper_ref = Warmup(Patience(2), 2)
        state, state_ref = nothing, nothing
        for loss = losses
            state = EarlyStopping.update(stopper, loss, state)
            state_ref = EarlyStopping.update(stopper_ref, loss, state_ref)
            @test message(stopper, state) == message(stopper_ref, state_ref)
        end
    end

    @testset "training" begin
        stopper = Warmup(PQ(), 3)
        is_training = map(x -> x%3 > 0, 1:length(losses))

        # Feed 2 training losses + 1 non-training to criteria with/without
        stop_time = stopping_time(stopper, losses, is_training)
        ref_stop_time = stopping_time(stopper.criterion, losses[3:end], is_training)

        # PQ only counts training loss updates
        @test round(stop_time/3, RoundUp) == ref_stop_time
    end

    @testset "integration" begin
        @test_criteria Warmup(Patience())
        @test_criteria Warmup(NumberSinceBest())
        @test_criteria Warmup(Patience(3) + InvalidValue())
    end
end

# # DEPRECATED

@testset "NotANumber" begin
    @test_deprecated criterion = NotANumber()
    @test stopping_time(criterion, losses) == 0
    N = 5
    losses2 = fill(123.4, N)
    @test all(reverse(eachindex(losses2))) do j
        losses2[j] = NaN
        stopping_time(criterion, losses2) == j
    end
    is_training = map(x -> x%3 > 0, 1:length(losses))
    @test stopping_time(criterion, losses, is_training) == 0
    for n = 1:2:length(losses)
        losses2 = copy(losses); losses2[n] = NaN
        @test stopping_time(criterion, losses2, is_training) == sum(!, is_training[1:n])
    end
end

true


text:
import TrajectoryOptimization: dynamics
using RigidBodyDynamics

f = Dynamics.dubins_dynamics!
n,m = 3,2
model = Model(f,n,m)
@test model.m == m
@test model.n == n

xdot = zeros(n)
x = rand(n)
u = rand(m)
@test evals(model) == 0
reset(model)
model.f(xdot,x,u)
@test evals(model) == 0
dynamics(model,xdot,x,u)
@test evals(model) == 1

acrobot = parse_urdf(Dynamics.urdf_doublependulum)
model = Model(acrobot)
@test evals(model) == 0
n,m = model.n, model.m
xdot = zeros(n)
x = rand(n)
u = rand(m)
dynamics(model,xdot,x,u)
@test evals(model) == 1
reset(model)
@test evals(model) == 0

ext: jl


text:
using MOCNeutronTransport
using Plots
pyplot()
T = Float64
Nsoln = 20
Nset = 1:10
f = in

p₁ = Point_2D(T, 0)
p₂ = Point_2D(T, 2)
p₃ = Point_2D(T, 2, 2)
p₄ = Point_2D(T, 3//2, 1//4)
p₅ = Point_2D(T, 3, 1)
p₆ = Point_2D(T, 1, 1)
tri6 = Triangle6_2D((p₁, p₂, p₃, p₄, p₅, p₆))
p = Point_2D(T, 0, 0)

soln = f(p, tri6; N=Nsoln)
err = [f(p, tri6; N = n) for n = Nset]
plot(Nset[1:length(Nset)-1],
     err[1:length(Nset)-1],
     title = "Newton-Raphson iterations vs Triangle6_2D ∈  error (Nsoln = $Nsoln)")
xlabel!("Newton-Raphson iterations")
ylabel!("∈  error")

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
# Written in Julia by André Ferrari
# April 2015
#
# Translation in julia of algorithm in section 2.2:
# "Generating Stationary Processes via Circulant embedding"

# Kroese, D.P. and Botev, Z.I. (2013).
# "Spatial Process Generation."
# Lectures on Stochastic Geometry, Spatial Statistics and Random Fields,
# Volume II: Analysis, Modeling and Simulation of Complex Structures,
# V. Schmidt (Ed.).  Springer-Verlag, Berlin.
#
# http://arxiv.org/abs/1308.0399
#

module GHF

export genghf

function genghf(tx::Vector, ty::Vector, rho::Function)
    m = length(tx)
    n = length(ty)
    rows = zeros(m,n)
    cols = zeros(m,n)

    for i = 1:n, j = 1:m
        rows[j,i] = rho(tx[i]-tx[1],ty[j]-ty[1])
        cols[j,i] = rho(tx[1]-tx[i],ty[j]-ty[1])
    end

    blkcirc_row = [[rows cols[:,end:-1:2]] ;
        [cols[end:-1:2,:] rows[end:-1:2,end:-1:2]]]

    lam = real(fft(blkcirc_row))/(2*m-1)/(2*n-1)

    ind0 = lam .< 0.0
    if (isempty(ind0)) && (abs(minimum(lam[ind0])) > 1e-15)
        error("Could not find positive definite embedding!")
    else
        lam[ind0] = 0.0
        lam = sqrt.(lam)
    end

    z = complex.(randn(2*m-1, 2*n-1), randn(2*m-1, 2*n-1))
    f = fft(lam .* z)
    f = f[1:m,1:n]
    field1 = real(f)
#   field2 = imag(f)  # if needed

end
end

max_stars_repo_name: andferrari/GHL

max_issues_repo_name: andferrari/GHL

max_forks_repo_name: andferrari/GHL


text:
using CSV
using Plots

function plot_gpu_temps(filename)
  data = CSV.File(filename, comment="#", delim=" ", ignorerepeated=true)

  x = size(data)[1]
  y = hcat(data.gtemp, data.pwr)

  p = plot(x, y,
    title = "GPU temperatures and power consumptions during idle and stress test",
    xlabel = "Seconds after start",
    ylabel = "Temp (C) and Power (W)",
    size = (1200, 800),
    label = ["Temperature", "Power"],
    )

  show(p)
  #display(p)
end


text:
# this file tests various aspects of the advection routines
using MagmaThermoKinematics
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
using Plots
using LinearAlgebra
using SpecialFunctions
using Test

#using WriteVTK

const CreatePlots = false      # easy way to deactivate plotting throughout

function test_HostRockVelocityFromDike(Dimension="2D", DikeType="ElasticDike", DikeAngle=[45])
  # test generating host velocity from various dikes, with different size/orientation/type in both 2D and 3DD

  if Dimension=="2D"
    # Model parameters
    W,H                     =   30.0,  30.0;                                # Width, Length, Height

    # Define grid
    Nx, Nz                  =   129, 129;                                     # resolution of coarse grid
    dx,dz                   =   W*1e3/(Nx-1), H*1e3/(Nz-1);                   # grid size [m]
    x,z                     =   0:dx:W*1e3, -H*1e3:dz:0;                      # 1D coordinate arrays
    coords                  =   collect(Iterators.product(x,z))               # generate coordinates from 1D coordinate vectors
    X,Z                     =   (x->x[1]).(coords), (x->x[2]).(coords);       # transfer coords to 3D arrays
    Grid, FullGrid, Spacing =   (x,z), (X,Z), (dx,dz);

    Hdike                   =   100.0;
    Wdike                   =   20000.0;
    T_in                    =   900.0;

    cen                     =   [W/2;-H/2].*1e3;
  elseif Dimension=="3D"
      # Model parameters
      W,L,H                 =   30., 40., 50.;                                    # Width, Length, Height

      # Define coarse grid
      Nx, Ny, Nz              =   65,65,65;                                                    # resolution of coarse grid
      dx,dy,dz                =   W*1e3/(Nx-1), L*1e3/(Ny-1), H*1e3/(Nz-1);                     # grid size [m]
      x,y,z                   =   0:dx:((Nx-1)*dx),  0:dy:((Ny-1)*dy), -((Nz-1)*dz):dz:0.;      # 1D coordinate arrays
      coords                  =   collect(Iterators.product(x,y,z))                             # generate coordinates from 1D coordinate vectors
      X,Y,Z                   =   (x->x[1]).(coords), (x->x[2]).(coords), (x->x[3]).(coords);   # transfer coords to 3D arrays
      Grid, FullGrid, Spacing =   (x,y,z), (X,Y,Z), (dx,dy,dz);
      cen                     =   [W/2;L/2; -H/2].*1e3;

      Hdike                   =   100.0;
      Wdike                   =   20000.0;
      T_in                    =   900.0;
  end

  # Create dike
  dike              =   Dike(W=Wdike, H=Hdike,Center=cen,Angle=DikeAngle,Type=DikeType,T=T_in);  # Specify dike

  # Compute velocity required to create space for dike
  Δ                 =   Hdike;          # max. dike opening (m)
  dt                =   1;              # time in which the dike opened fully

  Velocity          =   HostRockVelocityFromDike(Grid, FullGrid, Δ, dt,dike);          # compute velocity field

  if Dimension=="2D"
    Vel      =   Velocity[:];

    if CreatePlots
      Vx,Vz       =   Velocity[1],Velocity[2];
      p1          =   heatmap(x/1e3, z/1e3,      Vx',       aspect_ratio=1, xlims=(x[1]/1e3,x[end]/1e3), ylims=(z[1]/1e3,z[end]/1e3),   c=:inferno, title="2D Vx",  dpi=300, levels=30)
      p2          =   heatmap(x/1e3, z/1e3,      Vz',       aspect_ratio=1, xlims=(x[1]/1e3,x[end]/1e3), ylims=(z[1]/1e3,z[end]/1e3),   c=:inferno, title="Vz",  dpi=300, levels=30)

      #st=100; Xv=X[:]; Zv=Z[:];
      #quiver!(Xv[1:step:end]./1e3, Zv[1:step:end]./1e3, gradient=(Vx[1:step:end],Vz[1:step:end]), arrow = :arrow)

      plot(p1,p2);

      png("HostRockVelocity_$(Dimension)_$(DikeType)")
    end

  elseif Dimension=="3D"
    Vel      =   Velocity[:];

    if CreatePlots
      Vx,Vy,Vz    =   Velocity[1],Velocity[2],Velocity[3];
      p1          =   heatmap(x/1e3, z/1e3,      Vx[:,Int((Ny-1)/2),:]',       aspect_ratio=1, xlims=(x[1]/1e3,x[end]/1e3), ylims=(z[1]/1e3,z[end]/1e3),   c=:inferno, title="2D Vx",  dpi=300, levels=30)
      p2          =   heatmap(x/1e3, z/1e3,      Vz[:,Int((Ny-1)/2),:]',       aspect_ratio=1, xlims=(x[1]/1e3,x[end]/1e3), ylims=(z[1]/1e3,z[end]/1e3),   c=:inferno, title="Vz",  dpi=300, levels=30)

      plot(p1,p2);
      png("HostRockVelocity_$(Dimension)_$(DikeType)")

      # write this to a paraview VTK file, using the package WriteVTK.jl
      #vtkfile = vtk_grid("HostVelocity_3D", Vector(x/1e3), Vector(y/1e3), Vector(z/1e3)) # 3-D
      #vtkfile["Velocity"] = (Vx,Vy,Vz);
      #outfiles = vtk_save(vtkfile)
    end

  end

  return norm(Vel,2);        # return measure of Vel
end

function test_InjectDike(Dimension="2D", DikeType="ElasticDike", DikeAngle=[45], numDikeInjectionEvents=1; InterpolationMethod="Cubic", AdvectionMethod="RK2")
  # tests dike insertion in the domain including adding tracers

  if Dimension=="2D"
    # Model parameters
    W,H                     =   30.0,  30.0;                                # Width, Length, Height

    # Define grid
    Nx, Nz                  =   129, 129;                                     # resolution of coarse grid
    dx,dz                   =   W*1e3/(Nx-1), H*1e3/(Nz-1);                         # grid size [m]
    x,z                     =   0:dx:W*1e3, -H*1e3:dz:0;                            # 1D coordinate arrays
    coords                  =   collect(Iterators.product(x,z))               # generate coordinates from 1D coordinate vectors
    X,Z                     =   (x->x[1]).(coords), (x->x[2]).(coords);       # transfer coords to 3D arrays
    Grid, GridFull,Spacing  =   (x,z), (X,Z), (dx,dz);

    Hdike                   =   1000.0;
    Wdike                   =   20000.0;
    T_in                    =   900.0;

    cen                     =   [W/2;-H/2].*1e3;
  elseif Dimension=="3D"
      # Model parameters
      W,L,H                   =   30., 30., 30.;                                    # Width, Length, Height

      # Define coarse grid
      Nx, Ny, Nz              =   129,129,129;                                                    # resolution of coarse grid
      dx,dy,dz                =   W*1e3/(Nx-1), L*1e3/(Ny-1), H*1e3/(Nz-1);                     # grid size [m]
      x,y,z                   =   0:dx:((Nx-1)*dx),  0:dy:((Ny-1)*dy), -((Nz-1)*dz):dz:0.;      # 1D coordinate arrays
      coords                  =   collect(Iterators.product(x,y,z))                             # generate coordinates from 1D coordinate vectors
      X,Y,Z                   =   (x->x[1]).(coords), (x->x[2]).(coords), (x->x[3]).(coords);   # transfer coords to 3D arrays
      Grid, GridFull,Spacing  =   (x,y,z), (X,Y,Z), (dx,dy,dz);
      cen                     =   [W/2; L/2; -H/2].*1e3;

      Hdike                   =   1000.0;
      Wdike                   =   20000.0;
      T_in                    =   900.0;
  end

  # Create BG temperature structure
  GeoT                    =   20;
  T                       =   -Z./1e3.*GeoT;                                             # initial (linear) temperature profile

  # Create dike
  dike                    =   Dike(W=Wdike, H=Hdike,Center=cen,Angle=DikeAngle,Type=DikeType,T=T_in);  # Specify dike

  # Test the InsertDike routine, which modifies temperature and adds tracers
  nTr_dike                =   1000;
  Tracers                 =   StructArray{Tracer}(undef, 1)                                    # Initialize Tracers structure
  Tracers, Tnew, InjectVol, Velocity    =   InjectDike(Tracers, T, Grid, dike, nTr_dike,  InterpolationMethod=InterpolationMethod, AdvectionMethod=AdvectionMethod);           # Inject first dike

  for i=1:numDikeInjectionEvents-1
    T = Tnew;
    Tracers, Tnew, InjectVol, Velocity  =   InjectDike(Tracers, T, Grid, dike, nTr_dike, InterpolationMethod=InterpolationMethod, AdvectionMethod=AdvectionMethod);           # Inject more dikes
  end

  if Dimension=="2D"

    if CreatePlots
      Vx = Velocity[1];
      Vz = Velocity[2];

      Tr_coord    =   Tracers.coord; Tr_coord = hcat(Tr_coord...)';       # extract array with coordinates of tracers
      p1          =   heatmap(x/1e3, z/1e3,      T',     aspect_ratio=1, xlims=(x[1]/1e3,x[end]/1e3), ylims=(z[1]/1e3,z[end]/1e3),   c=:inferno, title="T",  dpi=300, levels=30)
      p2          =   scatter(Tr_coord[:,1]/1e3, Tr_coord[:,2]/1e3, zcolor = Tracers.T, m = (:inferno , 0.8, Plots.stroke(0.01, :black)), markersize=5.0, xlims=(x[1]/1e3,x[end]/1e3), ylims=(z[1]/1e3,z[end]/1e3),title="Tracers")
      p3          =   heatmap(x/1e3, z/1e3,      Vx',       aspect_ratio=1, xlims=(x[1]/1e3,x[end]/1e3), ylims=(z[1]/1e3,z[end]/1e3),   c=:inferno, title="Vx",  dpi=300, levels=30)
      p4          =   heatmap(x/1e3, z/1e3,      Vz',       aspect_ratio=1, xlims=(x[1]/1e3,x[end]/1e3), ylims=(z[1]/1e3,z[end]/1e3),   c=:inferno, title="Vz",  dpi=300, levels=30)

      plot(p1,p2,p3,p4);

      png("InsertDike_$(Dimension)_$(DikeType)")
    end

  elseif Dimension=="3D"

    if CreatePlots
      Vx = Velocity[1];
      Vz = Velocity[3];

      Tr_coord    =   Tracers.coord; Tr_coord = hcat(Tr_coord...)';       # extract array with coordinates of tracers
      p1          =   heatmap(x/1e3, z/1e3,     T[:,Int(ceil(Ny/2)),:]',       aspect_ratio=1, xlims=(x[1]/1e3,x[end]/1e3), ylims=(z[1]/1e3,z[end]/1e3),   c=:inferno, title="T",  dpi=300, levels=30)
      p2          =   scatter(Tr_coord[:,1]/1e3, Tr_coord[:,3]/1e3, zcolor = Tracers.T, m = (:inferno , 0.8, Plots.stroke(0.01, :black)), markersize=5.0,title="Tracers",xlims=(x[1]/1e3,x[end]/1e3), ylims=(z[1]/1e3,z[end]/1e3),)
      p3          =   heatmap(x/1e3, z/1e3,      Vx[:,Int(ceil(Ny/2)),:]',       aspect_ratio=1, xlims=(x[1]/1e3,x[end]/1e3), ylims=(z[1]/1e3,z[end]/1e3),   c=:inferno, title="Vx",  dpi=300, levels=30)
      p4          =   heatmap(x/1e3, z/1e3,      Vz[:,Int(ceil(Ny/2)),:]',       aspect_ratio=1, xlims=(x[1]/1e3,x[end]/1e3), ylims=(z[1]/1e3,z[end]/1e3),   c=:inferno, title="Vz",  dpi=300, levels=30)

      plot(p1,p2,p3,p4);
      png("InsertDike_$(Dimension)_$(DikeType)")

      # write this to a paraview VTK file, using the package WriteVTK.jl
      #vtkfile = vtk_grid("InsertDike_3D", Vector(x/1e3), Vector(y/1e3), Vector(z/1e3)) # 3-D
      #vtkfile["Temperature"] = (T);
      #vtkfile["Velocity"]    = (Velocity);
      #outfiles = vtk_save(vtkfile)
    end

  end

  return norm(T[:],2);
end

# ===================================================================================================

if 1==1

@testset "Dike_Velocity" begin
  @test test_HostRockVelocityFromDike("2D", "ElasticDike",[80    ])  ≈   2663.677375120158  atol=1e-8;
  @test test_HostRockVelocityFromDike("2D","SquareDike",  [80    ])  ≈   5286.539510870982  atol=1e-8;
  @test test_HostRockVelocityFromDike("3D","SquareDike",  [90; 90])  ≈  13114.877048604001  atol=1e-4;
  @test test_HostRockVelocityFromDike("3D","ElasticDike",[90; 45])   ≈   4762.014274270334  atol=1e-4;
end

# test dike structure
@testset "Dike_Struct" begin
  @test typeof(Dike(Center=[0; 0],Angle=[45],     T=900, Type="SquareDike",  W=1000, H=100 ))==Dike
  @test typeof(Dike(Center=[0; 0],Angle=[45; 90], T=900, Type="ElasticDike", W=1000, H=100 ))==Dike
  @test typeof(Dike(Center=[0; 0],Angle=[45; 90], T=800, Type="ElasticDike", Q=1e6, ΔP=1e7, E=1e10) )==Dike
end

# Volume of dike
@testset "Dike_Volume" begin
  @test volume_dike(Dike(Center=[0; 0],Angle=[45; 90], T=800, Type="ElasticDike", Q=1e6, ΔP=1e7, E=1e10))[1] ≈   2539.6349734808196 atol=1e-8;
  @test volume_dike(Dike(Center=[0; 0],Angle=[45],     T=900, Type="SquareDike",  W=1000, H=100 ))[1] ≈ 100000.0  atol=1e-8;
end

# Dike insertion algorithm
@testset "Dike_Inject" begin
  @test test_InjectDike("2D", "SquareDike", [80 ],1)    ≈   47525.46469221513 atol=1e-3;

  @test test_InjectDike("2D", "ElasticDike",[45 ],2, InterpolationMethod="Linear")    ≈   48448.71557320294  atol=1e-3;     # also tests what happens if we add 2 dikes
  @test test_InjectDike("2D", "ElasticDike",[45 ],2, InterpolationMethod="Quadratic") ≈   48770.72282124797  atol=1e-3;     # also tests what happens if we add 2 dikes
  @test test_InjectDike("2D", "ElasticDike",[45 ],2, InterpolationMethod="Cubic")     ≈   48782.180721943965 atol=1e-3;     # also tests what happens if we add 2 dikes

  @test test_InjectDike("3D", "ElasticDike",[80; 45])   ≈   519654.9176188356 atol=1e-3;
  @test test_InjectDike("3D", "SquareDike",[15; -30])   ≈   527521.5507505678  atol=1e-3;

end

end

#test_InjectDike("2D", "ElasticDike", [80 ],5, InterpolationMethod="Linear", AdvectionMethod="Euler")
#test_InjectDike("3D", "SquareDike",[80; 45], InterpolationMethod="Linear", AdvectionMethod="RK2")


text:
module PoissonSBM

using ..SBMs

export PoisSBM

import Distributions

type PoisSBM <: SBM

  G::SparseMatrixCSC
  g::Vector{Int}
  t::Vector{Vector{Int}}
  dir::Bool
  degcorr::Bool

end

# Just a helper "copy" constructor mainly for situations when the identified type has changed after a package reload
PoisSBM(x) = PoisSBM(x.G, copy(x.g), x.t, copy(x.dir), copy(x.degcorr))

#Each "ps" is an Any[] with the following fields (Dicts were REALLY slow...)
@enum PARAMFIELDS CL_COUNTS=1 NODE_INDEG=2 NODE_OUTDEG=3 CL_INDEG=4 CL_OUTDEG=5 CL_EDGE_COUNTS=6

function make_dummy( N=40, dir=false, degcorr=true )

  G = spzeros(UInt8,N,N)
  #3/4 axon, 1/4 dend; 1/2 g1
  t = [[1,2],[3,4]] #1&2 -> axon, 3&4 -> dend
  g = ones(Int,(N,))

  for i in eachindex(g)
    if     i > round(Int,7*N/8)  g[i] = 4
    elseif i > round(Int,6*N/8)  g[i] = 3
    elseif i > round(Int,3*N/8)  g[i] = 2
    else                         g[i] = 1
    end
  end

  p = [0  0  2  2;
       0  0  1  3;
       0  0  0  0;
       0  0  0  0];

  #NOTE: dummy model currently doesn't include much degree variation
  dists = [Distributions.Poisson(p[i,j]) for i in 1:4, j in 1:4];

  SBMs.fill_G!(G, g, dists, dir)

  PoisSBM(G,g,t,dir,degcorr)
end

function make_full_dummy( N=40, dir=false, degcorr=false )

  G = spzeros(UInt8,N,N)

  t = [[1,2]]
  g = zeros(Int,(N,))

  for i in eachindex(g)
    if i > round(Int,3*N/4)  g[i] = 2
    else                     g[i] = 1
    end
  end

  p = [1  5;
       2  1]

  #NOTE: dummy model currently doesn't include much degree variation
  dists = [Distributions.Poisson(p[i,j]) for i in 1:2, j in 1:2]

  SBMs.fill_G!(G, g, dists, dir)

  PoisSBM(G,g,t,dir,degcorr)
end

function SBMs.computeparams(sbm::PoisSBM, g=nothing)

  if g == nothing  g = sbm.g  end

  cluster_counts = count_groups(g)

  node_indeg, node_outdeg = compute_node_degrees(sbm, g)

  cl_indeg, cl_outdeg, cl_edges = group_degrees(node_indeg, node_outdeg, g)

  Any[cluster_counts, node_indeg, node_outdeg,
   cl_indeg, cl_outdeg, cl_edges]
end

function count_groups(g)

  max_group = maximum(g)

  counts = zeros(Int,max_group)

  for v in g  counts[v] += 1  end

  counts
end

function SBMs.updateparams!(sbm::PoisSBM, ps, old_g, g=nothing)

  if g == nothing  g = sbm.g  end

  update_group_counts!(ps, g)
  update_node_degrees!(sbm, ps, old_g, g)

  cl_indeg, cl_outdeg, cl_edges = group_degrees(ps[Int(NODE_INDEG)],
                                                ps[Int(NODE_OUTDEG)],g)

  ps[Int(CL_INDEG)]  = cl_indeg;
  ps[Int(CL_OUTDEG)] = cl_outdeg;
  ps[Int(CL_EDGE_COUNTS)] = cl_edges;

end

function update_group_counts!(ps,g)

  #recomputing is as fast as updating in this case
  ps[Int(CL_COUNTS)] = count_groups(g)

end

function update_node_degrees!(sbm::PoisSBM, ps, old_g, g)

  @assert typeof(old_g) == typeof(g)
  @assert size(old_g)   == size(g)

  node_indeg = ps[Int(NODE_INDEG)]
  node_outdeg = ps[Int(NODE_OUTDEG)]

  for i in eachindex(g)
    curr_group, old_group = g[i],old_g[i]

    if curr_group == old_group  continue  end

    senders, sender_vs = findnz(sbm.G[:,i])
    receivers, rec_vs  = findnz(sbm.G[i,:])

    for (s,sv) in zip(senders, sender_vs)
      node_outdeg[s,old_group]  -= sv
      node_outdeg[s,curr_group] += sv
    end

    for (r,rv) in zip(receivers,rec_vs)
      node_indeg[r,old_group] -= rv
      node_indeg[r,curr_group] += rv
    end

  end

end

function compute_node_degrees(sbm, g=nothing)

  if g == nothing g = sbm.g end

  G = sbm.G
  n = size(G,1)
  max_group = maximum(g)

  #assumes g is 1:max_group
  node_indeg  = zeros(Int,(n,max_group))
  node_outdeg = zeros(Int,(n,max_group))

  rows = rowvals(G)
  vals = nonzeros(G)
  for i in 1:n, j in nzrange(G,i)
    r = rows[j]; v = vals[j]

    node_indeg[i,g[r]] += v
    node_outdeg[r,g[i]] += v
  end

  node_indeg, node_outdeg
end

function group_degrees(node_indeg, node_outdeg, g)

  max_group = maximum(g)

  cl_edges = zeros(Int,max_group,max_group)

  sx,sy = size(node_outdeg)
  for j in 1:sy, i in 1:sx
    if node_outdeg[i,j] == 0 continue end
    cl_edges[g[i],j] += node_outdeg[i,j]
  end

  sum(cl_edges,1), sum(cl_edges,2), cl_edges
end

function SBMs.considermoves(sbm::PoisSBM, i, g=nothing, ps=nothing; forcemove=true)

  if g == nothing   g = sbm.g                end
  if ps == nothing  ps=computeparams(sbm,g)  end

  gi = g[i]

  gi_type = gi in sbm.t[1] ? sbm.t[1] : sbm.t[2]
  # g_type = collect(filter( x -> gi in x, sbm.t ))
  # @assert length(g_type) == 1 "Group exists within multilple types"

  moves = Int[]
  logliks = Float64[]

  for group in gi_type #g_type[1]

    if gi == group && forcemove  continue  end

    push!(moves,group)

    new_params = tweak_params(ps, i, gi, group)

    push!(logliks, loglikelihood(sbm, nothing, new_params))
  end

  moves, logliks
end

function tweak_params(ps, i, orig_group, new_group)

  cluster_counts = copy(ps[Int(CL_COUNTS)])

  cluster_counts[orig_group] -= 1
  cluster_counts[new_group]  += 1

  nbor_incls = ps[Int(NODE_INDEG)][i,:]
  nbor_outcls = ps[Int(NODE_OUTDEG)][i,:]

  cl_edges = copy(ps[Int(CL_EDGE_COUNTS)])

  for cl in eachindex(nbor_incls)
    cl_edges[orig_group,cl] -= nbor_outcls[cl]
    cl_edges[cl,orig_group] -= nbor_incls[cl]

    cl_edges[new_group,cl]  += nbor_outcls[cl]
    cl_edges[cl,new_group]  += nbor_incls[cl]
  end

  cl_indeg = sum(cl_edges,1)
  cl_outdeg = sum(cl_edges,2)

  Any[ cluster_counts, nothing, nothing,
       cl_indeg, cl_outdeg, cl_edges ]
end

function SBMs.loglikelihood(sbm::PoisSBM, g=nothing, ps=nothing)

  if g == nothing   g = sbm.g                  end
  if ps == nothing  ps = computeparams(sbm,g)  end

  cl_count    = ps[Int(CL_COUNTS)]
  cl_indeg    = ps[Int(CL_INDEG)]
  cl_outdeg   = ps[Int(CL_OUTDEG)]
  edge_counts = ps[Int(CL_EDGE_COUNTS)]

  ll = 0

  pairs = SBMs.assemble_pairs(sbm.t, sbm.dir)
  for (r,s) in pairs

    denom = 0

    if sbm.degcorr
      cl_out_r = cl_outdeg[r]
      cl_in_s  = cl_indeg[s]
      denom = cl_out_r * cl_in_s
    else
      n_r, n_s = cl_count[r], cl_count[s]
      denom = n_r*n_s
    end

    if denom == 0  continue  end

    m_rs = edge_counts[r,s]

    ll += m_rs*log(m_rs/(denom))
  end

  ll
end

#Not ready yet
#function compute_prob_matrix(sbm::PoisSBM, g=nothing, ps=nothing)
#
#  if g == nothing   g = sbm.g                  end
#  if ps == nothing  ps = computeparams(sbm,g)  end
#end

end #module PoissonSBM


text:
@withkw @trait immutable IFParameter
    τm::Float = 20ms
    τe::Float = 5ms
    τi::Float = 10ms
    Vt::Float = -50mV
    Vr::Float = -60mV
    El::Float = Vr
end

@withkw @trait type IF
    param::IFParameter = IFParameter()
    N::Int = 100
    v::Vector{Float} = param.Vr + rand(N) * (param.Vt - param.Vr)
    ge::Vector{Float} = zeros(N)
    gi::Vector{Float} = zeros(N)
    fire::Vector{Bool} = zeros(Bool, N)
    I::Vector{Float} = zeros(N)
    records::Dict = Dict()
end

@replace function integrate!(p::IF, param::IFParameter, dt::Float)
    @inbounds for i = 1:N
        v[i] += dt * (ge[i] + gi[i] - (v[i] - El) + I[i]) / τm
        ge[i] += dt * -ge[i] / τe
        gi[i] += dt * -gi[i] / τi
    end
    @inbounds for i = 1:N
        fire[i] = v[i] > Vt
        v[i] = ifelse(fire[i], Vr, v[i])
    end
end


text:
module TestArrayHelperFunctions

using FactCheck
using DataCubes
using DataCubes: wrap_array, gtake

facts("ArrayHelperFunctions tests") do
  @fact sum(@nalift([1,20,3])) --> x -> DataCubes.naop_eq(Nullable(24), x).value
  @fact sum(@nalift([1.0,20.0,3.0])) --> x -> DataCubes.naop_eq(Nullable(24.0), x).value
  @fact sum(@nalift([1,20,3,NA])) --> x -> DataCubes.naop_eq(Nullable(24), x).value
  @fact sum(@nalift([1.0,20.0,3.0,NA])) --> x -> DataCubes.naop_eq(Nullable(24.0), x).value
  @fact sum(@nalift([1,2,3,NA,50])) --> x -> DataCubes.naop_eq(Nullable(56), x).value
  @fact sum(@nalift([1.0,2.0,3.0,NA,50.0])) --> x -> DataCubes.naop_eq(Nullable(56.0), x).value
  @fact prod(@nalift([1,2,5])) --> x -> DataCubes.naop_eq(x, Nullable(10)).value
  @fact prod(@nalift([1.0,2.0,5.0])) --> x -> DataCubes.naop_eq(x, Nullable(10.0)).value
  @fact prod(@nalift([1,2,3,NA])) --> x -> DataCubes.naop_eq(x, Nullable(6)).value
  @fact prod(@nalift([1.0,2.0,3.0,NA])) --> x -> DataCubes.naop_eq(x, Nullable(6.0)).value
  @fact prod(@nalift([1,2,3,NA,15])) --> x -> DataCubes.naop_eq(x, Nullable(90)).value
  @fact prod(@nalift([1.0,2.0,3.0,NA,15.0])) --> x -> DataCubes.naop_eq(x, Nullable(90.0)).value
  @fact diff(nalift([1 2 3;6 5 4]),1,2) --> nalift([1 -4 -2;5 3 1])
  @fact diff(1.0*nalift([1 2 3;6 5 4]),1,2) --> 1.0*nalift([1 -4 -2;5 3 1])
  @fact collect(DataCubes.dropnaiter(@nalift([1,20,3,NA,50,NA,NA]))) --> [1,20,3,50]
  @fact collect(DataCubes.dropnaiter(@nalift([1.0,20.0,3.0,NA,50.0,NA,NA]))) --> [1.0,20.0,3.0,50.0]
  @fact DataCubes.simplify_array(map(x->x[1],collect(DataCubes.enum_dropnaiter(@nalift([1,20,3,NA,50,NA,NA]))))) --> map(x->x[1],DataCubes.wrap_array(collect(zip([1,2,3,5],@nalift([1,20,3,50])))))
  @fact DataCubes.simplify_array(map(x->x[1],collect(DataCubes.enum_dropnaiter(@nalift([1.0,20.0,3.0,NA,50.0,NA,NA]))))) --> map(x->x[1],DataCubes.wrap_array(collect(zip([1.0,2.0,3.0,5.0],@nalift([1.0,20.0,3.0,50.0])))))
  @fact DataCubes.simplify_array(map(x->x[2],collect(DataCubes.enum_dropnaiter(@nalift([1,20,3,NA,50,NA,NA]))))) --> map(x->x[2],DataCubes.wrap_array(collect(zip([1,2,3,5],[1,20,3,50]))))
  @fact DataCubes.simplify_array(map(x->x[2],collect(DataCubes.enum_dropnaiter(@nalift([1.0,20.0,3.0,NA,50.0,NA,NA]))))) --> map(x->x[2],DataCubes.wrap_array(collect(zip([1.0,2.0,3.0,5.0],[1.0,20.0,3.0,50.0]))))
  @fact collect(DataCubes.zip_dropnaiter(@nalift([1,2,3,NA,NA,4,5,7,8,9]), @nalift([NA,2,3,NA,10,14,NA,NA,10]))) --> [(2,2),(3,3),(4,14),(8,10)]
  @fact collect(DataCubes.zip_dropnaiter(@nalift([1.0,2.0,3.0,NA,NA,4.0,5.0,7.0,8.0,9.0]), @nalift([NA,2,3,NA,10,14,NA,NA,10]))) --> [(2.0,2),(3.0,3),(4.0,14),(8.0,10)]
  @fact cov(@nalift([1,2,35,NA,NA]),@nalift([10,NA,12,NA,5])).value --> roughly(cov([1,35],[10,12]))
  @fact cov(@nalift([1.0,2.0,35.0,NA,NA]),@nalift([10.0,NA,12.0,NA,5.0])).value --> roughly(cov([1.0,35.0],[10.0,12.0]))
  @fact cov(@nalift([1.0,2.0,35.0,NA,NA]),@nalift([10,NA,12,NA,5])).value --> roughly(cov([1,35],[10,12]))
  @fact cor(@nalift([1,2,35,NA,NA]),@nalift([10,NA,12,NA,5])).value --> roughly(cor([1,35],[10,12]))
  @fact cor(@nalift([1.0,2.0,35.0,NA,NA]),@nalift([10.0,NA,12.0,NA,5.0])).value --> roughly(cor([1,35],[10,12]))
  @fact median(@nalift([1,10,100,NA])).value --> roughly(10.0)
  @fact median(@nalift([1.0,10.0,100.0,NA])).value --> roughly(10.0)
  @fact median(@nalift([1,10,NA,100])).value --> roughly(10.0)
  @fact median(@nalift([1.0,10.0,NA,100.0])).value --> roughly(10.0)
  @fact middle(@nalift([1,10,NA,100])).value --> roughly(50.5)
  @fact middle(@nalift([1.0,10.0,NA,100.0])).value --> roughly(50.5)
  @fact middle(@nalift([1,10,NA,100])).value --> roughly(50.5)
  @fact middle(@nalift([1.0,10.0,NA,100.0])).value --> roughly(50.5)
  @fact quantile(@nalift([1,10,NA,100,2,3,NA]),0.0).value --> roughly(1.0)
  @fact quantile(@nalift([1.0,10.0,NA,100.0,2.0,3.0,NA]),0.0).value --> roughly(1.0)
  @fact quantile(@nalift([1,10,NA,100,2,3,NA]),0.25).value --> roughly(2.0)
  @fact quantile(@nalift([1.0,10.0,NA,100.0,2.0,3.0,NA]),0.25).value --> roughly(2.0)
  @fact quantile(@nalift([1,10,NA,100,2,3,NA]),0.5).value --> roughly(3.0)
  @fact quantile(@nalift([1.0,10.0,NA,100.0,2.0,3.0,NA]),0.5).value --> roughly(3.0)
  @fact quantile(@nalift([1,10,NA,100,2,3,NA]),0.75).value --> roughly(10.0)
  @fact quantile(@nalift([1.0,10.0,NA,100.0,2.0,3.0,NA]),0.75).value --> roughly(10.0)
  @fact quantile(@nalift([1,10,NA,100,2,3,NA]),1.0).value --> roughly(100.0)
  @fact quantile(@nalift([1.0,10.0,NA,100.0,2.0,3.0,NA]),1.0).value --> roughly(100.0)
  @fact nalift(quantile(@nalift([1,10,NA,100,2,3,NA]),[0.0,0.25,0.5,0.75,1.0])) --> nalift([1.0,2.0,3.0,10.0,100.0])
  @fact nalift(quantile(2.0*@nalift([1,10,NA,100,2,3,NA]),[0.0,0.25,0.5,0.75,1.0])) --> 2.0*nalift([1.0,2.0,3.0,10.0,100.0])
  @fact quantile(nalift([Nullable{Int}(),Nullable{Int}()]),0.5) --> isnull
  @fact cumsum(larr(a=reshape(1:10,2,5), axis1=darr(k=[:x,:y]))) --> larr(a=[1 3 5 7 9;3 7 11 15 19], axis1=darr(k=[:x,:y]))
  @fact cumsum(larr(a=1.0 * reshape(1:10,2,5), axis1=darr(k=[:x,:y]))) --> larr(a=1.0 * [1 3 5 7 9;3 7 11 15 19], axis1=darr(k=[:x,:y]))
  @fact cumsum(larr(a=reshape(1:10,2,5), axis1=darr(k=[:x,:y])),1) --> larr(a=[1 3 5 7 9;3 7 11 15 19], axis1=darr(k=[:x,:y]))
  @fact cumsum(larr(a=1.0*reshape(1:10,2,5), axis1=darr(k=[:x,:y])),1) --> larr(a=1.0*[1 3 5 7 9;3 7 11 15 19], axis1=darr(k=[:x,:y]))
  @fact cumsum(larr(a=reshape(1:10,2,5), axis1=darr(k=[:x,:y])),2) --> larr(a=[1 4 9 16 25;2 6 12 20 30], axis1=darr(k=[:x,:y]))
  @fact cumsum(larr(a=1.0*reshape(1:10,2,5), axis1=darr(k=[:x,:y])),2) --> larr(a=1.0*[1 4 9 16 25;2 6 12 20 30], axis1=darr(k=[:x,:y]))
  @fact cumprod(larr(a=reshape(1:10,2,5), axis1=darr(k=[:x,:y]))) --> larr(a=[1 3 5 7 9;2 12 30 56 90], axis1=darr(k=[:x,:y]))
  @fact cumprod(larr(a=1.0*reshape(1:10,2,5), axis1=darr(k=[:x,:y]))) --> larr(a=1.0*[1 3 5 7 9;2 12 30 56 90], axis1=darr(k=[:x,:y]))
  @fact cumprod(larr(a=reshape(1:10,2,5), axis1=darr(k=[:x,:y])),1) --> larr(a=[1 3 5 7 9;2 12 30 56 90], axis1=darr(k=[:x,:y]))
  @fact cumprod(larr(a=1.0*reshape(1:10,2,5), axis1=darr(k=[:x,:y])),1) --> larr(a=1.0*[1 3 5 7 9;2 12 30 56 90], axis1=darr(k=[:x,:y]))
  @fact cumprod(larr(a=reshape(1:10,2,5), axis1=darr(k=[:x,:y])),2) --> larr(a=[1 3 15 105 945;2 8 48 384 3840], axis1=darr(k=[:x,:y]))
  @fact cumprod(larr(a=1.0*reshape(1:10,2,5), axis1=darr(k=[:x,:y])),2) --> larr(a=1.0*[1 3 15 105 945;2 8 48 384 3840], axis1=darr(k=[:x,:y]))
  @fact cumsum(@nalift([1 2 3;4 5 NA])) --> nalift([1 2 3;5 7 3])
  @fact cumsum(@nalift([1.0 2.0 3.0;4.0 5.0 NA])) --> nalift([1.0 2.0 3.0;5.0 7.0 3.0])
  @fact cumprod(@nalift([1 2 3;4 5 NA])) --> nalift([1 2 3;4 10 3])
  @fact cumprod(@nalift([1.0 2.0 3.0;4.0 5.0 NA])) --> 1.0*nalift([1 2 3;4 10 3])
  @fact mean(larr(a=[1 2 3;4 5 6],b=[1.0 2.0 3.0;4.0 5.0 6.0], axis1=[:x,:y]))[:a].value --> roughly(3.5)
  @fact mean(larr(a=1.0*[1 2 3;4 5 6],b=1.0*[1.0 2.0 3.0;4.0 5.0 6.0], axis1=[:x,:y]))[:a].value --> roughly(3.5)
  @fact cor(darr(a=[1 2;3 4],b=[5 6;7 8]))[1,2].value --> roughly(1.0)
  @fact cor(darr(a=1.0*[1 2;3 4],b=1.0*[5 6;7 8]))[1,2].value --> roughly(1.0)
  @fact cor(darr(a=[1.0 2.0;3.0 4.0],b=[5 6;7 8]))[1,2].value --> roughly(1.0)
  @fact cor(larr(a=[1 2;3 4],b=[5 6;7 8]))[1,2].value --> roughly(1.0)
  @fact cor(larr(a=1.0*[1 2;3 4],b=1.0*[5 6;7 8]))[1,2].value --> roughly(1.0)
  @fact cor(larr(a=[1.0 2.0;3.0 4.0],b=[5 6;7 8]))[1,2].value --> roughly(1.0)
  @fact cov(larr([1,3,5]*1.0), larr(1.0*[2,4,6])).value --> roughly(4.0)
  @fact cor(larr([1,3,5]*1.0), larr(1.0*[2,4,6])).value --> roughly(1.0)
  @fact sum(@rap _ .* _ cor(larr(a=[1,3,5]*1.0,b=[1,2,3]), larr(a=1.0*[2,4,6],b=[5,4,3])) - larr(1.0*[1 -1;1 -1],axis=[:a,:b],axis=[:a,:b])).value < 0.00000001 --> true
  @fact sum(@rap _ .* _ cov(larr(a=[1,3,5]*1.0,b=[1,2,3]), larr(a=1.0*[2,4,6],b=[5,4,3])) - larr([4 -2;2 -1],axis=[:a,:b],axis=[:a,:b])).value --> roughly(0.0)
  @fact maximum(@nalift([1,10,NA,100,2,3,NA])).value --> roughly(quantile(@nalift([1,10,NA,100,2,3,NA]),1).value)
  @fact maximum(1.0*@nalift([1,10,NA,100,2,3,NA])).value --> roughly(1.0*quantile(@nalift([1,10,NA,100,2,3,NA]),1).value)
  @fact minimum(@nalift([1,10,NA,100,2,3,NA])).value --> roughly(quantile(@nalift([1,10,NA,100,2,3,NA]),0).value)
  @fact minimum(1.0*@nalift([1,10,NA,100,2,3,NA])).value --> roughly(1.0*quantile(@nalift([1,10,NA,100,2,3,NA]),0).value)
  @fact mmaximum(@nalift([1 NA 3;4 NA NA]),window=2,1,2) --> @nalift([1 4 3;4 NA 3])
  @fact mmaximum(@nalift([1.0 NA 3.0;4.0 NA NA]),window=2,1,2) --> @nalift([1.0 4.0 3.0;4.0 NA 3.0])
  @fact mmaximum(@nalift([1 NA 3;4 NA NA]),window=2,1,2) --> mmaximum(@nalift([1 NA 3;4 NA NA]),window=2,2,1)
  @fact_throws mmaximum(@nalift([1 NA 3;4 NA NA]),window=2,1,1)
  @fact mmaximum(@nalift([1 NA 3;4 NA NA]),window=2,1) --> @nalift([1 NA 3;4 NA 3])
  @fact mmaximum(1.0*@nalift([1 NA 3;4 NA NA]),window=2,1) --> 1.0*@nalift([1 NA 3;4 NA 3])
  @fact mmaximum(@nalift([1 NA 3;4 NA NA]),window=2,2) --> @nalift([1 1 3;4 4 NA])
  @fact mmaximum(1.0*@nalift([1 NA 3;4 NA NA]),window=2,2) --> 1.0*@nalift([1 1 3;4 4 NA])
  @fact mminimum(@nalift([1 NA 3;4 NA NA]),window=2,1) --> @nalift([1 NA 3;1 NA 3])
  @fact mminimum(1.0*@nalift([1 NA 3;4 NA NA]),window=2,1) --> 1.0*@nalift([1 NA 3;1 NA 3])
  @fact mminimum(@nalift([1 NA 3;4 NA NA]),2,window=2) --> @nalift([1 1 3;4 4 NA])
  @fact mminimum(1.0*@nalift([1 NA 3;4 NA NA]),2,window=2) --> 1.0*@nalift([1 1 3;4 4 NA])
  @fact mmedian(@nalift([1 NA 3;4 NA NA]),window=2,1) --> @nalift([1.0 NA 3.0;2.5 NA 3.0])
  @fact mmedian(1.0*@nalift([1 NA 3;4 NA NA]),window=2,1) --> @nalift([1.0 NA 3.0;2.5 NA 3.0])
  @fact mmedian(@nalift([1 NA 3;4 NA NA]),2,window=2) --> @nalift([1.0 1.0 3.0;4.0 4.0 NA])
  @fact mmedian(1.0*@nalift([1 NA 3;4 NA NA]),2,window=2) --> @nalift([1.0 1.0 3.0;4.0 4.0 NA])
  @fact mmiddle(@nalift([1 NA 3;4 NA NA]),window=2,1) --> @nalift([1.0 NA 3.0;2.5 NA 3.0])
  @fact mmiddle(1.0*@nalift([1 NA 3;4 NA NA]),window=2,1) --> @nalift([1.0 NA 3.0;2.5 NA 3.0])
  @fact mmiddle(@nalift([1 NA 3;4 NA NA]),window=2,2) --> @nalift([1.0 1.0 3.0;4.0 4.0 NA])
  @fact mmiddle(1.0*@nalift([1 NA 3;4 NA NA]),window=2,2) --> @nalift([1.0 1.0 3.0;4.0 4.0 NA])
  @fact mmiddle(@nalift([1 NA 3;4 NA NA]),1) --> @nalift([1.0 NA 3.0;2.5 NA 3.0])
  @fact mmiddle(1.0*@nalift([1 NA 3;4 NA NA]),1) --> @nalift([1.0 NA 3.0;2.5 NA 3.0])
  @fact mmiddle(@nalift([1 NA 3;4 NA NA]),2) --> @nalift([1.0 1.0 2.0;4.0 4.0 4.0])
  @fact mmiddle(1.0*@nalift([1 NA 3;4 NA NA]),2) --> @nalift([1.0 1.0 2.0;4.0 4.0 4.0])
  @fact mmean(@nalift([1 NA 3;4 NA NA]),window=2,1) --> @nalift([1.0 NA 3.0;2.5 NA 3.0])
  @fact mmean(1.0*@nalift([1 NA 3;4 NA NA]),window=2,1) --> @nalift([1.0 NA 3.0;2.5 NA 3.0])
  @fact mmean(@nalift([1 NA 3;4 NA NA]),2,window=2) --> @nalift([1.0 1.0 3.0;4.0 4.0 NA])
  @fact mmean(1.0*@nalift([1 NA 3;4 NA NA]),2,window=2) --> @nalift([1.0 1.0 3.0;4.0 4.0 NA])
  @fact mquantile(@nalift([1 NA 3;4 NA NA]),0.5,window=2,1) --> mmedian(@nalift([1 NA 3;4 NA NA]),window=2,1)
  @fact mquantile(1.0*@nalift([1 NA 3;4 NA NA]),0.5,window=2,1) --> mmedian(@nalift([1 NA 3;4 NA NA]),window=2,1)
  @fact mquantile(@nalift([1 NA 3;4 NA NA]),0.5,window=2,2) --> mmedian(@nalift([1 NA 3;4 NA NA]),2,window=2)
  @fact mquantile(1.0*@nalift([1 NA 3;4 NA NA]),0.5,window=2,2) --> mmedian(1.0*@nalift([1 NA 3;4 NA NA]),2,window=2)
  @fact mmedian(@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5)[8].value --> roughly(5.0)
  @fact mmedian(1.0*@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5)[8].value --> roughly(5.0)
  @fact mmean(@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5)[8].value --> roughly(17/3)
  @fact mmean(1.0*@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5)[8].value --> roughly(1.0*17/3)
  @fact mminimum(@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5)[8].value --> 4
  @fact mminimum(1.0*@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5)[8].value --> 4.0
  @fact mmaximum(@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5)[8].value --> 8
  @fact mmaximum(1.0*@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5)[8].value --> 8.0
  @fact mmiddle(@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5)[8].value --> roughly(6.0)
  @fact mmiddle(1.0*@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5)[8].value --> roughly(6.0)
  @fact cumsum(@nalift([1,2,NA,8,4,NA,NA,5,8])) --> nalift([1,3,3,11,15,15,15,20,28])
  @fact cumsum(1.0*@nalift([1,2,NA,8,4,NA,NA,5,8])) --> 1.0*nalift([1,3,3,11,15,15,15,20,28])
  @fact cumprod(@nalift([1,2,NA,8,4,NA,NA,5,8])) --> nalift([1,2,2,16,64,64,64,320,2560])
  @fact cumprod(1.0*@nalift([1,2,NA,8,4,NA,NA,5,8])) --> 1.0*nalift([1,2,2,16,64,64,64,320,2560])
  @fact cummin(@nalift([1,2,NA,8,4,NA,NA,5,8])) --> nalift([1,1,1,1,1,1,1,1,1])
  @fact cummin(1.0*@nalift([1,2,NA,8,4,NA,NA,5,8])) --> 1.0*nalift([1,1,1,1,1,1,1,1,1])
  @fact cummax(@nalift([1,2,NA,8,4,NA,NA,5,8])) --> nalift([1,2,2,8,8,8,8,8,8])
  @fact cummax(1.0*@nalift([1,2,NA,8,4,NA,NA,5,8])) --> 1.0*nalift([1,2,2,8,8,8,8,8,8])
  @fact msum(@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5) --> nalift([1,3,3,11,15,14,12,17,17])
  @fact msum(1.0*@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5) --> 1.0*nalift([1,3,3,11,15,14,12,17,17])
  @fact mprod(@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5) --> nalift([1,2,2,16,64,64,32,160,160])
  @fact mprod(1.0*@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5) --> 1.0*nalift([1,2,2,16,64,64,32,160,160])
  @fact mquantile(@nalift([1,2,NA,8,4,NA,NA,5,8]), 0.25, window=5, 1)[8].value --> roughly(4.5)
  @fact mquantile(1.0*@nalift([1,2,NA,8,4,NA,NA,5,8]), 0.25, window=5, 1)[8].value --> roughly(4.5)

  context("nullable array tests") do
    @fact maximum(@nalift([1,10,NA,100,2,3,NA])).value --> roughly(quantile(@nalift([1,10,NA,100,2,3,NA]),1).value)
    @fact maximum(1.0*@nalift([1,10,NA,100,2,3,NA])).value --> roughly(quantile(@nalift([1,10,NA,100,2,3,NA]),1).value)
    @fact minimum(@nalift([1,10,NA,100,2,3,NA])).value --> roughly(quantile(@nalift([1,10,NA,100,2,3,NA]),0).value)
    @fact minimum(1.0*@nalift([1,10,NA,100,2,3,NA])).value --> roughly(quantile(@nalift([1,10,NA,100,2,3,NA]),0).value)
    @fact mmaximum(@nalift([1 NA 3;4 NA NA]),window=2,1) --> @nalift([1 NA 3;4 NA 3])
    @fact mmaximum(1.0*@nalift([1 NA 3;4 NA NA]),window=2,1) --> 1.0*@nalift([1 NA 3;4 NA 3])
    @fact mmaximum(@nalift([1 NA 3;4 NA NA]),2,window=2) --> @nalift([1 1 3;4 4 NA])
    @fact mmaximum(1.0*@nalift([1 NA 3;4 NA NA]),2,window=2) --> 1.0*@nalift([1 1 3;4 4 NA])
    @fact mminimum(@nalift([1 NA 3;4 NA NA]),1,window=2) --> @nalift([1 NA 3;1 NA 3])
    @fact mminimum(1.0*@nalift([1 NA 3;4 NA NA]),1,window=2) --> 1.0*@nalift([1 NA 3;1 NA 3])
    @fact mminimum(@nalift([1 NA 3;4 NA NA]),2,window=2) --> @nalift([1 1 3;4 4 NA])
    @fact mminimum(1.0*@nalift([1 NA 3;4 NA NA]),2,window=2) --> 1.0*@nalift([1 1 3;4 4 NA])
    @fact mmedian(@nalift([1 NA 3;4 NA NA]),window=2,1) --> @nalift([1.0 NA 3.0;2.5 NA 3.0])
    @fact mmedian(1.0*@nalift([1 NA 3;4 NA NA]),window=2,1) --> 1.0*@nalift([1.0 NA 3.0;2.5 NA 3.0])
    @fact mmedian(@nalift([1 NA 3;4 NA NA]),2,window=2) --> @nalift([1.0 1.0 3.0;4.0 4.0 NA])
    @fact mmedian(1.0*@nalift([1 NA 3;4 NA NA]),2,window=2) --> 1.0*@nalift([1.0 1.0 3.0;4.0 4.0 NA])
    @fact mmiddle(@nalift([1 NA 3;4 NA NA]),window=2,1) --> @nalift([1.0 NA 3.0;2.5 NA 3.0])
    @fact mmiddle(1.0*@nalift([1 NA 3;4 NA NA]),window=2,1) --> 1.0*@nalift([1.0 NA 3.0;2.5 NA 3.0])
    @fact mmiddle(@nalift([1 NA 3;4 NA NA]),2,window=2) --> @nalift([1.0 1.0 3.0;4.0 4.0 NA])
    @fact mmiddle(1.0*@nalift([1 NA 3;4 NA NA]),2,window=2) --> 1.0*@nalift([1.0 1.0 3.0;4.0 4.0 NA])
    @fact mmean(@nalift([1 NA 3;4 NA NA]),window=2,1) --> @nalift([1.0 NA 3.0;2.5 NA 3.0])
    @fact mmean(1.0*@nalift([1 NA 3;4 NA NA]),window=2,1) --> @nalift([1.0 NA 3.0;2.5 NA 3.0])
    @fact mmean(@nalift([1 NA 3;4 NA NA]),2,window=2) --> @nalift([1.0 1.0 3.0;4.0 4.0 NA])
    @fact mmean(1.0*@nalift([1 NA 3;4 NA NA]),2,window=2) --> @nalift([1.0 1.0 3.0;4.0 4.0 NA])
    @fact mquantile(@nalift([1 NA 3;4 NA NA]),0.5,window=2) --> mmedian(@nalift([1 NA 3;4 NA NA]),window=2,1)
    @fact mquantile(1.0*@nalift([1 NA 3;4 NA NA]),0.5,window=2) --> mmedian(@nalift([1 NA 3;4 NA NA]),window=2,1)
    @fact mquantile(@nalift([1 NA 3;4 NA NA]),0.5,window=2,2) --> mmedian(@nalift([1 NA 3;4 NA NA]),2,window=2)
    @fact mquantile(1.0*@nalift([1 NA 3;4 NA NA]),0.5,window=2,2) --> mmedian(@nalift([1 NA 3;4 NA NA]),2,window=2)
    @fact mmedian(@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5)[8].value --> roughly(5.0)
    @fact mmedian(1.0*@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5)[8].value --> roughly(5.0)
    @fact mmean(@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5)[8].value --> roughly(17/3)
    @fact mmean(1.0*@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5)[8].value --> roughly(1.0*17/3)
    @fact mminimum(@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5)[8].value --> 4
    @fact mminimum(1.0*@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5)[8].value --> 4.0
    @fact mmaximum(@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5)[8].value --> 8
    @fact mmaximum(1.0*@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5)[8].value --> 8.0
    @fact mmiddle(@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5)[8].value --> roughly(6.0)
    @fact mmiddle(1.0*@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5)[8].value --> roughly(6.0)
    @fact cumsum(@nalift([1,2,NA,8,4,NA,NA,5,8])) --> nalift([1,3,3,11,15,15,15,20,28])
    @fact cumsum(1.0*@nalift([1,2,NA,8,4,NA,NA,5,8])) --> 1.0*nalift([1,3,3,11,15,15,15,20,28])
    @fact cumprod(@nalift([1,2,NA,8,4,NA,NA,5,8])) --> nalift([1,2,2,16,64,64,64,320,2560])
    @fact cumprod(1.0*@nalift([1,2,NA,8,4,NA,NA,5,8])) --> 1.0*nalift([1,2,2,16,64,64,64,320,2560])
    @fact cummin(@nalift([1,2,NA,8,4,NA,NA,5,8])) --> nalift([1,1,1,1,1,1,1,1,1])
    @fact cummin(1.0*@nalift([1,2,NA,8,4,NA,NA,5,8])) --> 1.0*nalift([1,1,1,1,1,1,1,1,1])
    @fact cummax(@nalift([1,2,NA,8,4,NA,NA,5,8])) --> nalift([1,2,2,8,8,8,8,8,8])
    @fact cummax(1.0*@nalift([1,2,NA,8,4,NA,NA,5,8])) --> 1.0*nalift([1,2,2,8,8,8,8,8,8])
    @fact msum(@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5) --> nalift([1,3,3,11,15,14,12,17,17])
    @fact msum(1.0*@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5) --> 1.0*nalift([1,3,3,11,15,14,12,17,17])
    @fact mprod(@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5) --> nalift([1,2,2,16,64,64,32,160,160])
    @fact mprod(1.0*@nalift([1,2,NA,8,4,NA,NA,5,8]), window=5) --> 1.0*nalift([1,2,2,16,64,64,32,160,160])
    @fact mquantile(@nalift([1,2,NA,8,4,NA,NA,5,8]), 0.25, window=5)[8].value --> roughly(4.5)
    @fact mquantile(1.0*@nalift([1,2,NA,8,4,NA,NA,5,8]), 0.25, window=5)[8].value --> roughly(4.5)
    @fact mmean(nalift(1:10)) --> nalift([1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5])
    @fact mmean(1.0*nalift(1:10)) --> nalift([1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5])
  end

  context("darr tests") do
    @fact mmaximum(@darr(a=[1 NA 3;4 NA NA]),window=2,1) --> @darr(a=[1 NA 3;4 NA 3])
    @fact mmaximum(@darr(a=[1 NA 3;4 NA NA]),2,window=2) --> @darr(a=[1 1 3;4 4 NA])
    @fact mminimum(@darr(a=[1 NA 3;4 NA NA]),window=2,1) --> @darr(a=[1 NA 3;1 NA 3])
    @fact mminimum(@darr(a=[1 NA 3;4 NA NA]),window=2,2) --> @darr(a=[1 1 3;4 4 NA])
    @fact mmedian(@darr(a=[1 NA 3;4 NA NA]),window=2,1) --> @darr(a=[1.0 NA 3.0;2.5 NA 3.0])
    @fact mmedian(@darr(a=[1 NA 3;4 NA NA]),2,window=2) --> @darr(a=[1.0 1.0 3.0;4.0 4.0 NA])
    @fact mmiddle(@darr(a=[1 NA 3;4 NA NA]),window=2,1) --> @darr(a=[1.0 NA 3.0;2.5 NA 3.0])
    @fact mmiddle(@darr(a=[1 NA 3;4 NA NA]),2,window=2) --> @darr(a=[1.0 1.0 3.0;4.0 4.0 NA])
    @fact mmean(@darr(a=[1 NA 3;4 NA NA]),window=2,1) --> @darr(a=[1.0 NA 3.0;2.5 NA 3.0])
    @fact mmean(@darr(a=[1 NA 3;4 NA NA]),window=2,2) --> @darr(a=[1.0 1.0 3.0;4.0 4.0 NA])
    @fact mquantile(@darr(a=[1 NA 3;4 NA NA]),0.5,window=2,1) --> mmedian(@darr(a=[1 NA 3;4 NA NA]),window=2,1)
    @fact mquantile(@darr(a=[1 NA 3;4 NA NA]),0.5,2,window=2) --> mmedian(@darr(a=[1 NA 3;4 NA NA]),window=2,2)
    @fact mmedian(@darr(a=[1,2,NA,8,4,NA,NA,5,8]), window=5)[8][:a].value --> roughly(5.0)
    @fact mmean(@darr(a=[1,2,NA,8,4,NA,NA,5,8]), window=5)[8][:a].value --> roughly(17/3)
    @fact mminimum(@darr(a=[1,2,NA,8,4,NA,NA,5,8]), window=5)[8][:a].value --> 4
    @fact mmaximum(@darr(a=[1,2,NA,8,4,NA,NA,5,8]), window=5)[:a][8].value --> 8
    @fact mmiddle(@darr(a=[1,2,NA,8,4,NA,NA,5,8]), window=5)[:a][8].value --> roughly(6.0)
    @fact cumsum(@darr(a=[1,2,NA,8,4,NA,NA,5,8])) --> darr(a=[1,3,3,11,15,15,15,20,28])
    @fact cumprod(@darr(a=[1,2,NA,8,4,NA,NA,5,8])) --> darr(a=[1,2,2,16,64,64,64,320,2560])
    @fact cummin(@darr(a=[1,2,NA,8,4,NA,NA,5,8])) --> darr(a=[1,1,1,1,1,1,1,1,1])
    @fact cummax(@darr(a=[1,2,NA,8,4,NA,NA,5,8])) --> darr(a=[1,2,2,8,8,8,8,8,8])
    @fact msum(@darr(a=[1,2,NA,8,4,NA,NA,5,8]), window=5) --> darr(a=[1,3,3,11,15,14,12,17,17])
    @fact mprod(@darr(a=[1,2,NA,8,4,NA,NA,5,8]), window=5) --> darr(a=[1,2,2,16,64,64,32,160,160])
    @fact mquantile(@darr(a=[1,2,NA,8,4,NA,NA,5,8]), 0.25, window=5, 1)[:a][8].value --> roughly(4.5)
    @fact maximum(@darr(a=[1,10,NA,100,2,3,NA]))[:a].value --> roughly(quantile(@darr(a=[1,10,NA,100,2,3,NA]),1)[:a].value)
    @fact minimum(@darr(a=[1,10,NA,100,2,3,NA]))[:a].value --> roughly(quantile(@darr(a=[1,10,NA,100,2,3,NA]),0)[:a].value)
    @fact mmaximum(@darr(a=[1 NA 3;4 NA NA]),window=2,1) --> @darr(a=[1 NA 3;4 NA 3])
    @fact mmaximum(@darr(a=[1 NA 3;4 NA NA]),window=2,2) --> @darr(a=[1 1 3;4 4 NA])
    @fact mminimum(@darr(a=[1 NA 3;4 NA NA]),window=2,1) --> @darr(a=[1 NA 3;1 NA 3])
    @fact mminimum(@darr(a=[1 NA 3;4 NA NA]),2,window=2) --> @darr(a=[1 1 3;4 4 NA])
    @fact mmedian(@darr(a=[1 NA 3;4 NA NA]),window=2,1) --> @darr(a=[1.0 NA 3.0;2.5 NA 3.0])
    @fact mmedian(@darr(a=[1 NA 3;4 NA NA]),2,window=2) --> @darr(a=[1.0 1.0 3.0;4.0 4.0 NA])
    @fact mmiddle(@darr(a=[1 NA 3;4 NA NA]),window=2,1) --> @darr(a=[1.0 NA 3.0;2.5 NA 3.0])
    @fact mmiddle(@darr(a=[1 NA 3;4 NA NA]),2,window=2) --> @darr(a=[1.0 1.0 3.0;4.0 4.0 NA])
    @fact mmean(@darr(a=[1 NA 3;4 NA NA]),window=2,1) --> @darr(a=[1.0 NA 3.0;2.5 NA 3.0])
    @fact mmean(@darr(a=[1 NA 3;4 NA NA]),2,window=2) --> @darr(a=[1.0 1.0 3.0;4.0 4.0 NA])
    @fact mquantile(@darr(a=[1 NA 3;4 NA NA]),0.5,window=2) --> mmedian(@darr(a=[1 NA 3;4 NA NA]),window=2,1)
    @fact mquantile(@darr(a=[1 NA 3;4 NA NA]),0.5,2,window=2) --> mmedian(@darr(a=[1 NA 3;4 NA NA]),window=2,2)
    @fact mmedian(@darr(a=[1,2,NA,8,4,NA,NA,5,8]), window=5)[8][:a].value --> roughly(5.0)
    @fact mmean(@darr(a=[1,2,NA,8,4,NA,NA,5,8]), window=5)[:a][8].value --> roughly(17/3)
    @fact mminimum(@darr(a=[1,2,NA,8,4,NA,NA,5,8]), window=5)[:a][8].value --> 4
    @fact mmaximum(@darr(a=[1,2,NA,8,4,NA,NA,5,8]), window=5)[8][:a].value --> 8
    @fact mmiddle(@darr(a=[1,2,NA,8,4,NA,NA,5,8]), window=5)[8][:a].value --> roughly(6.0)
    @fact cumsum(@darr(a=[1,2,NA,8,4,NA,NA,5,8])) --> darr(a=[1,3,3,11,15,15,15,20,28])
    @fact cumprod(@darr(a=[1,2,NA,8,4,NA,NA,5,8])) --> darr(a=[1,2,2,16,64,64,64,320,2560])
    @fact cummin(@darr(a=[1,2,NA,8,4,NA,NA,5,8])) --> darr(a=[1,1,1,1,1,1,1,1,1])
    @fact cummax(@darr(a=[1,2,NA,8,4,NA,NA,5,8])) --> darr(a=[1,2,2,8,8,8,8,8,8])
    @fact msum(@darr(a=[1,2,NA,8,4,NA,NA,5,8]), window=5) --> darr(a=[1,3,3,11,15,14,12,17,17])
    @fact mprod(@darr(a=[1,2,NA,8,4,NA,NA,5,8]), window=5) --> darr(a=[1,2,2,16,64,64,32,160,160])
    @fact mquantile(@darr(a=[1,2,NA,8,4,NA,NA,5,8]), 0.25, window=5)[:a][8].value --> roughly(4.5)
  end

  context("larr tests") do
    @fact maximum(@larr([1,10,NA,100,2,3,NA])).value --> roughly(quantile(@larr([1,10,NA,100,2,3,NA]),1).value)
    @fact minimum(@larr([1,10,NA,100,2,3,NA])).value --> roughly(quantile(@larr([1,10,NA,100,2,3,NA]),0).value)
    @fact mmaximum(@larr([1 NA 3;4 NA NA]),window=2,1) --> @larr([1 NA 3;4 NA 3])
    @fact mmaximum(@larr([1 NA 3;4 NA NA]),2,window=2) --> @larr([1 1 3;4 4 NA])
    @fact mminimum(@larr([1 NA 3;4 NA NA]),window=2,1) --> @larr([1 NA 3;1 NA 3])
    @fact mminimum(@larr([1 NA 3;4 NA NA]),2,window=2) --> @larr([1 1 3;4 4 NA])
    @fact mmedian(@larr([1 NA 3;4 NA NA]),window=2,1) --> @larr([1.0 NA 3.0;2.5 NA 3.0])
    @fact mmedian(@larr([1 NA 3;4 NA NA]),1) --> @larr([1.0 NA 3.0;2.5 NA 3.0])
    @fact mmedian(@larr([1 NA 3;4 NA NA]),2,window=2) --> @larr([1.0 1.0 3.0;4.0 4.0 NA])
    @fact mmedian(@larr([1 NA 3;4 NA NA]),2) --> @larr([1.0 1.0 2.0;4.0 4.0 4.0])
    @fact mmiddle(@larr([1 NA 3;4 NA NA]),window=2,1) --> @larr([1.0 NA 3.0;2.5 NA 3.0])
    @fact mmiddle(@larr([1 NA 3;4 NA NA]),window=2,2) --> @larr([1.0 1.0 3.0;4.0 4.0 NA])
    @fact mmiddle(@larr([1 NA 3 5 2 8;4 NA NA 7 2 NA]),2) --> larr([1.0 1.0 2.0 3.0 3.0 4.5;4.0 4.0 4.0 5.5 4.5 4.5])
    @fact mmiddle(@larr([1 NA 3 5 2 8;4 NA NA 7 2 NA]),2,rev=true) --> @larr([4.5 5.0 5.0 5.0 5.0 8.0;4.5 4.5 4.5 4.5 2.0 NA])
    @fact mmean(@larr([1 NA 3;4 NA NA]),window=2,1) --> @larr([1.0 NA 3.0;2.5 NA 3.0])
    @fact mmean(@larr([1 NA 3;4 NA NA]),2) --> @larr([1.0 1.0 2.0;4.0 4.0 4.0])
    @fact mquantile(@larr([1 NA 3;4 NA NA]),0.5,window=2,1) --> mmedian(@larr([1 NA 3;4 NA NA]),window=2,1)
    @fact mquantile(@larr([1 NA 3;4 NA NA]),0.5,window=2,2) --> mmedian(@larr([1 NA 3;4 NA NA]),2,window=2)
    @fact mmedian(@larr([1,2,NA,8,4,NA,NA,5,8]), window=5)[8].value --> roughly(5.0)
    @fact mmean(@larr([1,2,NA,8,4,NA,NA,5,8]), window=5)[8].value --> roughly(17/3)
    @fact mminimum(@larr([1,2,NA,8,4,NA,NA,5,8]), window=5)[8].value --> 4
    @fact mmaximum(@larr([1,2,NA,8,4,NA,NA,5,8]), window=5)[8].value --> 8
    @fact mmiddle(@larr([1,2,NA,8,4,NA,NA,5,8]), window=5)[8].value --> roughly(6.0)
    @fact cumsum(@larr([1,2,NA,8,4,NA,NA,5,8])) --> larr([1,3,3,11,15,15,15,20,28])
    @fact msum(@larr([1,2,NA,8,4,NA,NA,5,8])) --> larr([1,3,3,11,15,15,15,20,28])
    @fact cumprod(@larr([1,2,NA,8,4,NA,NA,5,8])) --> larr([1,2,2,16,64,64,64,320,2560])
    @fact mprod(@larr([1,2,NA,8,4,NA,NA,5,8])) --> larr([1,2,2,16,64,64,64,320,2560])
    @fact cummin(@larr([1,2,NA,8,4,NA,NA,5,8])) --> larr([1,1,1,1,1,1,1,1,1])
    @fact mminimum(@larr([1,2,NA,8,4,NA,NA,5,8])) --> larr([1,1,1,1,1,1,1,1,1])
    @fact cummax(@larr([1,2,NA,8,4,NA,NA,5,8])) --> larr([1,2,2,8,8,8,8,8,8])
    @fact mmaximum(@larr([1,2,NA,8,4,NA,NA,5,8])) --> larr([1,2,2,8,8,8,8,8,8])
    @fact msum(@larr([1,2,NA,8,4,NA,NA,5,8]), window=5) --> larr([1,3,3,11,15,14,12,17,17])
    @fact mprod(@larr([1,2,NA,8,4,NA,NA,5,8]), window=5) --> larr([1,2,2,16,64,64,32,160,160])
    @fact mquantile(@larr([1,2,NA,8,4,NA,NA,5,8]), 0.25, window=5, 1)[8].value --> roughly(4.5)
    @fact mquantile(@larr([1,2,NA,8,4,NA,NA,5,8]), 0.25, window=5)[8].value --> roughly(4.5)
    @fact mquantile(@larr([1,2,NA,8,4,NA,NA,5,8]), 0.25)[8].value --> roughly(2.0)
    @fact mquantile(@larr([1,2,NA,8,4,NA,NA,5,8]), 0.5)[8].value --> roughly(4.0)
    @fact mmedian(@larr([1,2,NA,8,4,NA,NA,5,8]))[8].value --> roughly(4.0)
  end

  context("reverse tests") do
    @fact mprod(@larr([1,2,NA,8,4,NA,NA,5,8]), window=5, rev=true) --> larr([64,64,32,160,160,40,40,40,8])
    @fact mmedian(@larr([1 NA 3;4 NA NA]),window=2,1, rev=true) --> @larr([2.5 NA 3.0;4.0 NA NA])
  end
  context("nafill tests") do
    @fact nafill(@larr([1 NA 3;4 NA NA]),window=1,2) --> @larr([1 NA 3;4 NA NA])
    @fact nafill(@larr([1 NA 3;4 NA NA]),window=2,2) --> @larr([1 1 3;4 4 NA])
    @fact nafill(@larr([1 NA 3;4 NA NA]),window=3,2) --> @larr([1 1 3;4 4 4])
    @fact nafill(@larr([1 NA 3;4 NA NA]),2) --> @larr([1 1 3;4 4 4])
    @fact nafill(@larr([1 NA 3;4 NA NA]),window=1,1) --> @larr([1 NA 3;4 NA NA])
    @fact nafill(@larr([1 NA 3;4 NA NA]),window=2,1) --> @larr([1 NA 3;4 NA 3])
    @fact nafill(@larr([1 NA 3;4 NA NA]),1) --> @larr([1 NA 3;4 NA 3])

    @fact nafill(1.0*@larr([1 NA 3;4 NA NA]),window=1,2) --> 1.0*@larr([1 NA 3;4 NA NA])
    @fact nafill(1.0*@larr([1 NA 3;4 NA NA]),window=2,2) --> 1.0*@larr([1 1 3;4 4 NA])
    @fact nafill(1.0*@larr([1 NA 3;4 NA NA]),window=3,2) --> 1.0*@larr([1 1 3;4 4 4])
    @fact nafill(1.0*@larr([1 NA 3;4 NA NA]),2) --> 1.0*@larr([1 1 3;4 4 4])
    @fact nafill(1.0*@larr([1 NA 3;4 NA NA]),window=1,1) --> 1.0*@larr([1 NA 3;4 NA NA])
    @fact nafill(1.0*@larr([1 NA 3;4 NA NA]),window=2,1) --> 1.0*@larr([1 NA 3;4 NA 3])
    @fact nafill(1.0*@larr([1 NA 3;4 NA NA]),1) --> 1.0*@larr([1 NA 3;4 NA 3])
  end
  #TODO describe has issues because I cannot put Nullable{Int} into Nullable{Float64} anymore in julia v0.6.
  #context("describe tests") do
  #  @fact wrap_array(extract(describe(@darr(a=[1,2,3,4,5])), Nullable(:a)).values) --> nalift([1.0,2.0,3.0,4.0,5.0,3.0,std([1,2,3,4,5]),5.0,0.0,0.0])
  #  @fact keys(extract(describe(@darr(a=[1,2,3,4,5])), Nullable(:a))) --> [:min,:q1,:med,:q3,:max,:mean,:std,:count,:nacount,:naratio]
  #  @fact describe(@darr(a=[1,2,3,NA,5]))[:nacount] --> nalift([Nullable(1)])
  #  @fact describe(@larr(a=[1,2,3,4,5])) --> describe(@darr(a=[1,2,3,4,5]))
  #  @fact describe(nalift([1,2,3,4,5])) --> describe(@larr(a=[1,2,3,4,5]))[1]
  #  @fact describe(@larr(a=[1,NA,3,4,NA], b=[1.0,2.0,3.0,4.0,5.0], axis1['x','y','z','u','v'])) --> describe(@darr(a=[1,NA,3,4,NA], b=[1.0,2.0,3.0,4.0,5.0]))
  #  @fact shift(@nalift([1,2,NA,4,5]), 2) --> @nalift([NA,4,5,NA,NA])
  #  @fact shift(1.0*@nalift([1,2,NA,4,5]), 2) --> 1.0*@nalift([NA,4,5,NA,NA])
  #  @fact shift(@nalift([1,2,NA,4,5]), -2) --> @nalift([NA,NA,1,2,NA])
  #  @fact shift(1.0*@nalift([1,2,NA,4,5]), -2) --> 1.0*@nalift([NA,NA,1,2,NA])
  #  @fact shift(nalift(reshape(1:20,4,5)),1,-2) --> extract(nalift(reshape(1:20,4,5)),2:5,-1:3)
  #  @fact shift(1.0*nalift(reshape(1:20,4,5)),1,-2) --> 1.0*extract(nalift(reshape(1:20,4,5)),2:5,-1:3)
  #  @fact shift(darr(a=[1 2 3;4 5 6]),1,1,isbound=true) --> darr(a=[5 6 6;5 6 6])
  #  @fact shift(darr(a=1.0*[1 2 3;4 5 6]),1,1,isbound=true) --> darr(a=1.0*[5 6 6;5 6 6])
  #  @fact shift(larr(a=[1 2 3;4 5 6],axis2=[:m,:n,:p]),1,1,isbound=true) --> larr(a=[5 6 6;5 6 6],axis2=[:m,:n,:p])
  #  @fact shift(larr(a=1.0*[1 2 3;4 5 6],axis2=[:m,:n,:p]),1,1,isbound=true) --> larr(a=1.0*[5 6 6;5 6 6],axis2=[:m,:n,:p])
  #end
  @fact DataCubes.mquantile_quickselect!(Float64,[30,10,40,20,50],0) --> 10.0
  @fact DataCubes.mquantile_quickselect!(Float64,[30,10,40,20,50],0.25) --> 20.0
  @fact DataCubes.mquantile_quickselect!(Float64,[30,10,40,20,50],0.5) --> 30.0
  @fact DataCubes.mquantile_quickselect!(Float64,[30,10,40,20,50],3//4) --> 40.0
  @fact DataCubes.mquantile_quickselect!(Float64,[30,10,40,20,50],1) --> 50.0
  @fact var(nalift([1,2,3,4,5])).value --> roughly(DataCubes.naop_mul(std(nalift([1,2,3,4,5])),std(nalift([1,2,3,4,5]))).value)
end

end


text:
import Base.getindex, Base.setindex!, Base.*

###################
### Bounds Type ###
###################

type Bounds{T}
    min_x::Float64
    max_x::Float64
    min_y::Float64
    max_y::Float64
end

# fill in the ellipse for the point type if its not provided
call(::Type{Bounds{LLA}}, min_x, max_x, min_y, max_y) = Bounds{add_param(LLA)}(min_x, max_x, min_y, max_y)
call(::Type{Bounds{LL}}, min_x, max_x, min_y, max_y) =  Bounds{add_param(LL)}(min_x, max_x, min_y, max_y)
call(::Type{Bounds{ECEF}}, min_x, max_x, min_y, max_y) = Bounds{add_param(ECEF)}(min_x, max_x, min_y, max_y)
call(::Type{Bounds{ENU}}, min_x, max_x, min_y, max_y) = Bounds{add_param(ENU)}(min_x, max_x, min_y, max_y)

# allow integer indexing
# TODO: learn the Julia convention for this, there must be a not horrible way...
Base.getindex{T}(bbox::Bounds{T}, idx::Int) = getfield(bbox, idx)
Base.setindex!{T}(bbox::Bounds{T}, val::Real, idx::Int) = setfield!(bbox, idx, Float64(val))  # why does Julia have reversed order for these?

# for degs / rads conversions
*{T}(bbox::Bounds{T}, val::Real) = Bounds{T}(bbox.min_x * val, bbox.max_x * val, bbox.min_y * val, bbox.max_y * val)

#############################
### Calculate LL bounds   ###
#############################

function call{T <: Union{LL, LLA}}(::Type{Bounds{T}}, X::Vector{T})

    # build both paths around the circle to the first two points
    bounds_hyps = init_hyps(X, true)

    # and convert both hypotheses
    pick_hyp(bounds_hyps, true)

end

################################
### Calculate other bounds   ###
################################

function call{T}(::Type{Bounds{T}}, X::Vector{T})

    max_x = max_y = -Inf
    min_x = min_y = Inf
    for Xp in X
        min_x, max_x = min(getX(Xp), min_x), max(getX(Xp), max_x)
        min_y, max_y = min(getY(Xp), min_y), max(getY(Xp), max_y)
    end
    return Bounds{T}(min_x, max_x, min_y, max_y)
end

# unspecified type
call{T}(::Type{Bounds}, X::Vector{T}) = Bounds{T}(X)

################################
### Transform Bounds objects ###
################################

# there's not an unambiguous conversion, but for now,
# returning the minimum bounds that contain all points contained
# by the input bounds

# renamed it because it constructs a Bounds type not a ENU Type.  Does it belong as a transform?
function call{T <: ENU, U <: Union{LL, LLA}, V <: Union{LL, LLA}}(::Type{Bounds{T}}, bounds::Bounds{U}, ll_ref::V = center(bounds))

    oT = add_param(add_param(T, U), V)

    max_x = max_y = -Inf
    min_x = min_y = Inf

    xs = [bounds.min_x, bounds.max_x]
    ys = [bounds.min_y, bounds.max_y]

    # AWBS - what's this about?
    if bounds.min_y < 0.0 < bounds.max_y
        push!(ys, 0.0)
    end

    ref_x = getX(ll_ref)

    if bounds.min_x < ref_x < bounds.max_x ||
       (bounds.min_x > bounds.max_x && !(bounds.min_x >= ref_x >= bounds.max_x))
        push!(xs, ref_x)
    end

    for x_ll in xs, y_ll in ys
        pt = geotransform(oT, U(x_ll, y_ll), ll_ref)
        x, y = getX(pt), getY(pt)

        min_x, max_x = min(x, min_x), max(x, max_x)
        min_y, max_y = min(y, min_y), max(y, max_y)
    end

    return Bounds{oT}(min_x, max_x, min_y, max_y)
end

#############################
### Convenience Functions ###
#############################

### Get center point of a Bounds region ###
function center{T}(bounds::Bounds{T})
    x_mid = (bounds.min_x + bounds.max_x) / 2
    y_mid = (bounds.min_y + bounds.max_y) / 2
    return T(x_mid, y_mid)
end

### Check whether a location is within bounds ###
function inBounds{T}(loc::T, bounds::Bounds{T})
    x, y = getX(loc), getY(loc)

    bounds.min_x <= x <= bounds.max_x &&
    bounds.min_y <= y <= bounds.max_y
end

# only for points that have passed the inBounds test
function onBounds{T}(loc::T, bounds::Bounds{T})
    x, y = getX(loc), getY(loc)

    x == bounds.min_x || x == bounds.max_x ||
    y == bounds.min_y || y == bounds.max_y
end

# only for points where inBounds(p1) != inBounds(p2)
# TODO: return actual altitude rather than zero
function boundaryPoint{T}(p1::T, p2::T, bounds::Bounds{T})
    x1, y1 = getX(p1), getY(p1)
    x2, y2 = getX(p2), getY(p2)

    x, y = Inf, Inf

    # Move x to x bound if segment crosses boundary
    if x1 < bounds.min_x < x2 || x1 > bounds.min_x > x2
        x = bounds.min_x
        y = y1 + (y2 - y1) * (bounds.min_x - x1) / (x2 - x1)
    elseif x1 < bounds.max_x < x2 || x1 > bounds.max_x > x2
        x = bounds.max_x
        y = y1 + (y2 - y1) * (bounds.max_x - x1) / (x2 - x1)
    end

    p3 = T(x, y)
    inBounds(p3, bounds) && return p3

    # Move y to y bound if segment crosses boundary
    if y1 < bounds.min_y < y2 || y1 > bounds.min_y > y2
        x = x1 + (x2 - x1) * (bounds.min_y - y1) / (y2 - y1)
        y = bounds.min_y
    elseif y1 < bounds.max_y < y2 || y1 > bounds.max_y > y2
        x = x1 + (x2 - x1) * (bounds.max_y - y1) / (y2 - y1)
        y = bounds.max_y
    end

    p3 = T(x, y)
    inBounds(p3, bounds) && return p3

    error("Failed to find boundary point.")
end

#############################
### Helpers               ###
#############################

# bound theta in rads, so on output -pi <= theta < pi
bound_theta(theta::Real) = theta - floor((theta+pi) / (2*pi)) * 2*pi

# bound theta in degs, so on output -pi <= theta < pi
bound_thetad(theta::Real) = theta - floor((theta+180) / (360)) * 360

# function to start bounding boxes for the vector of LLA points X
function init_hyps{T <: Union{LLA, LL}}(X::Vector{T}, degs::Bool=true)

    bbox_hyps = [Bounds{T}(NaN, NaN, NaN, NaN),
                 Bounds{T}(NaN, NaN, NaN, NaN)]  # two hypotheses
    sc = degs ? pi / 180 : 1.0

    # check for points
    if (length(X) > 0)
        # unsmart method - build both paths around the circle to the first two points
        cnt_min = length(X)+1
        for dim = 1:2

            # start point
            ts = bound_theta(sc * (X[1][dim]))
            cnt = min(2,length(X))

            # find a second unique point
            while (cnt <= length(X)) && (ts == bound_theta(sc * X[cnt][dim]))
                cnt += 1
            end
            te = (cnt <= length(X)) ? bound_theta(sc * X[cnt][dim]) : ts + eps(ts)

            # build the two paths possible paths between points
            bbox_hyps[1][(dim-1)*2 + 1] = ts
            bbox_hyps[1][(dim-1)*2 + 2] = te > ts ? te : te + 2 * pi

            bbox_hyps[2][(dim-1)*2 + 1] = te
            bbox_hyps[2][(dim-1)*2 + 2] = te > ts ? ts + 2 * pi : ts

            cnt_min = min(cnt_min, cnt)  # where to start
        end

        if (degs)
            bbox_hyps[1] *= 180.0 / pi
            bbox_hyps[2] *= 180.0 / pi
        end

        # now include the rest of the points
        updatebounds!(bbox_hyps, X, cnt_min+1, degs)

    end

    return bbox_hyps # make a matrix out of it

end

# function to update a bounding box for an extra point
# updates bbox_h1 and bbox_h1
function updatebounds!{T <: Union{LLA, LL}}(bbox_hyps::Vector{Bounds{T}}, X::Vector{T}, first::Int=1, degs::Bool=true)

    sc = degs ? pi / 180 : 1.0
    bbox_hyps[1] *= sc
    bbox_hyps[2] *= sc

    # now add everything else top the competing paths
    for i = first:length(X)

        # x
        btheta = bound_theta(sc * getX(X[i]))
        updatebounds_worker!(bbox_hyps[1], btheta, 1)
        updatebounds_worker!(bbox_hyps[2], btheta, 1)

        # y
        btheta = bound_theta(sc * getY(X[i]))
        updatebounds_worker!(bbox_hyps[1], btheta, 2)
        updatebounds_worker!(bbox_hyps[2], btheta, 2)
    end

    # put it back in degs
    bbox_hyps[1] *= (1 / sc)
    bbox_hyps[2] *= (1 / sc)

end

# worker function for updating bounds
# input theta should be -pi <= theta < pi
function updatebounds_worker!{T <: Union{LLA, LL}}(bbox::Bounds{T}, theta::Real, dim::Int)

    # indexing shortcuts
    idx1 = (dim-1)*2+1
    idx2 = (dim-1)*2+2

    # clockwise distance from the lower bound
    dlower = mod(theta - bbox[idx1], 2*pi)  # clockwise distance from lower bound

    # anticlockwise distance between bounds
    dbound = bbox[idx2] - bbox[idx1]

    if (dlower > dbound)

        # clockwise distance from the upper bound
        dupper = mod(theta - bbox[idx2], 2*pi)

        # convert distance from the lower bound to anticlockwise
        dlower = 2*pi - dlower

        # pick the minimum adjustement to the bounds
        if (dlower < dupper)
            bbox[idx1] -= dlower
        else
            bbox[idx2] += dupper
        end
    end
end

# function to turn bounding boxes for each hypothesis into a Bounding box type
function pick_hyp{T <: Union{LLA, LL}}(bbox_hyps::Vector{Bounds{T}}, degs::Bool=true)

    wrap = degs ? 360.0 : 2*pi

    # x
    idx = 1 + (bbox_hyps[1].max_x - bbox_hyps[1].min_x >=  bbox_hyps[2].max_x - bbox_hyps[2].min_x)
    d1b = [bbox_hyps[idx].min_x, bbox_hyps[idx].max_x]
    d1b -= floor((mean(d1b) + wrap / 2) / (wrap)) * wrap  # keep it pretty

    # y
    idx = 1 + (bbox_hyps[1].max_y - bbox_hyps[1].min_y >=  bbox_hyps[2].max_y - bbox_hyps[2].min_y)
    d2b = [bbox_hyps[idx].min_y, bbox_hyps[idx].max_y]
    d2b -= floor((mean(d2b) + wrap / 2) / (wrap)) * wrap  # keep it pretty

    # and return the best combo
    return Bounds{T}(d1b[1], d1b[2], d2b[1], d2b[2])

end


text:
const initialized_modules = Set{String}(["Base", "Main"])

const imported_packages_PkgCompiler = Set{String}()
const dependency_tree = Dict{String, Set{String}}()
const what_to_import = Dict{String, Any}()

init_packages = false

function add_to_dependency_tree(parent::String, child::String)
    if !haskey(dependency_tree, parent)
        dependency_tree[parent] = Set{String}()
    end
    push!(dependency_tree[parent], child)
end

function _init_module_pkgcompiler(into::Module, mod::String)
    uuidkey = identify_package(into, mod)
    if uuidkey === nothing
        return
    end
    modul::Module = require(uuidkey)
    parent = "$(modul)"
    if !in(parent, initialized_modules)
        if parent != "PkgCompiler"
            push!(initialized_modules, parent)
            if haskey(dependency_tree, mod)
                for mod2 in dependency_tree[mod]
                    # @info("Try to init $(mod2)")
                    if mod2 !== nothing
                        _init_module_pkgcompiler(into, mod2)
                    end
                    if !in(mod2, initialized_modules) && haskey(what_to_import, mod2)
                        # @info("PkgCompiler inits submodule \"$(mod2)\"")
                        ccall(:jl_init_restored_modules, Cvoid, (Any,), what_to_import[mod2])
                    end
                end
            end
            if haskey(what_to_import, parent)
                # @info("PkgCompiler inits \"$(modul)\"")
                ccall(:jl_init_restored_modules, Cvoid, (Any,), what_to_import[parent])
            end
        end
    end
    nothing
end

function require(into::Module, mod::Symbol)
    uuidkey = identify_package(into, String(mod))
    # Core.println("require($(PkgId(into)), $mod) -> $uuidkey")
    if uuidkey === nothing
        where = PkgId(into)
        if where.uuid === nothing
            throw(ArgumentError("""
                Package $mod not found in current path:
                - Run `import Pkg; Pkg.add($(repr(String(mod))))` to install the $mod package.
                """))
        else
            s = """
            Package $(where.name) does not have $mod in its dependencies:
            - If you have $(where.name) checked out for development and have
              added $mod as a dependency but haven't updated your primary
              environment's manifest file, try `Pkg.resolve()`.
            - Otherwise you may need to report an issue with $(where.name)"""

            uuidkey = identify_package(PkgId(string(into)), String(mod))
            uuidkey === nothing && throw(ArgumentError(s))

            # fall back to toplevel loading with a warning
            if !(where in modules_warned_for)
                @warn string(
                    full_warning_showed[] ? "" : s, "\n",
                    string("Loading $(mod) into $(where.name) from project dependency, ",
                           "future warnings for $(where.name) are suppressed.")
                ) _module = nothing _file = nothing _group = nothing
                push!(modules_warned_for, where)
            end
            full_warning_showed[] = true
        end
    end
    if _track_dependencies[]
        push!(_require_dependencies, (into, binpack(uuidkey), 0.0))
    end
    _init_module_pkgcompiler(into, String(mod))
    return require(uuidkey)
end

# these return either the array of modules loaded from the path / content given
# or an Exception that describes why it couldn't be loaded
# and it reconnects the Base.Docs.META
function _include_from_serialized(path::String, depmods::Vector{Any})
    sv = ccall(:jl_restore_incremental, Any, (Cstring, Any), path, depmods)
    if isa(sv, Exception)
        return sv
    end
    restored = sv[1]
    if !isa(restored, Exception)
        for M in restored::Vector{Any}
            M = M::Module
            if isdefined(M, Base.Docs.META)
                push!(Base.Docs.modules, M)
            end
            if parentmodule(M) === M
                register_root_module(M)
            end
        end
    end
    if isassigned(sv, 2)
        if init_packages
            ccall(:jl_init_restored_modules, Cvoid, (Any,), sv[2])
        else
            modul = restored[1]
            if length(sv[2]) > 0
                # @info "To future init $(modul)"
                parent = "$(parentmodule(modul))"
                what_to_import[parent] = sv[2]
            end
        end
    end
    return restored
end


text:
module RocketScanOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: scan()" begin

    println("Testing: operator scan()")

    run_proxyshowcheck("Scan", scan(Any, +, 0))
    run_proxyshowcheck("ScanNoSeed", scan(+))

    run_testset([
        (
            source      = from(1:5) |> scan(Int, +, 1),
            values      = @ts([ 2, 4, 7, 11, 16, c ]),
            source_type = Int
        ),
        (
            source      = from(1:5) |> scan(+),
            values      = @ts([ 1, 3, 6, 10, 15, c ]),
            source_type = Int
        ),
        (
            source      = from(1:5) |> scan(Vector{Int}, (d, c) -> [ c..., d ], Int[]),
            values      = @ts([ [1], [1, 2], [1, 2, 3], [1, 2, 3, 4], [1, 2, 3, 4, 5], c ]),
            source_type = Vector{Int}
        ),
        (
            source      = completed(Int) |> scan(+),
            values      = @ts(c),
            source_type = Int
        ),
        (
            source      = completed(Int) |> scan(Int, +, 2),
            values      = @ts(c),
            source_type = Int
        ),
        (
            source      = faulted("e") |> scan(+),
            values      = @ts(e("e")),
        ),
        (
            source = never() |> scan(+),
            values = @ts()
        )
    ])

end

end


text:
using Polynomials

"""
    y = evalRPoly(roots, x, k=1)

Compute the value of a polynomial which is given in terms of its `roots`.
"""
function evalRPoly(roots, x, k=1)
    f = poly(roots) * k
    return f(x)
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
@kwdef mutable struct SDMHandler <: AbstractHandler
    id = 0
end

@register(SDMHandler, SDM, [Number] => [DBit])

function (handler::SDMHandler)(netlist::Netlist,
                               inputs::Vector{Variable},
                               outputs::Vector{Variable})
    outstring = """
        $stdcomment
        // BEGIN sdm$(handler.id)
        sdm #(
                .BIT_WIDTH(BIT_WIDTH),
                .INT_WIDTH(INT_WIDTH)
            ) sdm$(handler.id) (
                .CLK(CLK),
                .nRST(nRST),
                .x($(inputs[1].name)),
                .y($(outputs[1].name))
            );
        // END sdm$(handler.id)
        \n"""

    handler.id += 1

    return outstring
end


text:
#=
Conversion functions between any two of the color systems.

All inputs and outputs are three floats in the range [0.0...1.0]
(with the exception of I and Q, which covers a slightly larger range
 also with the exception of R, G and B, which range from 0 to 255
 also with the exception of H, which range from 0 to 360).

Inputs outside the valid range may cause exceptions or invalid outputs.

Supported color systems:
RGB: Red, Green, Blue components
YIQ: Luminance, Chrominance (used by composite video signals)
HLS: Hue, Luminance, Saturation
HSV: Hue, Saturation, Value
=#
# macro
const one_third = 1/3
const one_sixth = 1/6
const two_third = 2/3
# Util Functions
function _v(m1, m2, hue)
    hue = hue - floor(hue)
    if hue < one_sixth
        return m1+(m2-m1)*hue*6
    end
    if hue < .5
        return m2
    end
    if hue < two_third
        return m1+(m2-m1)*(two_third-hue)*6
    end
    return m1
end

#=
YIQ: used by composite video signals (linear combinations of RGB)
Y: perceived grey level (0.0 == black, 1.0 == white)
I, Q: color components
=#
function rgb2yiq(r, g, b)
    r1,g1,b1 = [r,g,b]/255
    return .3*r1+.59*g1+.11*b1, .6*r1-.28*g1-.32*b1, .21*r1-.52*g1+.31*b1
end

function yiq2rgb(y, i, q)
    r,g,b = y+.948262*i+.624013*q, y-.276066*i-.639810*q, y-1.105450*i+1.729860*q
    r = r<0 ? 0 : (r>1 ? 1 : r)
    g = g<0 ? 0 : (g>1 ? 1 : g)
    b = b<0 ? 0 : (b>1 ? 1 : b)
    return trunc(Int,r*255), trunc(Int,g*255), trunc(Int,b*255)
end

#=
HLS: Hue, Luminance, Saturation
H: position in the spectrum(0~360)
L: color lightness
S: color saturation
=#
function rgb2hls(r, g, b)
    r1,g1,b1 = [r,g,b]/255
    maxc,minc = max([r1,g1,b1]...), min([r1,g1,b1]...)
    l = (maxc+minc)/2
    if minc == maxc
        return .0, l, .0
    end
    s = l<=.5 ? (maxc-minc)/(maxc+minc) : (maxc-minc)/(2-maxc-minc)
    rc,gc,bc = [(maxc-r1),(maxc-g1),(maxc-b1)]/(maxc-minc)
    h = r1==maxc ? (bc-gc) : (g1==maxc ? (2+rc-bc) : (4+gc-rc))
    h0 = h/6
    h = h0 - floor(h0)
    return h*360, l, s
end

function hls2rgb(h, l, s)
    if s == 0
        return trunc(Int,l), trunc(Int,l), trunc(Int,l)
    end
    m2 = l<=.5 ? (l*(1+s)) : (l+s-l*s)
    m1 = 2*l - m2
    h /= 360
    return trunc(Int,_v(m1,m2,h+one_third)*255), trunc(Int,_v(m1,m2,h)*255), trunc(Int,_v(m1,m2,h-one_third)*255)
end

#=
HSV: Hue, Saturation, Value
H: position in the spectrum(0~360)
S: color saturation ("purity")
V: color brightness
=#
function rgb2hsv(r, g, b)
    r1,g1,b1 = [r,g,b]/255
    maxc = max([r1,g1,b1]...)
    minc = min([r1,g1,b1]...)
    if maxc == minc
        return .0, .0, maxc
    end
    s = (maxc-minc)/maxc
    rc,gc,bc = [(maxc-r1),(maxc-g1),(maxc-b1)]/(maxc-minc)
    h = r1==maxc ? (bc-gc) : (g1==maxc ? (2+rc-bc) : (4+gc-rc))
    h0 = h / 6
    h = h0 - floor(h0)
    return h*360, s, maxc
end

function hsv2rgb(h, s, v)
    h /= 360
    if s == 0
        return trunc(Int,v*255), trunc(Int,v*255), trunc(Int,v*255)
    end
    i = trunc(Int,h*6)
    f = h*6 - i
    p,q,t = [(1-s),(1-s*f),(1-s*(1-f))]*v
    i %= 6
    if i == 0
        return trunc(Int,v*255), trunc(Int,t*255), trunc(Int,p*255)
    elseif i == 1
        return trunc(Int,q*255), trunc(Int,v*255), trunc(Int,p*255)
    elseif i == 2
        return trunc(Int,p*255), trunc(Int,v*255), trunc(Int,t*255)
    elseif i == 3
        return trunc(Int,p*255), trunc(Int,q*255), trunc(Int,v*255)
    elseif i == 4
        return trunc(Int,t*255), trunc(Int,p*255), trunc(Int,v*255)
    else
        return trunc(Int,v*255), trunc(Int,p*255), trunc(Int,q*255)
    end
end

function yiq2hls(y, i, q)
    r,g,b = yiq2rgb(y,i,q)
    return rgb2hls(r,g,b)
end

function yiq2hsv(y, i, q)
    r,g,b = yiq2rgb(y,i,q)
    return rgb2hsv(r,g,b)
end

function hls2yiq(h, l, s)
    r,g,b = hls2rgb(h,l,s)
    return rgb2yiq(r,g,b)
end

function hls2hsv(h, l, s)
    r,g,b = hls2rgb(h,l,s)
    return rgb2hsv(r,g,b)
end

function hsv2yiq(h,s,v)
    r,g,b = hsv2rgb(h,s,v)
    return rgb2yiq(r,g,b)
end

function hsv2hls(h,s,v)
    r,g,b = hsv2rgb(h,s,v)
    return rgb2hls(r,g,b)
end


text:
#! /usr/bin/julia

# Rosetta Code, Digital root/Multiplicative digital root

const bs = 10
const excnt = 5

function digitalmultroot{S<:Integer,T<:Integer}(n::S, bs::T=10)
    -1 < n && 1 < bs || throw(DomainError())
    ds = n
    pers = 0
    while bs <= ds
        ds = prod(digits(ds, bs))
        pers += 1
    end
    return (pers, ds)
end

println("Testing Multiplicative Digital Root.\n")
for i in [123321, 7739, 893, 899998]
    (pers, ds) = digitalmultroot(i, bs)
    print(@sprintf("%8d", i))
    print(" has persistence ", pers)
    println(" and digital root ", ds)
end

dmr = zeros(Int, bs, excnt)
hasroom = trues(bs)
dex = ones(Int, bs)

i = 0
while any(hasroom)
    (pers, ds) = digitalmultroot(i, bs)
    ds += 1
    if hasroom[ds]
        dmr[ds, dex[ds]] = i
        dex[ds] += 1
        if dex[ds] > excnt
            hasroom[ds] = false
        end
    end
    i += 1
end

println("\n MDR:    First ", excnt, " numbers having this MDR")
for (i, d) in enumerate(0:(bs-1))
    print(@sprintf("%4d: ", d))
    println(join([@sprintf("%6d", dmr[i, j]) for j in 1:excnt], ","))
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
module ShockTube

using Reexport
@reexport using Unitful
@reexport using PyThermo

export shockjump!, shockjump, driverpressure!, driverpressure, shockcalc!, shockcalc

γ(f::PyThermo.Chemical)  = isentropic_exponent(f)

# Shock jump conditions
function shockjump!(gas, Mach)
    α1 = (γ(gas) + 1) / (γ(gas) - 1)
    PR = (Mach^2 * (1 + α1) - 1) / α1
    TR = PR * (PR + α1) / (1 + α1 * PR)
    u2 = (soundspeed(gas) * (α1 - 1) * (PR - 1) /
          #-----------------------------
           √((1 + α1) * (1 + α1 * PR)))
    gas.P = gas.P * PR
    gas.T = gas.T * TR
    return gas, uconvert(u"m/s", u2)
end

shockjump(gas, Mach) = shockjump!(copy(gas), Mach)

function driverpressure!(driver, driven, Ms)
    γ1, γ4 = γ(driven), γ(driver)
    a1, a4 = soundspeed(driven), soundspeed(driver)
    driver.P = driven.P * (1 + 2γ1 / (γ1 + 1) * (Ms^2 - 1)) *
                    (1 + a1 / a4 * (γ4 - 1) / (γ1 + 1) * (1/Ms - Ms)) ^ (2γ4 / (1 - γ4))
    return driver
end
driverpressure(driver, driven, Ms) = driverpressure!(copy(driver), driven, Ms)

function shockreflect!(shocked, Ms)
    γ2 = γ(shocked)
    lhs = Ms/(Ms^2 - 1) * sqrt(1 + 2(γ2 - 1)/(γ2 + 1)^2 * (Ms^2 - 1) * (γ2 + 1/Ms^2))
    Mr = (1 + sqrt(1 +4*lhs^2))/(2lhs)
    shockjump!(shocked, Mr)
end
shockreflect(shocked, Ms) = shockreflect!(copy(shocked), Ms)

function shockcalc!(driver, driven, Ms)
    shocked, u2 = shockjump(driven, Ms)
    driverpressure!(driver, driven, Ms)
    reflected, _ = shockreflect(shocked, Ms)
    return (driver = driver, driven = driven, shocked=shocked, reflected=reflected, u2=u2)
end
shockcalc(driver, driven, Ms) = shockcalc!(copy(driver), driven, Ms)

end # module


text:
"""
```
init_subspec!(m::Li2020)
```
initializes alternative subspecs of the Li2020 model
"""
function init_subspec!(m::Li2020)
    if subspec(m) == "ss0"
        return
    else
        error("Subspec $(subspec(m)) has not been defined.")
    end
end


text:
using CSV
using DataFrames
using DelimitedFiles
using Glob
using Plots; pgfplotsx()

# Create plots
function plot_gmres_jacobi()

	# Plot labels
	labels = ["Parallel Jacobi", "Serial Jacobi"]

	# Benchmark folders
	folders = ["1e6_1000-samples", "1e7_100-samples"]

	# Plot 1
	# Left: Solver runtimes ilu0 and parilu0
	# Right: Speedup parilu0 over ilu0

	# Iterate over matrix size number of samples folders
	for (ifolder, folder) in enumerate(folders)

		# Get folders with csv files
		csvfolder = joinpath(folder, readdir(folder)[1], "csv")

		# Create plots
		plt1 = plot(ylabel="Runtime in seconds", xlabel="Number of threads",
					legend=:topright, framestyle=:box)
		plt2 = plot(ylabel="Parallel speedup", xlabel="Number of threads",
					legend=:none, framestyle=:box)

		# Runtimes csv files
		csvfile_ser = joinpath(csvfolder,
			"runtimes-solver_par_csr_dgks_jacobi.csv")
		csvfile_par = joinpath(csvfolder,
			"runtimes-solver_par_csr_dgks_parjacobi.csv")

		# Add runtimes from all files to the first plot
		df_par = DataFrame!(CSV.File(csvfile_par, header=true))
		df_ser = DataFrame!(CSV.File(csvfile_ser, header=true))
		y_par = convert(Vector, df_par[!, 5])
		y_ser = convert(Vector, df_ser[!, 5])
		x = convert(Vector, df_ser[!, 1])
		plot!(plt1, x, y_par, label=labels[1])
		plot!(plt1, x, y_ser, label=labels[2])

		# Compute speedups for solver and add to second plot
		y = y_ser ./ y_par
		plot!(plt2, x, y, label=labels[1])

		# Plot runtimes and speedups side by side
		plot(plt1, plt2, layout=(1,2), size=(530,530/2))
		savefig(string("gmres_jacobi_", folder, "_solver.png"))
		savefig(string("gmres_jacobi_", folder, "_solver.tex"))

	end # for

	# Plot 2
	# Left: Precon runtimes ilu0 and parilu0
	# Right: Speedup parilu0 over ilu0

	# Iterate over matrix size number of samples folders
	for (ifolder, folder) in enumerate(folders)

		# Get folders with csv files
		csvfolder = joinpath(folder, readdir(folder)[1], "csv")

		# Create plots
		plt1 = plot(ylabel="Runtime in seconds", xlabel="Number of threads",
					legend=:none, framestyle=:box)
		plt2 = plot(ylabel="Parallel speedup", xlabel="Number of threads",
					legend=:bottomright, framestyle=:box)

		# Runtimes csv files
		csvfile_ser = joinpath(csvfolder,
			"runtimes-precon_par_csr_dgks_jacobi.csv")
		csvfile_par = joinpath(csvfolder,
			"runtimes-precon_par_csr_dgks_parjacobi.csv")

		# Add runtimes from all files to the first plot
		df_par = DataFrame!(CSV.File(csvfile_par, header=true))
		df_ser = DataFrame!(CSV.File(csvfile_ser, header=true))
		y_par = convert(Vector, df_par[!, 5])
		y_ser = convert(Vector, df_ser[!, 5])
		x = convert(Vector, df_ser[!, 1])
		plot!(plt1, x, y_par, label=labels[1])
		plot!(plt1, x, y_ser, label=labels[2])

		# Compute speedups for solver and add to second plot
		y = y_ser ./ y_par
		plot!(plt2, x, y, label=labels[1])
		plot!(plt2, [1], [0], label=labels[2])

		# Plot runtimes and speedups side by side
		plot(plt1, plt2, layout=(1,2), size=(530,530/2))
		savefig(string("gmres_jacobi_", folder, "_precon.png"))
		savefig(string("gmres_jacobi_", folder, "_precon.tex"))

	end # for

end # function

max_stars_repo_name: jkrch/JuliaParGMRES-Plots

max_issues_repo_name: jkrch/JuliaParGMRES-Plots

max_forks_repo_name: jkrch/JuliaParGMRES-Plots


text:
using Gtk

win = GtkWindow("Julius", 400, 200)
"""
GtkWindowLeaf(name="", parent, width-request=-1, height-request=-1,
 visible=TRUE, sensitive=TRUE, app-paintable=FALSE, can-focus=FALSE,
 has-focus=FALSE, is-focus=FALSE, can-default=FALSE, has-default=FALSE,
 receives-default=FALSE, composite-child=FALSE, style, events=0,
 no-show-all=FALSE, has-tooltip=FALSE, tooltip-markup=NULL,
 tooltip-text=NULL, window, double-buffered=TRUE, halign=GTK_ALIGN_FILL,
 valign=GTK_ALIGN_FILL, margin-left=0, margin-right=0, margin-top=0,
margin-bottom=0, margin=0, hexpand=FALSE, vexpand=FALSE, hexpand-set=FALSE,
 vexpand-set=FALSE, expand=FALSE, border-width=0, resize-mode=GTK_RESIZE_QUEUE,
 child, type=GTK_WINDOW_TOPLEVEL, title="My window", role=NULL, resizable=TRUE,
 modal=FALSE, window-position=GTK_WIN_POS_NONE, default-width=-1,
 default-height=-1, destroy-with-parent=FALSE,
hide-titlebar-when-maximized=FALSE, icon, icon-name=NULL,
 screen, type-hint=GDK_WINDOW_TYPE_HINT_NORMAL, skip-taskbar-hint=FALSE,
 skip-pager-hint=FALSE, urgency-hint=FALSE, accept-focus=TRUE,
focus-on-map=TRUE, decorated=TRUE, deletable=TRUE,
gravity=GDK_GRAVITY_NORTH_WEST, transient-for, attached-to,
 opacity=1.000000, has-resize-grip=TRUE, resize-grip-visible=TRUE,
 application, ubuntu-no-proxy=FALSE, is-active=FALSE,
has-toplevel-focus=FALSE, startup-id, mnemonics-visible=TRUE,
 focus-visible=TRUE, )
"""

#b = GtkButton("Talk to me")
#push!(win,b)
#function on_button_clicked(w)
#  println("The button has been clicked")
#end
#signal_connect(on_button_clicked, b, "clicked")

vbox = GtkBox(:v)  # :h makes a horizontal layout, :v a vertical layout
push!(win, vbox)
cancel = GtkButton("Cancel")
ok = GtkButton("OK")
push!(vbox, cancel)
push!(vbox, ok)
showall(win)

visible(win)

visible(win, false)

visible(win)

visible(win, true)

max_stars_repo_name: ashwani-rathee/Chatbots

max_issues_repo_name: ashwani-rathee/Chatbots

max_forks_repo_name: ashwani-rathee/Chatbots


text:
latexbracket(args...; kwargs...) = process_latexify(args...; kwargs..., env=:bracket)

function _latexbracket(x; kwargs...)
    latexstr = LaTeXString( "\\[\n" * latexraw(x; kwargs...) * "\\]\n")
    COPY_TO_CLIPBOARD && clipboard(latexstr)
    return latexstr
end

_latexbracket(args...; kwargs...) = latexbracket(args; kwargs...)


text:
#!/usr/bin/env julia
#
# Jerzy Wroczyński (nr. indeksu 250075)
# Obliczenia Naukowe, Laboratorium
# Lista 3, Zadanie 6. metoda bisekcji, funkcja `f_2`
#

include("iterative-methods.jl")
using .IterativeMethods

fl = IterativeMethods.FL

a = fl(-0.5)
b = fl(0.5)

if length(ARGS) > 0
    global a = fl(-0.29)
    global b = fl(1)
    println("warning: using different interval")
end

r, fr, k, err =  IterativeMethods.mbisekcji(
    x -> x * exp(-x),
    a,
    b,
    fl(10^(-5)),
    fl(10^(-5))
)

println(
    "r = ", r, "\n",
    "f(r) = ", fr, "\n",
    "k = ", k, "\n",
    err == 0 ? "" : string("error code: ", err)
)

max_stars_repo_name: jerry-sky/academic-notebook

max_issues_repo_name: jerry-sky/academic-notebook

max_forks_repo_name: jerry-sky/academic-notebook


text:
module ConvNet

include("base.jl")
include("layer.jl")
include("convolution_layer.jl")
include("dropout_layer.jl")
include("dropconnect_layer.jl")
include("fully_connected_layer.jl")
include("input_layer.jl")
include("loss_layers.jl")
include("maxout_layer.jl")
include("net.jl")
include("nonlinear_layers.jl")
include("sgd.jl")
include("trainer.jl")

end # module


text:
using Subsystems
using Test
# using TimerOutputs

# macro timed_testset(name::String, block)
#     # copied from https://github.com/KristofferC/Tensors.jl/blob/master/test/runtests.jl#L8
#     return quote
#         @timeit "$($(esc(name)))" begin
#             @testset "$($(esc(name)))" begin
#                 $(esc(block))
#             end
#         end
#     end
# end

@testset "Subsystems.jl" begin
    # include("test_UIO.jl")
    include("test_ldti.jl")
    include("test_UMVE.jl")
    include("test_LSS.jl")
end


text:
module HierarchicalUtils

using DataStructures
using OrderedCollections: OrderedDict

abstract type NodeType end
struct LeafNode <: NodeType end
struct InnerNode <: NodeType end

const PairVec = Vector{<:Pair}

NodeType(::Type{T}) where T = error("Define NodeType(::Type{$T}) to be either LeafNode() or InnerNode()")
NodeType(x::T) where T = NodeType(T)

isleaf(n) = isleaf(NodeType(n), n)
isleaf(::LeafNode, _) = true
isleaf(::InnerNode, n) = nchildren(n) == 0

children(n) = children(NodeType(n), n)
children(::LeafNode, _) = ()
children(_, ::T) where T =
    error("Define children(n::$T) to return a collection of children (one of allowed types)")

printchildren(n) = children(n)

# TODO incorporate setfield
# set_children(n::T, chs::U) where {T, U} = error("Define set_children(n::$T, chs::$U) where chs are new children to use PreOrder maps")

nodeshow(io, x) = Base.show(io, x)
nodecommshow(io, x) = print(io)

nchildren(n) = nchildren(NodeType(n), n)
nchildren(::LeafNode, n) = 0
nchildren(::InnerNode, n) = length(children(n))
nprintchildren(n) = nchildren(n)

export NodeType, LeafNode, InnerNode
export children, nchildren, isleaf
export nodeshow, nodecommshow, printchildren

include("utilities.jl")

include("statistics.jl")
export nnodes, nleafs, treeheight

include("traversal_encoding.jl")
export encode_traversal, walk, list_traversal, find_traversal

include("printing.jl")
export printtree

include("iterators.jl")
export NodeIterator, LeafIterator, TypeIterator, PredicateIterator, MultiIterator
export PreOrder, PostOrder, LevelOrder
export traverse!

include("maps.jl")
export treemap!, leafmap!, typemap!, predicatemap!, treemap

include("definitions.jl")
export @hierarchical_dict, @hierarchical_vector, @primitives

end # module