text:
module Test_linsol

using DescriptorSystems
using LinearAlgebra
using Polynomials
using Test

println("Test_linsol")
@testset "linsol " begin

@testset "grsol" begin

m = 1; mf = 0
sysgf = dss(ones(1,1));
@time sysx, info, sysgen = grsol(sysgf, mf, solgen = true)
isysg = 1:m; isysf = m+1:m+mf;
mnull = size(sysgen,2)-mf; i0 = 1:mf; inull = mf+1:mf+mnull
@test gnrank(sysgf[:,isysg]*sysx-sysgf[:,isysf]) == 0 &&
      gnrank(sysgf[:,isysg]*sysgen[:,i0]-sysgf[:,isysf], atol = 1.e-7) == 0 &&
      gnrank(sysgf[:,isysg]*sysgen[:,inull], atol = 1.e-7) == 0 &&
      gnrank(sysgf[:,isysg]*(sysgen[:,i0]+sysgen[:,inull]*rdss(3,mnull,mf))-sysgf[:,isysf]) == 0 &&
      info.nrank == 1 && info.rdeg == Int[]

m = 1; mf = 0
sysg = dss(ones(m,m)); sysf = dss(ones(m,mf));
@time sysx, info, sysgen = grsol(sysg, sysf, solgen = true)
isysg = 1:m; isysf = m+1:m+mf;
mnull = size(sysgen,2)-mf; i0 = 1:mf; inull = mf+1:mf+mnull
@test gnrank(sysg*sysx-sysf) == 0 &&
      gnrank(sysg*sysgen[:,i0]-sysf, atol = 1.e-7) == 0 &&
      gnrank(sysg*sysgen[:,inull], atol = 1.e-7) == 0 &&
      gnrank(sysg*(sysgen[:,i0]+sysgen[:,inull]*rdss(3,mnull,mf))-sysf) == 0 &&
      info.nrank == 1 && info.rdeg == Int[]

m = 5; mf = 2
sysgf = dss(rand(m,m+mf));
@time sysx, info, sysgen = grsol(sysgf, mf, solgen = true)
isysg = 1:m; isysf = m+1:m+mf;
mnull = size(sysgen,2)-mf; i0 = 1:mf; inull = mf+1:mf+mnull
@test gnrank(sysgf[:,isysg]*sysx-sysgf[:,isysf], atol=1.e-7) == 0 &&
      gnrank(sysgf[:,isysg]*(sysgen[:,i0]+sysgen[:,inull]*rdss(3,mnull,mf))-sysgf[:,isysf], atol = 1.e-7) == 0 &&
      info.nrank == 5 && info.rdeg == [0, 0] && info.nr == 0

m = 5; mf = 2
sysg = dss(rand(m,m)); sysf = dss(rand(m,mf));
@time sysx, info, sysgen = grsol(sysg, sysf, solgen = true)
isysg = 1:m; isysf = m+1:m+mf;
mnull = size(sysgen,2)-mf; i0 = 1:mf; inull = mf+1:mf+mnull
@test gnrank(sysg*sysx-sysf, atol=1.e-7) == 0 &&
      gnrank(sysg*(sysgen[:,i0]+sysgen[:,inull]*rdss(3,mnull,mf))-sysf, atol = 1.e-7) == 0 &&
      info.nrank == 5 && info.rdeg == [0, 0] && info.nr == 0

p = 5; m = 8; mf = 2
sysg = dss(rand(p,m)); sysf = dss(rand(p,mf));
@time sysx, info, sysgen = grsol(sysg, sysf, solgen = true)
isysg = 1:m; isysf = m+1:m+mf;
mnull = size(sysgen,2)-mf; i0 = 1:mf; inull = mf+1:mf+mnull
@test gnrank(sysg*sysx-sysf, atol=1.e-7) == 0 &&
      gnrank(sysg*(sysgen[:,i0]+sysgen[:,inull]*rdss(3,mnull,mf))-sysf, atol = 1.e-7) == 0 &&
      info.nrank == 5 && info.rdeg == [0, 0]  && info.nr == 0

# Wang and Davison Example (1973)
s = Polynomial([0, 1],'s');
gn = [ s+1 s+3 s^2+3*s; s+2 s^2+2*s 0 ];
gd = [s^2+3*s+2 s^2+3*s+2 s^2+3*s+2 ;
s^2+3*s+2 s^2+3*s+2 s^2+3*s+2];
m = 3; mf = 2;
sysg = dss(gn,gd,minimal = true, atol = 1.e-7);
sysf = dss([1 0;0 1.]);

@time sysx, info, sysgen = grsol(sysg, sysf, solgen = true)
isysg = 1:m; isysf = m+1:m+mf;
mnull = m - info.nrank; i0 = 1:mf; inull = mf+1:mf+mnull
@test gnrank(sysg*sysx-sysf, atol=1.e-7) == 0 &&
      gnrank(sysg*sysgen[:,i0]-sysf, atol = 1.e-7) == 0 &&
      gnrank(sysg*sysgen[:,inull], atol = 1.e-7) == 0 &&
      gnrank(sysg*(sysgen[:,i0]+sysgen[:,inull]*rdss(3,mnull,mf))-sysf, atol = 1.e-7) == 0 &&
      info.nrank == 2 && info.rdeg == [0, 0] && info.nr == 3

evals=[-1, -2, -3];
@time sysx, info, sysgen = grsol(sysg, sysf, poles = evals, solgen = true)
isysg = 1:m; isysf = m+1:m+mf;
mnull = m - info.nrank; i0 = 1:mf; inull = mf+1:mf+mnull
@test gnrank(sysg*sysx-sysf, atol=1.e-7) == 0 &&
      gnrank(sysg*sysgen[:,i0]-sysf, atol = 1.e-7) == 0 &&
      gnrank(sysg*sysgen[:,inull], atol = 1.e-7) == 0 &&
      gnrank(sysg*(sysgen[:,i0]+sysgen[:,inull]*rdss(3,mnull,mf))-sysf, atol = 1.e-7) == 0 &&
      info.nrank == 2 && info.rdeg == [0, 0] && info.nr == 3 &&
      sort(gpole(sysx)) ≈ sort(evals)

sdeg = -1;
@time sysx, info, sysgen = grsol(sysg, sysf, sdeg = sdeg, solgen = true)
isysg = 1:m; isysf = m+1:m+mf;
mnull = m - info.nrank; i0 = 1:mf; inull = mf+1:mf+mnull
@test gnrank(sysg*sysx-sysf, atol=1.e-7) == 0 &&
      gnrank(sysg*sysgen[:,i0]-sysf, atol = 1.e-7) == 0 &&
      gnrank(sysg*sysgen[:,inull], atol = 1.e-7) == 0 &&
      gnrank(sysg*(sysgen[:,i0]+sysgen[:,inull]*rdss(3,mnull,mf))-sysf, atol = 1.e-7) == 0 &&
      info.nrank == 2 && info.rdeg == [0, 0] && info.nr == 3 &&
      all(real(gpole(sysx)) .< sdeg*(1-eps()^(1/info.nr)))

@time sysx, info, sysgen = grsol(sysg, sysf; mindeg = true, solgen = true)
isysg = 1:m; isysf = m+1:m+mf;
mnull = m - info.nrank; i0 = 1:mf; inull = mf+1:mf+mnull
@test gnrank(sysg*sysx-sysf, atol=1.e-7) == 0 &&
      gnrank(sysg*sysgen[:,i0]-sysf, atol = 1.e-7) == 0 &&
      gnrank(sysg*sysgen[:,inull], atol = 1.e-7) == 0 &&
      gnrank(sysg*(sysgen[:,i0]+sysgen[:,inull]*rdss(3,mnull,mf))-sysf, atol = 1.e-7) == 0 &&
      info.nrank == 2 && info.rdeg == [0, 0] && info.nr == 3

# Gao & Antsaklis (1989)
s = Polynomial([0, 1],'s');
Gn = [s-1 s-1];
Gd = [s*(s+1) s*(s+2)];
Fn = [s-1 s-1];
Fd = [(s+1)*(s+3) (s+1)*(s+4)];
m = 2; mf = 2;
sysg = dss(Gn,Gd,minimal = true, atol = 1.e-7);
sysf = dss(Fn,Fd,minimal = true, atol = 1.e-7);

@time sysx, info, sysgen = grsol(sysg, sysf; mindeg = true, solgen = true)
mnull = m - info.nrank; i0 = 1:mf; inull = mf+1:mf+mnull
@test gnrank(sysg*sysx-sysf, atol=1.e-7) == 0 &&
      gnrank(sysg*sysgen[:,i0]-sysf, atol = 1.e-7) == 0 &&
      gnrank(sysg*sysgen[:,inull], atol = 1.e-7) == 0 &&
      gnrank(sysg*(sysgen[:,i0]+sysgen[:,inull]*rdss(3,mnull,mf))-sysf, atol = 1.e-7) == 0 &&
      info.nrank == 1 && info.rdeg == [0, 0] && info.nr == 1

s = Polynomial([0, 1],'s');
g = [s^2+s+1 4*s^2+3*s+2 2*s^2-2;
    s 4*s-1 2*s-2;
    s^2 4*s^2-s 2*s^2-2*s];
f = g[:,1:2];
m = 3; mf = 2;
sysg = dss(g,minimal = true, atol = 1.e-7);
sysf = dss(f,minimal = true, atol = 1.e-7);
@time sysx, info, sysgen = grsol(sysg, sysf; atol = 1.e-7, mindeg = true, solgen = true)
mnull = m - info.nrank; i0 = 1:mf; inull = mf+1:mf+mnull
@test gnrank(sysg*sysx-sysf, atol=1.e-7) == 0 &&
      gnrank(sysg*sysgen[:,i0]-sysf, atol = 1.e-7) == 0 &&
      gnrank(sysg*sysgen[:,inull], atol = 1.e-7) == 0 &&
      gnrank(sysg*(sysgen[:,i0]+sysgen[:,inull]*rdss(3,mnull,mf))-sysf, atol = 1.e-7) == 0 &&
      info.nrank == 2 && info.rdeg == [0, 0] && info.nr == 0

z = Polynomial([0, 1],'s');
gn = [z^4-z^3/2-16*z^2-29/2*z+18 z^4+5*z^3-z^2-11*z+6 11/2*z^3+15*z^2+7/2*z-12
    -3*z^2+12 z^3-z^2-4*z+4 z^3+2*z^2-4*z-8;
    z^4-z^3/2-19*z^2-23/2*z+24 z^4+6*z^3-3*z^2-12*z+8 13/2*z^3+16*z^2-z/2-16];
gd = [z^4+5/2*z^3+2*z^2+z/2 z^4+5/2*z^3+2*z^2+z/2 z^4+5/2*z^3+2*z^2+z/2;
z^4+5/2*z^3+2*z^2+z/2 z^4+5/2*z^3+2*z^2+z/2 z^4+5/2*z^3+2*z^2+z/2;
z^4+5/2*z^3+2*z^2+z/2 z^4+5/2*z^3+2*z^2+z/2 z^4+5/2*z^3+2*z^2+z/2];
fn = gn[:,1:2];
fd = gd[:,1:2];
m = 3; mf = 2;
sysg = dss(gn,gd,minimal = true, atol = 1.e-7);
sysf = dss(fn,fd,minimal = true, atol = 1.e-7);

@time sysx, info, sysgen = grsol(sysg, sysf; atol = 1.e-7, mindeg = true, solgen = true)
mnull = m - info.nrank; i0 = 1:mf; inull = mf+1:mf+mnull
@test gnrank(sysg*sysx-sysf, atol=1.e-7) == 0 &&
      gnrank(sysg*sysgen[:,i0]-sysf, atol = 1.e-7) == 0 &&
      gnrank(sysg*sysgen[:,inull], atol = 1.e-7) == 0 &&
      gnrank(sysg*(sysgen[:,i0]+sysgen[:,inull]*rdss(3,mnull,mf))-sysf, atol = 1.e-7) == 0 &&
      info.nrank == 2 && info.rdeg == [0, 0] && info.nr == 0

fast = true; Ty = Complex{Float64}; Ty = Float64
mgf = 6; n = 5; p = 3; mmax = 3; mfmax = 2; mgfmax = mmax+mfmax;
m = 3; mf = 2;
mindeg = false;

# random examples
for Ty in (Float64, Complex{Float64})

#sysgf = rss(n,p,mgf,T = Ty,disc=false);

for fast in (true, false)

for mindeg in (true, false)

      # continuous, standard
      sysg = rss(n,p,m,T = Ty, disc=false);
      sysx0 = rss(3,m,mf,T = Ty, disc=false);
      sysf = sysg*sysx0;
      @time sysx, info, sysgen = grsol(sysg, sysf; atol = 1.e-7, mindeg = mindeg, solgen = true, fast = fast)
      mnull = m - info.nrank; i0 = 1:mf; inull = mf+1:mf+mnull
      @test gnrank(sysg*sysx-sysf, atol = 1.e-7) == 0 &&
            gnrank(sysg*sysgen[:,i0]-sysf, atol = 1.e-7) == 0 &&
            gnrank(sysg*sysgen[:,inull], atol = 1.e-7) == 0 &&
            gnrank(sysg*(sysgen[:,i0]+sysgen[:,inull]*rdss(3,mnull,mf))-sysf, atol = 1.e-7) == 0

      # discrete, descriptor
      sysg = rdss(n,p,m,T = Ty, disc=true);
      sysx0 = rss(3,m,mf,T = Ty, disc=true);
      sysf = sysg*sysx0;
      @time sysx, info, sysgen = grsol(sysg, sysf; atol = 1.e-7, mindeg = mindeg, solgen = true, fast = fast)
      mnull = m - info.nrank; i0 = 1:mf; inull = mf+1:mf+mnull
      @test gnrank(sysg*sysx-sysf, atol = 1.e-7) == 0 &&
            gnrank(sysg*sysgen[:,i0]-sysf, atol = 1.e-7) == 0 &&
            gnrank(sysg*sysgen[:,inull], atol = 1.e-7) == 0 &&
            gnrank(sysg*(sysgen[:,i0]+sysgen[:,inull]*rdss(3,mnull,mf,disc=true))-sysf, atol = 1.e-7) == 0

      # discrete, polynomial
      sysg = rdss(n,p,m,T = Ty, disc=true);
      sysx0 = rdss(0,m,mf,T = Ty, disc=true,id=[3*ones(Int,1);2*ones(Int,1)]);
      sysf = sysg*sysx0;
      @time sysx, info, sysgen = grsol(sysg, sysf; atol = 1.e-7, mindeg = mindeg, solgen = true, fast = fast)
      mnull = m - info.nrank; i0 = 1:mf; inull = mf+1:mf+mnull
      @test gnrank(sysg*sysx-sysf, atol = 1.e-7) == 0 &&
            gnrank(sysg*sysgen[:,i0]-sysf, atol = 1.e-7) == 0 &&
            gnrank(sysg*sysgen[:,inull], atol = 1.e-7) == 0 &&
            gnrank(sysg*(sysgen[:,i0]+sysgen[:,inull]*rdss(3,mnull,mf, disc=true))-sysf, atol = 1.e-7) == 0

      # continuous, improper
      sysg = rdss(n,p,m,T = Ty, disc=false, id=[2*ones(Int,1);2*ones(Int,1)]);
      sysx0 = rdss(3,m,mf,T = Ty, disc=false,id=[3*ones(Int,1);2*ones(Int,1)]);
      sysf = sysg*sysx0;
      @time sysx, info, sysgen = grsol(sysg, sysf; atol = 1.e-7, mindeg = mindeg, solgen = true, fast = fast)
      #@time sysx, info, sysgen = grsol(gir([sysg sysf],atol=1.e-7), mf; atol = 1.e-6, mindeg = mindeg, solgen = true, fast = fast)
      mnull = m - info.nrank; i0 = 1:mf; inull = mf+1:mf+mnull
      @test gnrank(sysg*sysx-sysf, atol = 1.e-7) == 0 &&
            gnrank(sysg*sysgen[:,i0]-sysf, atol = 1.e-7) == 0 &&
            gnrank(sysg*sysgen[:,inull], atol = 1.e-7) == 0 &&
            gnrank(sysg*(sysgen[:,i0]+sysgen[:,inull]*rdss(3,mnull,mf))-sysf, atol = 1.e-7) == 0

for m = p:mgfmax
      mf = mgfmax-m
      # continuous, standard
      sysgf = rss(n,p,mgfmax,T = Ty, disc=false);
      @time sysx, info, sysgen = grsol(sysgf, mf; atol = 1.e-7, mindeg = mindeg, solgen = true, fast = fast)
      isysg = 1:m; isysf = m+1:m+mf;
      mnull = m - info.nrank; i0 = 1:mf; inull = mf+1:mf+mnull
      @test gnrank(sysgf[:,isysg]*sysx-sysgf[:,isysf], atol = 1.e-7) == 0 &&
            gnrank(sysgf[:,isysg]*sysgen[:,i0]-sysgf[:,isysf], atol = 1.e-7) == 0 &&
            gnrank(sysgf[:,isysg]*sysgen[:,inull], atol = 1.e-7) == 0 &&
            gnrank(sysgf[:,isysg]*(sysgen[:,i0]+sysgen[:,inull]*rdss(3,mnull,mf))-sysgf[:,isysf], atol = 1.e-7) == 0

      # discrete, descriptor
      sysgf = rdss(n,p,mgfmax,T = Ty, disc=true);
      @time sysx, info, sysgen = grsol(sysgf, mf; atol = 1.e-7, mindeg = mindeg, solgen = true, fast = fast)
      isysg = 1:m; isysf = m+1:m+mf;
      mnull = m - info.nrank; i0 = 1:mf; inull = mf+1:mf+mnull
      @test gnrank(sysgf[:,isysg]*sysx-sysgf[:,isysf], atol = 1.e-7) == 0 &&
            gnrank(sysgf[:,isysg]*sysgen[:,i0]-sysgf[:,isysf], atol = 1.e-7) == 0 &&
            gnrank(sysgf[:,isysg]*sysgen[:,inull], atol = 1.e-7) == 0 &&
            gnrank(sysgf[:,isysg]*(sysgen[:,i0]+sysgen[:,inull]*rdss(3,mnull,mf, disc=true))-sysgf[:,isysf], atol = 1.e-7) == 0

      # discrete, polynomial
      sysgf = rdss(0,p,mgfmax,T = Ty, disc=false,id=[3*ones(Int,1);2*ones(Int,1)]);
      @time sysx, info, sysgen = grsol(sysgf, mf; atol = 1.e-7, mindeg = mindeg, solgen = true, fast = fast)
      isysg = 1:m; isysf = m+1:m+mf;
      mnull = m - info.nrank; i0 = 1:mf; inull = mf+1:mf+mnull
      @test gnrank(sysgf[:,isysg]*sysx-sysgf[:,isysf], atol = 1.e-7) == 0 &&
            gnrank(sysgf[:,isysg]*sysgen[:,i0]-sysgf[:,isysf], atol = 1.e-7) == 0 &&
            gnrank(sysgf[:,isysg]*sysgen[:,inull], atol = 1.e-7) == 0 &&
            gnrank(sysgf[:,isysg]*(sysgen[:,i0]+sysgen[:,inull]*rdss(3,mnull,mf))-sysgf[:,isysf], atol = 1.e-7) == 0

      # continuous, improper
      sysgf = rdss(n,p,mgfmax,T = Ty, disc=true,id=[3*ones(Int,1);2*ones(Int,1)]);
      @time sysx, info, sysgen = grsol(sysgf, mf; atol = 1.e-7, mindeg = mindeg, solgen = true, fast = fast)
      isysg = 1:m; isysf = m+1:m+mf;
      mnull = m - info.nrank; i0 = 1:mf; inull = mf+1:mf+mnull
      @test gnrank(sysgf[:,isysg]*sysx-sysgf[:,isysf], atol = 1.e-7) == 0 &&
            gnrank(sysgf[:,isysg]*sysgen[:,i0]-sysgf[:,isysf], atol = 1.e-7) == 0 &&
            gnrank(sysgf[:,isysg]*sysgen[:,inull], atol = 1.e-7) == 0 &&
            gnrank(sysgf[:,isysg]*(sysgen[:,i0]+sysgen[:,inull]*rdss(3,mnull,mf, disc=true))-sysgf[:,isysf], atol = 1.e-7) == 0
end

end # mindeg
end # fast
end # Ty

end # grsol

@testset "glsol" begin

p = 1; pf = 0
sysgf = dss(ones(1,1));
@time sysx, info, sysgen = glsol(sysgf, pf, solgen = true)
isysg = 1:p; isysf = p+1:p+pf;
pnull = p-info.nrank; i0 = 1:pf; inull = pf+1:pf+pnull
@test gnrank(sysx*sysgf[isysg,:]-sysgf[isysf,:]) == 0 &&
      gnrank(sysgen[i0,:]*sysgf[isysg,:]-sysgf[isysf,:], atol = 1.e-7) == 0 &&
      gnrank(sysgen[inull,:]*sysgf[isysg,:], atol = 1.e-7) == 0 &&
      gnrank((sysgen[i0,:]+rdss(3,pf,pnull)*sysgen[inull,:])*sysgf[isysg,:]-sysgf[isysf,:]) == 0 &&
      info.nrank == 1 && info.rdeg == Int[]

p = 1; pf = 0
sysg = dss(ones(p,p)); sysf = dss(ones(pf,p));
@time sysx, info, sysgen = glsol(sysg, sysf, solgen = true)
pnull = p-info.nrank; i0 = 1:pf; inull = pf+1:pf+pnull
@test gnrank(sysx*sysg-sysf) == 0 &&
      gnrank(sysgen[i0,:]*sysg-sysf, atol = 1.e-7) == 0 &&
      gnrank(sysgen[inull,:]*sysg, atol = 1.e-7) == 0 &&
      gnrank((sysgen[i0,:]+rdss(3,pf,pnull)*sysgen[inull,:])*sysg-sysf) == 0 &&
      info.nrank == 1 && info.rdeg == Int[]

p = 5; pf = 2
sysgf = dss(rand(p+pf,p));
@time sysx, info, sysgen = glsol(sysgf, pf, solgen = true)
isysg = 1:p; isysf = p+1:p+pf;
pnull = p-info.nrank; i0 = 1:pf; inull = pf+1:pf+pnull
@test gnrank(sysx*sysgf[isysg,:]-sysgf[isysf,:], atol = 1.e-7) == 0 &&
      gnrank(sysgen[i0,:]*sysgf[isysg,:]-sysgf[isysf,:], atol = 1.e-7) == 0 &&
      gnrank(sysgen[inull,:]*sysgf[isysg,:], atol = 1.e-7) == 0 &&
      gnrank((sysgen[i0,:]+rdss(3,pf,pnull)*sysgen[inull,:])*sysgf[isysg,:]-sysgf[isysf,:], atol = 1.e-7) == 0 &&
      info.nrank == 5 && info.rdeg == [0, 0] && info.nl == 0

p = 5; pf = 2
sysg = dss(rand(p,p)); sysf = dss(rand(pf,p));
@time sysx, info, sysgen = glsol(sysg, sysf, solgen = true)
pnull = p-info.nrank; i0 = 1:pf; inull = pf+1:pf+pnull
@test gnrank(sysx*sysg-sysf, atol = 1.e-7) == 0 &&
      gnrank(sysgen[i0,:]*sysg-sysf, atol = 1.e-7) == 0 &&
      gnrank(sysgen[inull,:]*sysg, atol = 1.e-7) == 0 &&
      gnrank((sysgen[i0,:]+rdss(3,pf,pnull)*sysgen[inull,:])*sysg-sysf, atol = 1.e-7) == 0 &&
      info.nrank == 5 && info.rdeg == [0, 0] && info.nl == 0

m = 5; p = 8; pf = 2
sysg = dss(rand(p,m)); sysf = dss(rand(pf,m));
@time sysx, info, sysgen = glsol(sysg, sysf, solgen = true)
pnull = p-info.nrank; i0 = 1:pf; inull = pf+1:pf+pnull
@test gnrank(sysx*sysg-sysf, atol = 1.e-7) == 0 &&
      gnrank(sysgen[i0,:]*sysg-sysf, atol = 1.e-7) == 0 &&
      gnrank(sysgen[inull,:]*sysg, atol = 1.e-7) == 0 &&
      gnrank((sysgen[i0,:]+rdss(3,pf,pnull)*sysgen[inull,:])*sysg-sysf, atol = 1.e-7) == 0 &&
      info.nrank == 5 && info.rdeg == [0, 0]  && info.nl == 0

# Wang and Davison Example (1973)
s = Polynomial([0, 1],'s');
gn = [ s+1 s+2; s+3 s^2+2*s; s^2+3*s 0 ];
gd = [s^2+3*s+2 s^2+3*s+2; s^2+3*s+2 s^2+3*s+2; s^2+3*s+2 s^2+3*s+2];
p = 3; pf = 2;
sysg = dss(gn,gd,minimal = true, atol = 1.e-7);
sysf = dss([1 0;0 1.]);

@time sysx, info, sysgen = glsol(sysg, sysf, solgen = true)
pnull = info.pl; i0 = 1:pf; inull = pf+1:pf+pnull
@test gnrank(sysx*sysg-sysf, atol = 1.e-7) == 0 &&
      gnrank(sysgen[i0,:]*sysg-sysf, atol = 1.e-7) == 0 &&
      gnrank(sysgen[inull,:]*sysg, atol = 1.e-7) == 0 &&
      gnrank((sysgen[i0,:]+rdss(3,pf,pnull)*sysgen[inull,:])*sysg-sysf, atol = 1.e-7) == 0 &&
      info.nrank == 2 && info.rdeg == [0, 0] && info.nl == 3

evals=[-1, -2, -3];
@time sysx, info, sysgen = glsol(sysg, sysf, poles = evals, solgen = true)
pnull = info.pl; i0 = 1:pf; inull = pf+1:pf+pnull
@test gnrank(sysx*sysg-sysf, atol = 1.e-7) == 0 &&
      gnrank(sysgen[i0,:]*sysg-sysf, atol = 1.e-7) == 0 &&
      gnrank(sysgen[inull,:]*sysg, atol = 1.e-7) == 0 &&
      gnrank((sysgen[i0,:]+rdss(3,pf,pnull)*sysgen[inull,:])*sysg-sysf, atol = 1.e-7) == 0 &&
      info.nrank == 2 && info.rdeg == [0, 0] && info.nl == 3 && info.pl == 1 &&
      sort(gpole(sysx)) ≈ sort(evals)

sdeg = -1;
@time sysx, info, sysgen = glsol(sysg, sysf, sdeg = sdeg, solgen = true)
pnull = info.pl; i0 = 1:pf; inull = pf+1:pf+pnull
@test gnrank(sysx*sysg-sysf, atol = 1.e-7) == 0 &&
      gnrank(sysgen[i0,:]*sysg-sysf, atol = 1.e-7) == 0 &&
      gnrank(sysgen[inull,:]*sysg, atol = 1.e-7) == 0 &&
      gnrank((sysgen[i0,:]+rdss(3,pf,pnull)*sysgen[inull,:])*sysg-sysf, atol = 1.e-7) == 0 &&
      info.nrank == 2 && info.rdeg == [0, 0] && info.nl == 3 && info.pl == 1 &&
      all(real(gpole(sysx)) .< sdeg*(1-eps()^(1/info.nl)))

@time sysx, info, sysgen = glsol(sysg, sysf; mindeg = true, solgen = true)
pnull = info.pl; i0 = 1:pf; inull = pf+1:pf+pnull
@test gnrank(sysx*sysg-sysf, atol = 1.e-7) == 0 &&
      gnrank(sysgen[i0,:]*sysg-sysf, atol = 1.e-7) == 0 &&
      gnrank(sysgen[inull,:]*sysg, atol = 1.e-7) == 0 &&
      gnrank((sysgen[i0,:]+rdss(3,pf,pnull)*sysgen[inull,:])*sysg-sysf, atol = 1.e-7) == 0 &&
      info.nrank == 2 && info.rdeg == [0, 0] && info.nl == 3 && info.pl == 1

end # glsol

end # linsol

end # module


text:
#Text provided under a Creative Commons Attribution license, CC-BY.  All code is made available under the FSF-approved BSD-3 license.  (c) Lorena A. Barba, Gilbert F. Forsyth 2017. Thanks to NSF for support via CAREER award #1149784.
#[@LorenaABarba](https://twitter.com/LorenaABarba)
# Translation to Julia by [@miguelraz_](https://twitter.com/miguelraz_)

## 12 steps to Navier–Stokes
# ======
# ***

# Hello! Welcome to the **12 steps to Navier–Stokes**. This is a practical module that is used in the beginning of an interactive
# Computational Fluid Dynamics (CFD) course taught by [Prof. Lorena Barba](http://lorenabarba.com) since Spring 2009 at Boston University.
# The course assumes only basic programming knowledge (in any language) and of course some
# foundation in partial differential equations and fluid mechanics. The practical module was inspired by the
# ideas of Dr. Rio Yokota, who was a post-doc in Barba's lab, and has been refined by Prof. Barba and
# her students over several semesters teaching the course. The course is taught entirely using Python (but this version is in Julia)
# and students who don't know Julia just learn as we work through the module.

# This [Pluto notebook](https://github.com/fonsp/Pluto.jl) will lead you through the first step of programming
# your own Navier–Stokes solver in Julia from the ground up.  We're going to dive right in.
# Don't worry if you don't understand everything that's happening at first, we'll cover it in
# detail as we move forward and you can support your learning with the videos of
# [Prof. Barba's lectures on YouTube](http://www.youtube.com/playlist?list=PL30F4C5ABCE62CB61).

# For best results, after you follow this notebook, prepare your own code for Step 1, either as a Julia script or in a clean Pluto notebook.

# To execute this Notebook, we assume you have invoked the notebook server using: `using Pluto; Pluto.run()`.

## Step 1: 1-D Linear Convection
# -----
# ***

# The 1-D Linear Convection equation is the simplest, most basic model that can be used to learn something about CFD.
# It is surprising that this little equation can teach us so much! Here it is:

# $$\frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0$$

# With given initial conditions (understood as a *wave*), the equation represents the propagation of that
# initial *wave* with speed $c$, without change of shape. Let the initial condition be $u(x,0)=u_0(x)$.
# Then the exact solution of the equation is $u(x,t)=u_0(x-ct)$.

# We discretize this equation in both space and time, using the Forward Difference scheme for the time derivative
# and the Backward Difference scheme for the space derivative. Consider discretizing the spatial coordinate
# $x$ into points that we index from $i=1$ to $N$, and stepping in discrete time intervals of size $\Delta t$.

# From the definition of a derivative (and simply removing the limit), we know that:

# $$\frac{\partial u}{\partial x}\approx \frac{u(x+\Delta x)-u(x)}{\Delta x}$$

# Our discrete equation, then, is:

# $$\frac{u_i^{n+1}-u_i^n}{\Delta t} + c \frac{u_i^n - u_{i-1}^n}{\Delta x} = 0 $$

# Where $n$ and $n+1$ are two consecutive steps in time, while $i-1$ and $i$ are two neighboring points of the discretized $x$ coordinate. If there are given initial conditions, then the only unknown in this discretization is $u_i^{n+1}$.  We can solve for our unknown to get an equation that allows us to advance in time, as follows:

# $$u_i^{n+1} = u_i^n - c \frac{\Delta t}{\Delta x}(u_i^n-u_{i-1}^n)$$

# Now let's try implementing this in Julia.

# We'll start by importing a few libraries to help us out.

# * `Plots` is a 2D plotting library that we will use to plot our results

# Remember: comments in python are denoted by the pound sign
using Plots                        #here we load numpy

# Now let's define a few variables; we want to define an evenly spaced grid of points within
# a spatial domain that is 2 units of length wide, i.e., $x_i\in(0,2)$.
# We'll define a variable `nx`, which will be the number of grid points we want and `dx` will be the distance between any pair of adjacent grid points.

nx = 41  # try changing this number from 41 to 81 and Run All ... what happens?
dx = 2 / (nx-1)
nt = 25    #nt is the number of timesteps we want to calculate
dt = .025  #dt is the amount of time each timestep covers (delta t)
c = 1      #assume wavespeed of c = 1

# We also need to set up our initial conditions. The initial velocity $u_0$ is given as
# $u = 2$ in the interval $0.5 \leq x \leq 1$  and $u = 1$ everywhere else in $(0,2)$ (i.e., a hat function).

# Here, we use the function `ones()` defining an array which is `nx` elements long with every value equal to 1.

u = ones(nx)      #makes a vector of size nx filled with only 1.0s
u[Int.((.5/dx):(1/dx))] .= 2 # #setting u = 2 between 0.5 and 1 as per our I.C.s

println(u)

# Now let's take a look at those initial conditions using a Plots.jl plot.
# We've imported the `Plots` plotting library `pyplot` and the plotting function is called `plot`,
# so we'll call `pyplot.plot`. To learn about the myriad possibilities of Plots.jl,
# explore the [Documentation](https://docs.juliaplots.org/) of example plots.

# Here, we use the syntax for a simple 2D plot: `plot(x,y)`, where the `x` values are evenly distributed grid points:

using Plots
Plots.pyplot()
pyplot.plot(range(start = 0, stop = 2, length = nx), u);

# Why doesn't the hat function have perfectly straight sides? Think for a bit.

# Now it's time to implement the discretization of the convection equation using a finite-difference scheme.

# For every element of our array `u`, we need to perform the operation $u_i^{n+1} = u_i^n - c \frac{\Delta t}{\Delta x}(u_i^n-u_{i-1}^n)$

# We'll store the result in a new (temporary) array `un`, which will be the solution $u$ for the next time-step.
# We will repeat this operation for as many time-steps as we specify and then we can see how far the wave has convected.

# We first initialize our placeholder array `un` to hold the values we calculate for the $n+1$ timestep, using once again the Julia function `ones()`.

# Then, we may think we have two iterative operations: one in space and one in time (we'll learn differently later),
# so we'll start by nesting one loop inside the other. Note the use of the nifty `enumerate` function.
# When we write: `for i in nx` we will iterate through the `u` array.

un = ones(nx) #initialize a temporary array

for n in nt  #loop for values of n from 0 to nt, so it will run nt times
    un = copy(u) ##copy the existing values of u into un
    for (elem, i) in enumerate(nx)
    # for i in 1:length(nx) # Try commenting the previous line and uncommenting this one to see what happens!
        u[i] = un[i] - c * dt / dx * (un[i] - un[i-1])
    end
end

# **Note**—We will learn later that the code as written above is quite inefficient, and there are better ways to write this, Julia-style. But let's carry on.

# Now let's try plotting our `u` array after advancing in time.
plot(range(start = 0, end = 2, length = nx), u)

# OK! So our hat function has definitely moved to the right, but it's no longer a hat.  **What's going on?**

## Learn More
# -----
# ***

# For a more thorough explanation of the finite-difference method, including topics like the truncation error,
# order of convergence and other details, watch **Video Lessons 2 and 3** by Prof. Barba on YouTube.

# TODO
# from IPython.display import YouTubeVideo
# YouTubeVideo('iz22_37mMkk')

# TODO
# YouTubeVideo('xq9YTcv-fQg')
#
# For a careful walk-through of the discretization of the linear convection equation with finite differences
# (and also the following steps, up to Step 4), watch **Video Lesson 4** by Prof. Barba on YouTube.

# TODO
# YouTubeVideo('y2WaK7_iMRI')

## Last but not least

# **Remember** to rewrite Step 1 as a fresh Julia script or in *your own* Pluto notebook and
# then experiment by changing the discretization parameters.
# Once you have done this, you will be ready for [Step 2](./02_Step_2.jl).

ext: jl

max_stars_repo_name: mkitti/CFDJulia

max_issues_repo_name: mkitti/CFDJulia

max_forks_repo_name: mkitti/CFDJulia


text:
using DashBootstrapComponents

input_group = html_div([
    dbc_inputgroup([dbc_inputgrouptext("Large"), dbc_input()], size = "lg"),
    html_br(),
    dbc_inputgroup([dbc_inputgrouptext("Default"), dbc_input()]),
    html_br(),
    dbc_inputgroup([dbc_inputgrouptext("Small"), dbc_input()], size = "sm"),
]);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
# TODO - overflow
# TODO - fast mult as http://maths-people.anu.edu.au/~brent/pd/rpb232.pdf

# a pragmatic library for doing modular arithmetic.

# the aim here is not to encapsulate a large amount of theory, or to
# describe the theoretical relationships between different structures,
# but to enable arithmetic on various types, motivated largely by the
# practical needs of crypto code.

# incomplete; pull requests welcome

# much thanks to Andreas Noack Jensen
# https://groups.google.com/d/msg/julia-users/Ui977brdzAU/u4rWiDeJv-MJ
# https://gist.github.com/andreasnoackjensen/9132511 (but since this
# changes things significanty, mistakes are mine...)

module IntModN

import Base: show, showcompact, zero, one, inv, real, abs, convert,
       promote_rule, length, getindex, setindex!, start, done, next,
       rand, rand!, print, map, leading_zeros, divrem, endof, bits,
       div, rem, mod, trailing_zeros, /, ==, <=, <, -, ^, +, *, &,
       |, $, %, <<, >>>, >>
import Base.Random.AbstractRNG

using AutoTypeParameters

export ZModN, ZField, ZRing, ZR, ZF, GF2, @zring, @zfield,
       ZPoly, ZP, X, P, order, modulus, factor, itype,
       GF2Poly, GF2P, GF2X,
       FModN, FRing, FR,
       extended_euclidean

# we generally support promotion into types here (from integers) and
# conversion back.  but that's all.  users can add other promotion /
# conversion if required.

# using Integer as the base type here as: (1) some numeric type seems
# to be necessary to access many generic numerical functions; (2) need
# Real to use comparisons; (3) have integer-like accuracy behaviour
# (no need to worry about rounding etc).

# but we need an additional type so that we can intercept things like
# automatic promotion to floats from integers.

abstract Residue <: Integer

/{I<:Integer,M<:Residue}(i::I, m::M) = convert(M, i) / m
/{I<:Integer,M<:Residue}(m::M, i::I) = m / convert(M, i)
inv{M<:Residue}(m::M) = error("no inverse defined")

# --- integers modulo n

# J is the type of N.  this may seem superfluous, but N itself is frozen, so its type is not always available (see AutoTypeParameters)

abstract ZModN{J<:Integer, N, I<:Integer} <: Residue

# we need two types here because prime moduli have a faster inverse
# (via euler's theorem) (note - we do NOT test for primality).

# you can construct these directly, but code assumes that the integer is
# already reduced mod N, that capacities are ok, etc.

immutable ZRing{J<:Integer, N, I<:Integer} <: ZModN{J,N,I}
    i::I
    ZRing(i) = (@assert isa(i, J); new(i))
end

# this assumes N is prime and uses a faster inverse
immutable ZField{J<:Integer, N, I<:Integer} <: ZModN{J,N,I}
    i::I
    ZField(i) = (@assert isa(i, J); new(i))
end

# but typically you would use one of these constructors (ZF or ZR)

function prepare_z{I<:Integer}(n::Int, i::I)
    n > 0 || error("modulus ($n) too small")
    n <= typemax(typeof(i)) || error("modulus ($n) too large for $(typeof(i))")
    mod(i, n)
end

for (f, Z) in ((:ZR, :ZRing), (:ZF, :ZField))
    @eval $f{J<:Integer, I<:Integer}(n, ::Type{J}, i, ::Type{I}) = $Z{J, freeze(convert(J, n)), I}(prepare_z(n, convert(I, i)))
    @eval $f{I<:Integer}(n::Int, i, ::Type{I}) = $Z{typeof(n), freeze(n), I}(prepare_z(n, convert(I, i)))
    @eval $f(n::Int, i::Integer) = $Z{typeof(n), freeze(n), typeof(i)}(prepare_z(n, i))
    @eval $f(n::Integer, i::Integer) = $Z{typeof(n), freeze(n), typeof(i)}(mod(i, n))
    @eval $f(n) = i::Integer -> $f(n, i)  # used to construct constructors
end
GF2 = ZF(2)

modulus{J<:Integer, N, I<:Integer}(::Type{ZModN{J, N, I}}) = convert(J, thaw(eval, N))
modulus{T<:ZModN}(::Type{T}) = modulus(super(T))  # jameson type chaina
modulus{T<:ZModN}(::T) = modulus(T)
order{T<:ZModN}(::Type{T}) = modulus(T)
order{T<:ZModN}(::T) = modulus(T)
itype{J, N, I<:Integer}(::Type{ZModN{J, N, I}}) = I
itype{T<:ZModN}(::Type{T}) = itype(super(T))  # jameson type chain
itype{T<:ZModN}(::T) = itype(T)
ntype{J, N, I<:Integer}(::Type{ZModN{J, N, I}}) = J
ntype{T<:ZModN}(::Type{T}) = ntype(super(T))  # jameson type chain
ntype{T<:ZModN}(::T) = ntype(T)

for (f, Z) in ((:ZR, :ZRing), (:ZF, ZField))
    # julia's comversion is pretty damn dumb at times
    @eval convert{J<:Integer,N,I<:Integer}(::Type{$Z{J,N,I}}, a::$Z{J,N,I}) = a
    # this promotion needed because linalg code has direct compariosn with 0
    @eval promote_rule{J<:Integer, N, I<:Integer}(::Type{$Z{J,N,I}}, ::Type{Int}) = $Z{J,N,I}
    @eval convert{J<:Integer,N,I<:Integer}(::Type{$Z{J,N,I}}, i::Integer) = $f(convert(J, N), convert(I, i))
end
# conversion to int, but no promotion
for T in (:Int, :UInt)
    @eval convert(::Type{$T}, z::ZModN) = convert($T, z.i)
end

# all the methods that make these numbers

for Z in (:ZRing, :ZField)
    @eval zero{J,N,I}(::Type{$Z{J,N,I}}) = $Z{J,N,I}(zero(I))
    @eval one{J,N,I}(::Type{$Z{J,N,I}}) = $Z{J,N,I}(one(I))
end
zero{Z<:ZModN}(z::Z) = zero(Z)
one{Z<:ZModN}(z::Z) = one(Z)

showcompact{J,N,I}(io::IO, z::ZModN{J,N,I}) = showcompact(io, z.i)
print{J,N,I}(io::IO, z::ZModN{J,N,I}) = print(io, "$(z.i) mod $N")
show{J,N,I<:Integer}(io::IO, z::ZField{J,N,I}) = print(io, "ZField{$J,$N,$I}($(z.i))")
show{J,N,I<:Integer}(io::IO, z::ZRing{J,N,I}) = print(io, "ZRing{$J,$N,$I}($(z.i))")

rand_{J,N,I<:Integer}(T::Type{ZRing{J,N,I}}) = ZR(N, J, rand(one(I):convert(I,modulus(T))), I)
rand_{J,N,I<:Integer}(T::Type{ZField{J,N,I}}) = ZF(N, J, rand(one(I):convert(I,modulus(T))), I)
function rand!{T<:ZModN}(::Type{T}, a::AbstractArray)
    for i in 1:length(a)
        a[i] = rand_(T)
    end
    a
end

rand(::Type{ZModN}, dims...) = error("use more specific type")

function rand{T<:ZModN}(::Type{T}, dims...)
    if length(dims) > 0
        rand!(T, Array(T, dims))
    else
        rand_(T)
    end
end

if VERSION >= v"0.4-"
    rand{J<:Integer,N,I<:Integer}(rng::AbstractRNG, T::Type{ZRing{J,N,I}}) = ZR(N, J, rand(rng, one(I):convert(I,N)), I)
    rand{J<:Integer,N,I<:Integer}(rng::AbstractRNG, T::Type{ZField{J,N,I}}) = ZF(N, J, rand(rng, one(I):convert(I,N)), I)
    function rand!{T<:ZModN}(rng::AbstractRNG, ::Type{T}, a::AbstractArray)
        for i in 1:length(a)
            a[i] = rand(rng, T)
        end
        a
    end
    rand(rng::AbstractRNG, ::Type{ZModN}, dims::Int...) =
    error("use more specific type")
    rand{T<:ZModN}(rng::AbstractRNG, ::Type{T}, dims::Int...) =
    rand!(rng, T, Array(T, dims))
end

# equal N from types, so we don't need to test explicitly
=={T,N,I}(a::ZModN{T,N,I}, b::ZModN{T,N,I}) = a.i == b.i
<={T,N,I}(a::ZModN{T,N,I}, b::ZModN{T,N,I}) = a.i <= b.i
<{T,N,I}(a::ZModN{T,N,I}, b::ZModN{T,N,I}) = a.i < b.i

real{T,N,I}(a::ZModN{T,N,I}) = a
abs{T,N,I}(a::ZModN{T,N,I}) = a

liftz(c, i, f, a) = c(convert(i, mod(f(a.i), modulus(a))))
liftz(c, i, f, a, b) = c(convert(i, mod(f(a.i, b.i), modulus(a))))
for Z in (ZRing, ZField)
    @eval -{J,N,I}(a::$Z{J,N,I}) = liftz($Z{J,N,I}, I, -, a)
    @eval ^{J,N,I}(a::$Z{J,N,I}, p::Int) = liftz($Z{J,N,I}, I, x -> powermod(x, p, modulus(a)), a)
    for op in (:+, :-, :*, :&, :|, :$)
        @eval $op{J,N,I}(a::$Z{J,N,I}, b::$Z{J,N,I}) = liftz($Z{J,N,I}, I, $op, a, b)
    end
end

# http://en.wikipedia.org/wiki/Extended_Euclidean_algorithm#Residue_integers
function extended_euclidean{I<:Integer}(a::I, n::I)
    t::I, newt::I = zero(I), one(I)
    r::I, newr::I = n, a
    while newr != zero(I)
        q::I = div(r, newr)
        t, newt = newt, t - q * newt
        r, newr = newr, r - q * newr
    end
    r <= 1 || error("$a is not invertible mod $n")
    t < 0 ? t + n : t
end

inv{J,N,I}(a::ZRing{J,N,I}) = ZRing{J,N,I}(extended_euclidean(a.i, convert(I, modulus(a))))
inv{J,N,I}(a::ZField{J,N,I}) = a^(modulus(a)-2)

function /{T,N,I}(a::ZModN{T,N,I}, b::ZModN{T,N,I})
    if a == b
        one(a)
    elseif a == zero(a)
        a
    elseif b == zero(b)
        error("division by zero: $a/$b")
    else
        a * inv(b)
    end
end

# macros to rewrite literals mod n
# (these rewrite ALL ints, including args to literal ZF(...), etc)

macro zring(n::Int, code)
    :($(esc(rewrite_ints(n, x -> ZR(n, x), code))))
end

macro zfield(n::Int, code)
    :($(esc(rewrite_ints(n, x -> ZF(n, x), code))))
end

function rewrite_ints(n, c, val)
    if isa(val, Integer)
        :($(c(val)))
    else
        if isa(val, Expr) && val.head != :line
            val.args = [rewrite_ints(n, c, a) for a in val.args]
        end
        val
    end
end

# --- supertype for all integer polynomial representations

# I is the type that [] returns when indexing a coefficient
abstract Poly{I<:Integer} <: Residue

endof(a::Poly) = length(a)
start(::Poly) = 1
done(a::Poly, i) = i > length(a)
next(a::Poly, i) = (a[i], i+1)

length{P<:Poly}(a::P) = error("define a length() for type $P")
getindex{P<:Poly}(a::P, i) = error("define a getindex() for type $P")

# delegate modulus to the coefficient type

modulus{T}(::Poly{T}) = modulus(T)
modulus{T}(::Type{Poly{T}}) = modulus(T)   # is this useful?

# generic comparison

function cmp{P<:Poly}(a::P, b::P, eq)
    if length(a) < length(b)
        true
    elseif length(a) > length(b)
        false
    else
        for i in length(a):-1:1
            if a[i] < b[i]
                return true
            elseif a[i] > b[i]
                return false
            end
        end
        eq
    end
end
<={P<:Poly}(a::P, b::P) = cmp(a, b, true)
<{P<:Poly}(a::P, b::P) = cmp(a, b, false)

# exotic polynomials may require a suitable show_coeff and/or show_mod to
# display well

show_coeff(io::IO, i::Integer) = show(io, i)
show_coeff{T,N,I<:Integer}(io::IO, z::ZModN{T,N,I}) = show(io, convert(I, z))

show_mod{I<:Integer}(io::IO, p::Poly{I}) = nothing

function print_no_mod{I<:Integer}(io::IO, p::Poly{I})
    n = length(p)
    if n <= 0
        show_coeff(io, zero(I))
    else
        for j = n:-1:1
            pj = p[j]
            if pj != zero(I)
                if j == n
                    pj < 0 && print(io, "-")
                else
                    pj < 0 ? print(io," - ") : print(io," + ")
                end
                magpj = abs(pj)
                if j == 1 || magpj != one(I)
                    show_coeff(io, magpj)
                end
                if j > 1
                    print(io, :x)
                    if j > 2
                        print(io, '^', j-1)
                    end
                end
            end
        end
    end
end

function print{I<:Integer}(io::IO, p::Poly{I})
    print_no_mod(io, p)
    show_mod(io, p)
end

function show{I<:Integer}(io::IO, p::Poly{I})
    print(io, "ZP($I")
    for i in length(p):-1:1
        print(io, ",")
        show_coeff(io, p[i])
    end
    print(io, ")")
end

# --- arbitrary (but "integer") coeff polynomials

immutable ZPoly{I<:Integer} <: Poly{I}
    a::Vector{I}
    ZPoly(a) = length(a) == 0 || a[end] != zero(I) ? new(a) : error("zero leading coeff")
end

immutable FZP{I<:Integer}
    a::Tuple
end

# constructors

function prepare_p{N}(coeffs::Vector{N})
    start = endof(coeffs)
    while start > 0 && coeffs[start] == zero(N)
        start -= 1
    end
    start == endof(coeffs) ? coeffs : coeffs[1:start]
end

ZP() = error("provide at least one (zero?) coeff")
ZP(c::Function) = error("provide at least one (zero?) coeff")
# the vector constructors are in storage order
ZP(t::Type, coeffs::Vector) = ZP(t[map(c -> convert(t, c), coeffs)...])
ZP(c::Function, coeffs::Vector) = ZP(map(c, coeffs))
ZP{I<:Integer}(coeffs::Vector{I}) = ZPoly{I}(prepare_p(coeffs))
# the direct (tuple) constructors are in natural order
ZP(t::Type, coeffs...) = ZP(reverse(t[map(c -> convert(t, c), coeffs)...]))
ZP(c::Function, coeffs...) = ZP(reverse([map(c, coeffs)...]))
ZP{I<:Integer}(coeffs::I...) = ZP(reverse([coeffs...]))

# nice syntax for creating polynomials.
# use x=X(ZF(2)) or x=X(GF2) or x=X(ZField{2,Int}) and then x^2+3 etc
X(c::Function) = ZP([c(0), c(1)])
X{I<:Integer}(::Type{I}) = ZP([zero(I), one(I)])

# this promotion to write x^2+1 etc
convert{I<:Integer}(::Type{ZPoly{I}}, p::ZPoly{I}) = p
promote_rule{I<:Integer}(::Type{ZPoly{I}}, ::Type{Int}) = ZPoly{I}
convert{I<:Integer}(::Type{ZPoly{I}}, i::Integer) = ZP([convert(I, i)])

# can be accessed and modified (but don't do that in public) as arrays
getindex(a::ZPoly, i) = getindex(a.a, i)
setindex!{T,I<:Integer}(a::ZPoly{T}, v::I, i) = setindex!(a.a, convert(T, v), i)
length(a::ZPoly) = length(a.a)

show_mod{Z<:ZModN}(io::IO, p::ZPoly{Z}) = print(io, " mod $(modulus(Z))")
showcompact(io::IO, p::ZPoly) = showcompact(io, p.a)

# number methods

zero{T}(::Type{ZPoly{T}}) = ZPoly{T}(T[])
zero{Z<:ZPoly}(z::Z) = zero(Z)
one{T}(::Type{ZPoly{T}}) = ZPoly{T}([one(T)])
one{Z<:ZPoly}(z::Z) = one(Z)

=={T}(a::ZPoly{T}, b::ZPoly{T}) = a.a == b.a
-{T}(a::ZPoly{T}) = ZPoly{T}(-a.a)

function trailing_zeros{T}(p::ZPoly{T})
    n = 1
    while n < length(p) && p[n] == zero(T)
        n += 1
    end
    n - 1
end

# a and b same length; discard leading zeroes
function _truncate{T}(f, a::Vector{T}, b::Vector{T})
    ZERO, la = zero(T), length(a)
    copied = false
    for i in la:-1:1
        x::T = f(a[i], b[i])
        if x != ZERO && !copied
            result = Array(T, i)
            copied = true
        end
        if copied
            result[i] = x
        end
    end
    if !copied
        result = T[]
    end
    result
end

# big longer than small; copy extra values
function _skip{T}(f, big::Vector{T}, small::Vector{T})
    lb, ls = length(big), length(small)
    result = Array(T, lb)
    for i in lb:-1:(ls+1)
        result[i] = big[i]
    end
    for i in ls:-1:1
        result[i] = f(big[i], small[i])
    end
    result
end

function +{T}(a::ZPoly{T}, b::ZPoly{T})
    la, lb = length(a), length(b)
    if la == 0
        b
    elseif lb == 0
        a
    elseif la == lb
        ZPoly{T}(_truncate(+, a.a, b.a))
    elseif la > lb
        ZPoly{T}(_skip(+, a.a, b.a))
    else
        ZPoly{T}(_skip(+, b.a, a.a))
    end
end

function -{T}(a::ZPoly{T}, b::ZPoly{T})
    la, lb = length(a), length(b)
    if la == lb
        if a == b
            zero(ZPoly{T})
        else
            ZPoly{T}(_truncate(-, a.a, b.a))
        end
    elseif la > lb
        ZPoly{T}(_skip(-, a.a, b.a))
    else
        # in this case, b is strictly bigger than a, and leading
        # coeffs must be negated.  so there's no chance to truncate
        # and we cannot shift past the start.  so do it here.
        result = Array(T, lb)
        for i = lb:-1:(la+1)
            result[i] = -b.a[i]
        end
        for i = la:-1:1
            result[i] = a[i] - b[i]
        end
        ZPoly{T}(result)
    end
end

function *{T}(a::ZPoly{T}, b::ZPoly{T})
    big, small = length(a) > length(b) ? (a, b) : (b, a)
    if length(small) == 0  # no need to test big
        small
    elseif length(small) == 1 && small[1] == one(T)
        big
    else
        result = zeros(T, length(big) + length(small) - 1)
        for i in 1:length(small)
            z = small[i]  # enumerate here uses strange amount of memory
            if z != zero(T)
                for j in 1:length(big)
                    result[i+j-1] += big[j] * z
                end
            end
        end
        ZP(result)
    end
end

function _divrem{T}(a::Vector{T}, b::Vector{T})
    la, lb = length(a), length(b)
    lb > 0 || error("division by zero polynomial")
    if la == 0
        (a, a)
    elseif lb == 1 && b[1] == one(T)
        (a, T[])
    elseif lb > la
        (T[], a)
    elseif a == b
        ([one(T)], T[])
    else
        shift = la - lb  # non-negative
        rem, div = copy(a), zeros(T, shift + 1)
        for s in 0:shift
            factor = rem[end-s] / b[end]
            if factor != zero(T)
                div[end-s] = factor
                for i in 1:lb
                    rem[i+shift-s] -= factor * b[i]
                end
            end
        end
        div, rem
    end
end

function divrem{T}(a::ZPoly{T}, b::ZPoly{T})
    map(ZP, _divrem(a.a, b.a))  # ZP discards leading zeroes
end

div{T}(a::ZPoly{T}, b::ZPoly{T}) = ZP(_divrem(a.a, b.a)[1])
/{T}(a::ZPoly{T}, b::ZPoly{T}) = ZP(_divrem(a.a, b.a)[1])
rem{T}(a::ZPoly{T}, b::ZPoly{T}) = ZP(_divrem(a.a, b.a)[2])
mod{T}(a::ZPoly{T}, b::ZPoly{T}) = ZP(_divrem(a.a, b.a)[2])
%{T}(a::ZPoly{T}, b::ZPoly{T}) = ZP(_divrem(a.a, b.a)[2])
# ^ will come from power_by_squaring in stdlib
degree{T}(a::ZPoly{T}) = length(a) - 1

# defining bitwise ops seems odd, but makes sense for GF2, especially
# when using the compact binary repn below.  for N>2 we apply the
# operation to each coefficient pair in turn.

function same_length{T}(a::Vector{T}, b::Vector{T})
    big, small = length(a) > length(b) ? (a, b) : (b, a)
    pad = length(big) - length(small)
    small = vcat(small, zeros(T, pad))
    big, small
end

map{P<:ZPoly}(f::Function, a::P, b::P) = ZP(map(f, same_length(a.a, b.a)...))
for op in (:&, :|, :$)
    @eval $op{P<:ZPoly}(a::P, b::P) = map($op, a, b)
end
<<{T}(a::ZPoly{T}, n::Int) = ZPoly{T}(vcat(zeros(T, n), a.a))
>>>{T}(a::ZPoly{T}, n::Int) = ZPoly{T}(a.a[n+1:end])
>>{T}(a::ZPoly{T}, n::Int) = a >>> n

# --- GF(2) polynomials encoded as bits

# the type of Poly{} is the type returned by []
immutable GF2Poly{U<:Unsigned} <: Poly{ZField{Int,2,Int}}
    i::U
end

# constructors
GF2P{U<:Unsigned}(p::U) = GF2Poly{U}(p)
GF2P{S<:Signed}(p::S) = GF2P(unsigned(p))
GF2X{U<:Unsigned}(::Type{U}) = GF2Poly{U}(one(U) + one(U))
GF2X() = GF2X(UInt)

modulus(::Type{GF2Poly}) = 2
modulus(::GF2Poly) = 2

convert{I<:Integer, U<:Unsigned}(::Type{GF2Poly{U}}, i::I) = GF2P(convert(U, i % 2))
promote_rule{I<:Unsigned, U<:Unsigned}(::Type{GF2Poly{U}}, ::Type{I}) = GF2Poly{U}
promote_rule{U<:Unsigned}(::Type{GF2Poly{U}}, ::Type{Int}) = GF2Poly{U}

# promote to ZPoly mainly for equality tests (convert below)
promote_rule{U<:Unsigned,J<:Integer,I<:Integer}(::Type{GF2Poly{U}}, ::Type{ZPoly{ZField{J,2,I}}}) = ZPoly{ZField{J,2,I}}

# no setindex! as immutable
getindex{U<:Unsigned}(a::GF2Poly{U}, i) = a.i & (one(U) << (i-1)) == 0 ? GF2(0) : GF2(1)
length{U<:Unsigned}(a::GF2Poly{U}) = 8 * sizeof(U) - leading_zeros(a.i)

# for display we first convert to to ZPoly.  this is not efficient,
# but we probably need the interop anyway, and who cares if printing
# is efficient?

convert{J<:Integer, I<:Integer, U<:Unsigned}(::Type{ZPoly{ZField{J,2,I}}}, p::GF2Poly{U}) = _convert(ZPoly{ZField{J,2,I}}, p)
convert{J<:Integer, I<:Integer, U<:Unsigned}(::Type{GF2Poly{U}}, p::ZPoly{ZField{J,2,I}}) = _convert(GF2Poly{U}, p)
convert{U<:Unsigned}(::Type{GF2Poly{U}}, p::GF2Poly{U}) = p

function _convert{T, F}(::Type{T}, from::F)
    result, mask = zero(T), one(T)
    for _ in 1:length(from)
        if from & one(F) != zero(F)
            result += mask
        end
        from >>>= 1
        mask <<= 1
    end
    result
end

#print{U<:Unsigned}(io::IO, p::GF2Poly{U}) = print(io, convert(ZPoly{ZField{2,Int}}, p))
#print_no_mod{U<:Unsigned}(io::IO, p::GF2Poly{U}) = print_no_mod(io, convert(ZPoly{ZField{2,Int}}, p))
show_mod{U<:Unsigned}(io::IO, p::GF2Poly{U}) = print(io, " mod 2")
show{U<:Unsigned}(io::IO, p::GF2Poly{U}) = print(io, "GF2Poly{$U}($(p.i))")

bits(p::GF2Poly) = bits(p.i)

zero{U<:Unsigned}(::Type{GF2Poly{U}}) = GF2P(zero(U))
zero{U<:Unsigned}(::GF2Poly{U}) = GF2P(zero(U))
one{U<:Unsigned}(::Type{GF2Poly{U}}) = GF2P(one(U))
one{U<:Unsigned}(::GF2Poly{U}) = GF2P(one(U))

-{U<:Unsigned}(a::GF2Poly{U}) = a
for (name, op) in ((:+, :$), (:-, :$), (:|, :|), (:&, :&), (:$, :$))
    @eval $name{U<:Unsigned}(a::GF2Poly{U}, b::GF2Poly{U}) = GF2Poly{U}($op(a.i, b.i))
end
for op in (:(==), :<=, :<)
    @eval $op{U<:Unsigned}(a::GF2Poly{U}, b::GF2Poly{U}) = $op(a.i, b.i)
end
<<{U<:Unsigned}(a::GF2Poly{U}, n::Int) = GF2Poly(a.i << n)
>>>{U<:Unsigned}(a::GF2Poly{U}, n::Int) = GF2Poly(a.i >>> n)
>>{U<:Unsigned}(a::GF2Poly{U}, n::Int) = a >>> n
leading_zeros{U<:Unsigned}(a::GF2Poly{U}) = leading_zeros(a.i)

function *{U<:Unsigned}(a::GF2Poly{U}, b::GF2Poly{U})
    big, small = a > b ? (a, b) : (b, a)
    if small == zero(GF2Poly{U})  # no need to test big
        small
    elseif small == one(GF2Poly{U})
        big
    else
        result = zero(GF2Poly{U})
        while small != zero(GF2Poly{U}) && big != zero(GF2Poly{U})
            if small & one(GF2Poly{U}) == one(GF2Poly{U})
                result += big
            end
            small >>>= 1
            big <<= 1
        end
        result
    end
end

function divrem{U<:Unsigned}(a::GF2Poly{U}, b::GF2Poly{U})
    ONE, ZERO = one(GF2Poly{U}), zero(GF2Poly{U})
    b != ZERO || error("division by zero polynomial")
    if a == ZERO
        (a, a)
    elseif b == ONE
        (a, ZERO)
    elseif a == b
        (ONE, ZERO)
    else
        shift = leading_zeros(b) - leading_zeros(a)
        if shift < 0
            (ZERO, a)
        else
            rem::GF2Poly{U}, div::GF2Poly{U} = a, ZERO
            factor = ONE << shift
            b = b << shift  # no overflow (b *= factor)
            mask = ONE << (8 * sizeof(U) - leading_zeros(rem) - 1)  # msb(rem)
            for _ in shift:-1:0
                if rem & mask != ZERO
                    div += factor
                    rem -= b
                end
                factor >>>= 1
                b >>>= 1
                mask >>>= 1
            end
            div, rem
        end
    end
end

div{U<:Unsigned}(a::GF2Poly{U}, b::GF2Poly{U}) = divrem(a, b)[1]
/{U<:Unsigned}(a::GF2Poly{U}, b::GF2Poly{U}) = divrem(a, b)[1]
rem{U<:Unsigned}(a::GF2Poly{U}, b::GF2Poly{U}) = divrem(a, b)[2]
mod{U<:Unsigned}(a::GF2Poly{U}, b::GF2Poly{U}) = divrem(a, b)[2]
%{U<:Unsigned}(a::GF2Poly{U}, b::GF2Poly{U}) = divrem(a, b)[2]

degree{U<:Unsigned}(p::GF2Poly{U}) = 8*sizeof(U) - leading_zeros(p)

end

#ZField{ATP ZP(IntModN.ZField{2,Int64},1,0,0,0,1,1,0,1,1),
#       IntModN.ZPoly{IntModN.ZField{2,Int64}}}(x^6 + x^5 + x^3 + x^2 + x + 1 mod 2)
#ZField{ATP ZP(IntModN.ZField{2,Int64},1,0,0,0,1,1,0,1,1),
#       IntModN.ZPoly{IntModN.ZField{2,Int64}}}(x^6 + x^4 + x + 1 mod 2)


text:
@testset "FourierNeuralOperator" begin
    m = FourierNeuralOperator()

    𝐱, 𝐲 = rand(Float32, 2, 1024, 5), rand(Float32, 1024, 5)
    @test size(m(𝐱)) == size(𝐲)

    loss(𝐱, 𝐲) = sum(abs2, 𝐲 .- m(𝐱)) / size(𝐱)[end]
    data = [(𝐱, 𝐲)]
    Flux.train!(loss, params(m), data, Flux.ADAM())
end

@testset "MarkovNeuralOperator" begin
    m = MarkovNeuralOperator()

    𝐱, 𝐲 = rand(Float32, 1, 64, 64, 5), rand(Float32, 1, 64, 64, 5)
    @test size(𝐱) == size(m(𝐱)) == size(𝐲)

    loss(𝐱, 𝐲) = sum(abs2, 𝐲 .- m(𝐱)) / size(𝐱)[end]
    data = [(𝐱, 𝐲)]
    Flux.train!(loss, params(m), data, Flux.ADAM())
end


text:
"Counts difference between to strings"
function distance(a, b)
    length(a) != length(b) && throw(ArgumentError("Different lengths"))
    sum(map(!=, a, b))
end

max_stars_repo_name: koladilip/julia-exercism

max_issues_repo_name: koladilip/julia-exercism

max_forks_repo_name: koladilip/julia-exercism


text:
if length(ARGS) != 0
    N = parse(Int,ARGS[1])
else
    N = 1
end

@show typeof(N)

println("N = $N")


text:
using Random; import Future
using Base
using Plots

include("rendering.jl")

function run_agent!(env::TE, agent::TA, num_episodes::Int, rng::AbstractRNG)::Array{Float64, 1} where {TE<:AbstractEnvironment,TA<:AbstractAgent}
    returns = zeros(num_episodes)
    step::Int = 0
    for episode in 1:num_episodes
        env_reset!(env, rng)
        new_episode!(agent, rng)
        step = 0

        while !is_terminal(env)
            act!(agent, env, rng)
            returns[episode] += env.reward
        end
    end
    return returns
end

function eval_agent(env::TE, agent::TA, num_episodes::Int, rng::AbstractRNG)::Array{Float64, 1} where {TE<:AbstractEnvironment,TA<:AbstractAgent}
    returns = zeros(num_episodes)
    step::Int = 0
    for episode in 1:num_episodes
        env_reset!(env, rng)
        new_episode!(agent, rng)
        step = 0

        while !is_terminal(env)
            act!(agent, env, rng, false)  # set training flag to false
            returns[episode] += env.reward
        end
    end
    return returns
end

function eval_policy(env::TE, policy, γ, num_episodes::Int, rng::AbstractRNG)::Array{Float64, 1} where {TE<:AbstractEnvironment}
    returns = zeros(num_episodes)
    step::Int = 0
    for episode in 1:num_episodes
        env_reset!(env, rng)
        # new_episode!(agent, rng)
        step = 0

        while !is_terminal(env)
            get_action!(policy, env.state, rng)
            r = step!(env, policy.action, rng)
            returns[episode] += γ^step * r
            step += 1
        end
    end
    return returns
end

function run_agent!(env::AbstractEnvironment, agent::AbstractAgent, num_episodes::Int, rng::AbstractRNG, draw::Bool)::Array{Float64, 1}
    returns = zeros(num_episodes)
    step::Int = 0
    for episode in 1:num_episodes
        env_reset!(env, rng)
        new_episode!(agent, rng)
        step = 0
        if draw
            plt = plot(env)
            display(plt);
            gui()
        end

        while !is_terminal(env)
            act!(agent, env, rng)
            returns[episode] += env.reward
            if draw
                plt = plot(env)
                display(plt);
                gui()
            end
        end
    end
    return returns
end

function run_agent(env::AbstractEnvironment, agent::AbstractAgent, num_trials::Int, num_episodes::Int, rng::MersenneTwister)::Array{Float64, 2}
    returns = zeros(num_episodes, num_trials)
    agents = [clone(agent) for i in 1:num_trials]
    envs = [clone(env) for i in 1:num_trials]

    for trial in 1:num_trials
        for episode in 1:num_episodes
            env_reset!(envs[trial], rng)
            new_episode!(agents[trial], rng)

            while !is_terminal(envs[trial])
                act!(agents[trial], envs[trial], rng)
                returns[episode, trial] += envs[trial].reward
            end
        end
    end

    return returns
end

function run_agent_parallel(env::AbstractEnvironment, agent::AbstractAgent, num_trials::Int, num_episodes::Int, rng::MersenneTwister)::Array{Float64, 2}
    agents = [clone(agent) for i in 1:num_trials]
    envs = [clone(env) for i in 1:num_trials]

    rngs = [rng; accumulate(Future.randjump, fill(big(10)^20, num_trials-1), init=rng)]

    returns = zeros(num_episodes, num_trials)

    Threads.@threads for trial in 1:num_trials
        for episode in 1:num_episodes
            env_reset!(envs[trial], rngs[trial])
            new_episode!(agents[trial], rngs[trial])

            while !is_terminal(envs[trial])
                act!(agents[trial], envs[trial], rngs[trial])
                returns[episode, trial] += envs[trial].reward
            end
        end
    end

    return returns
end

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
using QuantumOptics: NLevelBasis, nlevelstate
using Test, IonSim
const pc = IonSim.PhysicalConstants
using Suppressor

@suppress_err begin

@testset "ions -- Ca40" begin
    C = Ca40()
    # test for required fields
    @test mass(C) == pc.m_ca40
    @test typeof(level_structure(C)) == OrderedDict{String,NamedTuple}
    @test selected_level_structure(C) == level_structure(C)
    @test typeof(matrix_elements(C)) == OrderedDict{Tuple,Function}
    @test selected_matrix_elements(C) == matrix_elements(C)
    @test ismissing(ion_number(C))
    @test ismissing(ion_position(C))
    @test typeof(stark_shift(C)) == OrderedDict{String,Real}

    # test ==
    C1 = Ca40()
    @test C1 == C

    # test inner constructors for Ca40
    C1 = Ca40(["S-1/2", "D-1/2"])
    @test collect(keys(C1.selected_level_structure)) == ["S-1/2", "D-1/2"]
    C1copy = copy(C1)
    @test selected_level_structure(C1copy) == selected_level_structure(C1)

    # make sure print/show don't throw any errors
    @suppress print(C); show(C)

    # test matrix_element
    E = 100randn(); γ = 100randn(); ϕ = 100randn()
    @test matrix_element(C, ["S-1/2", "D-1/2"], E, γ, ϕ) == C.matrix_elements[("S-1/2", "D-1/2")](E, γ, ϕ)

    # Test matrix_element function for a specific collecion of parameters that have been
    # pre-evaluated independently. Hopefully this will fail if the function is constructed
    # improperly
    @test C.selected_matrix_elements[("S-1/2", "D-1/2")](1, 0, 45) ≈ 4.2248 rtol=1e-4
    @test sum(C.selected_matrix_elements[("S-1/2", "D-1/2")].(0, 90, 0:0.1:90)) == 0

    # Test zero_stark_shift
    C.stark_shift["S-1/2"] = 10
    zero_stark_shift(C)
    @test sum(values(C.stark_shift)) == 0

    # make sure improper indexing of Ca40 yields an AssertionError
    @test_throws AssertionError C[""]

    # test indexing
    @test C1["S-1/2"].data == ComplexF64[1; 0]
    @test C1["D-1/2"].data == ComplexF64[0; 1]

    # test set properties
    @test_throws AssertionError C.selected_level_structure = []
    C.mass = 0  # shouldn't be allowed
    @test C.mass == pc.m_ca40
    C.selected_level_structure = ["S-1/2", "D-1/2"]
    @test C.selected_level_structure == C1.selected_level_structure

    # test get properties
    warning = "ion has not been added to a configuration"
    @test_logs (:warn, warning) C.number
    @test_logs (:warn, warning) C.position
end

@testset "ions -- general" begin

    # test specific case for Lande g-factor
    @test gJ(2, 5/2) == 6/5

    # test specific pre-evaluated case for zeeman_shift
    val = 4.198873488465285e6
    @test zeeman_shift(1e-4, 2, 5/2, 5/2) ≈ val
    @test zeeman_shift(1e-4, (2, 5/2, 5/2)) ≈ val
    @test zeeman_shift(B=1e-4, l=2, j=5/2, mⱼ=5/2) ≈ val
end
end  # end suppress


text:
using Documenter, BlockNLPModels
makedocs(sitename = "My Documentation")


text:
function test_memory_data_layer(backend::Backend, T, eps)
  println("-- Testing Memory Data Layer on $(typeof(backend)){$T}...")

  ############################################################
  # Prepare Data for Testing
  ############################################################
  batch_size = 6
  tensor_dim = abs(rand(Int)) % 6 + 1
  data_dim = tuple(rand(1:6, tensor_dim)...)
  println("    > $data_dim")

  data = rand(T, data_dim..., 9)
  mean_data = rand(T, data_dim..., 1)
  mean_blob = make_blob(backend, mean_data)

  ############################################################
  # Setup
  ############################################################

  # batch size is determined by
  layer = MemoryDataLayer(data = Array[data], tops = [:data], batch_size=batch_size,
      transformers=[(:data, DataTransformers.SubMean(mean_blob=mean_blob))])
  state = setup(backend, layer, Blob[], Blob[])

  data_idx = map(x->1:x, data_dim)
  layer_data = Array{eltype(data)}(tuple(data_dim..., batch_size))

  data_aug = cat(tensor_dim+1, data, data)
  data_aug .-= mean_data
  forward(backend, state, Blob[])
  copy!(layer_data, state.blobs[1])
  @test all(-eps .< layer_data - data_aug[data_idx..., 1:1+batch_size-1] .< eps)
  @test state.epoch == 0

  forward(backend, state, Blob[])
  copy!(layer_data, state.blobs[1])
  @test all(-eps .< layer_data - data_aug[data_idx..., batch_size+1:2batch_size] .< eps)
  @test state.epoch == 1

  forward(backend, state, Blob[])
  copy!(layer_data, state.blobs[1])
  @test all(-eps .< layer_data - data_aug[data_idx..., 2batch_size+1:3batch_size] .< eps)
  @test state.epoch == 2

  shutdown(backend, state)
end
function test_memory_data_layer(backend::Backend)
  test_memory_data_layer(backend, Float32, 1e-5)
  test_memory_data_layer(backend, Float64, 1e-10)
end

if test_cpu
  test_memory_data_layer(backend_cpu)
end
if test_gpu
  test_memory_data_layer(backend_gpu)
end


text:
export train

function train(params, opt::Adam; maxiter=2000)
    # mask used to fix first two elements
    msk = [false, false, true, true, true, true, true, true, true, true]
    pp = params[:,msk]
    for i=1:maxiter
        g = NiLang.AD.gradient(Val(1), embedding_loss, (0.0, params))[2][:,msk]
        update!(pp, g, opt)
        view(params, :, msk) .= pp
        if i%1000 == 0
            println("Step $i, loss = $(embedding_loss(0.0, params)[1])")
        end
    end
    params
end

using Optim
function train(params, opt::NewtonTrustRegion)
    # mask used to fix first two elements
    msk = [false, false, true, true, true, true, true, true, true, true]
    i = Ref(1)
    function f(x)
        vec(view(params,:,msk)) .= x
        l = embedding_loss(0.0, params)[1]
        println("Step $(i[]), loss = $l")
        return l
    end
    function g!(G, x)
        vec(view(params,:,msk)) .= x
        G .= vec(NiLang.AD.gradient(Val(1), embedding_loss, (0.0, params))[2][:,msk])
    end
    function h!(H, x)
        vec(view(params,:,msk)) .= x
        nm = sum(msk)*size(params, 1)
        i[] += 1
        H .= reshape(reshape(get_hessian(params), size(params)..., size(params)...)[:,msk,:,msk], nm, nm)
    end
    optimize(f, g!, h!, vec(params[:,msk]), opt)
    params
end


text:
using Soss

m = @model begin
    x ~ Normal() |> iid(2)
    xxt = x * x'
end;

s = rand(m(), 1000)


text:
using Coverage
import DataFrames
import DelimitedFiles

thisFolder = Pkg.dir("DecisionTrees")
thisFolder = joinpath(thisFolder, "..", "memHist", "2018 06 24c")
isdir(thisFolder)
mr = analyze_malloc(thisFolder)
mr2 = map(x->[x.bytes,x.filename,x.linenumber], mr);
mr3 = deepcopy(mr2[1])
for i = 2:length(mr2)
	mr3 = hcat(mr3, mr2[i])
end
mr4 = permutedims(mr3, [2,1]);
# mr4=readcsv("C:\\temp\\alloc.csv")

function list_of_files(;fldr=pwd(),ext=".jl")
    li = Vector{Vector{String}}()
    for (root, dirs, files) in walkdir(fldr)
        # println("Directories in $root")
        # for dir in dirs
        #   println(joinpath(root, dir)) # path to directories
        # end
        # println("Files in $root")
        for file in files
            thisf = joinpath(root, file)
            if !occursin(".git", thisf) && (lowercase(thisf[end - length(ext) + 1:end]) == lowercase(ext))
                # println(thisf) # path to files
                fn_only = splitdir(thisf)[2]
                wo_fldr = replace(thisf, fldr => "")
                if wo_fldr[1] == '\\' || wo_fldr[1] == '/'
                    wo_fldr = wo_fldr[2:end]
                end
                push!(li, [thisf,wo_fldr,fn_only])
            end
        end
    end
    return li
end

# add new column
newcolumn = similar(mr4[:,1])
newcolumn2 = similar(mr4[:,1])
mr5 = hcat(mr4, newcolumn, newcolumn2)
mr5[:,4] .= "file not found"
mr5[:,5] .= "code not found"

# lookup filenames of *.mem
filelist = unique(mr4[:,2])
filelist2 = deepcopy(filelist)
for i = 1:length(filelist)
    fi = filelist[i][1:end - 4]
    fi2 = splitdir(fi)[2]
    filelist2[i] = fi2
    # @show isfile(fi)
end

# generate list of all *.jl files
file_list_of_julia_code = list_of_files(fldr=Pkg.dir("DecisionTrees"))
file_list_of_julia_code2 = map(x->x[3], file_list_of_julia_code)

filelocations = Dict{String,String}()
for ii = 1:size(filelist2, 1)
    fi = filelist2[ii]
    fiMem = filelist[ii]
    idx = indexin([fi], file_list_of_julia_code2)
    if idx[1] == nothing
        @warn("Oops: this one was not found")
        @show fi
    else
        # @show file_list_of_julia_code[idx[1]][1]
        # @show fi
        filelocations[fiMem] = file_list_of_julia_code[idx[1]][1]
    end
end

# add location
for i = 1:size(mr5, 1)
    if haskey(filelocations, mr5[i,2])
        mr5[i,4] = filelocations[mr5[i,2]]
    end
end

codeDict = Dict{String,Vector{String}}()
for fi in values(filelocations)
    if isfile(fi)
        codeDict[fi] = readlines(fi)
    end
end

for memline = 1:size(mr5, 1)
    # memline=2118
    fi = mr5[memline,4]
    lineNr = mr5[memline,3]
    # read the file
    if haskey(codeDict, fi)
        rl = codeDict[fi]
        if length(rl) >= lineNr
            codeOfThatLine = rl[lineNr]
            mr5[memline,5] = codeOfThatLine
        end
    end
end

fnOnly = map(x->x[2], splitdir.(mr5[:,4]))
mr6 = hcat(mr5, fnOnly)

# sort data
srtp = sortperm(convert(Vector{Int}, mr6[:,1]), rev=true)

# define header
hdr = ["Allocations" "MemFile" "LineNumber" "CodeFile" "CodeInThatLine" "FileName Only"]
mr7 = vcat(hdr, mr6)

# change columns
tmp = mr7[:,4]
mr7[:,4] = mr7[:,3]
mr7[:,3] = tmp

DelimitedFiles.writedlm("C:\\temp\\allocations.csv",mr7,',');

#=
include(joinpath(Pkg.dir("DecisionTrees"),"testing","allocation_coverage_analysis.jl")) =#


text:
# Script to sample the sensitivity of the C3a residual for the Alternate pathway -

# Run the model for the ensemble of solutions -
include("SolveBalances.jl")

# some global parameters -
BIG = 1e10
SMALL = 1e-6

using PyCall
@pyimport numpy as np

function calculate_error_obj(T,X,MEASURED_ARRAY_C3a,MEASURED_ARRAY_C5a)

  obj_array = BIG*ones(2)

  # Need to interpolate the simulation onto the experimental time scale -
  number_of_measurements = 20
  IC3a = np.interp(MEASURED_ARRAY_C3a[:,1],T,X[:,19])
  IC5a = np.interp(MEASURED_ARRAY_C5a[:,1],T,X[:,15])

  # Compute the error values -
  MC3a = 0.11002*MEASURED_ARRAY_C3a[:,2]
  MC5a = 0.00009615*MEASURED_ARRAY_C5a[:,2]
  # error_C3a = sum((MC3a - IC3a).^2)
  # error_C5a = sum((MC5a - IC5a).^2)

  # Compute the scaled IC3a and IC5a -
  min_IC3a = minimum(IC3a)
  max_IC3a = maximum(IC3a)
  min_IC5a = minimum(IC5a)
  max_IC5a = maximum(IC5a)

  min_MC3a = minimum(MC3a)
  max_MC3a = maximum(MC3a)
  min_MC5a = minimum(MC5a)
  max_MC5a = maximum(MC5a)

  delta_C3a = max_IC3a - min_IC3a
  delta_C5a = max_IC5a - min_IC5a

  # if (delta_C3a<SMALL)
  #   delta_C3a = SMALL
  # end
  #
  # if (delta_C5a<SMALL)
  #   delta_C5a = SMALL
  # end

  shape_C3a = (IC3a - min_IC3a)/(delta_C3a)
  shape_C5a = (IC5a - min_IC5a)/(delta_C5a)

  shape_MC3a = (MC3a - min_MC3a)/(max_MC3a - min_MC3a)
  shape_MC5a = (MC5a - min_MC5a)/(max_MC5a - min_MC5a)

  # Compute the scale error -
  scale_error_C3a = (max_IC3a - max_MC3a)/(max_MC3a)
  scale_error_C5a = (max_IC5a - max_MC5a)/(max_MC5a)

  #error_C3a = sum((shape_C3a - shape_MC3a).^2)+(scale_error_C3a)^2
  #error_C5a = sum((shape_C5a - shape_MC5a).^2)+(scale_error_C5a)^2
  error_C3a = sum((IC3a - MC3a).^2)
  error_C5a = sum((IC5a - MC5a).^2)

  # Compute the objective array -
  obj_array[1] = error_C3a
  obj_array[2] = error_C5a

  value = sum(obj_array);
  if (value>10)
    value = 10
  end

  return (value)
end

function main()

  # Load the data -
  MEASURED_ARRAY_C3a = readdlm("./data/Shaw2015_Fig2e_C3a.txt",',')
  MEASURED_ARRAY_C5a = readdlm("./data/Shaw2015_Fig3c_C5a.txt",',')

  # Load the paremeter samples -
  pc_array_full = readdlm("./sens/Samples.txt")
  (nrows,ncols) = size(pc_array_full)

  # for each sample, grab the parameters, solve the model and compute the error -
  data_dictionary = Dict()
  initial_condition_array = zeros(19)
  initial_condition_array[1] = 1.0      # 1 zymosan
  initial_condition_array[2] = 1.9      # 2 C4
  initial_condition_array[3] = 0.322    # 3 C2
  initial_condition_array[8] = 7.57     # 4 C3
  initial_condition_array[14] = 0.195   # 5 C5
  initial_condition_array[17] = 2.23    # 6 Factor H
  initial_condition_array[18] = 0.417   # 7 C4BP
  data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;

  tStart = 0.0;
  tStop = 25.0;
  tStep = 0.1;

  # initialize the error_array -
  error_array = zeros(nrows)

  for sample_index = 1:nrows

    # Set the parameters -
    parameter_array = pc_array_full[sample_index,:]
    data_dictionary["PARAMETER_ARRAY"] = parameter_array

    # Run the model -
    (T,X) = SolveBalances(tStart,tStop,tStep,data_dictionary)

    # Calculate the error -
    error_array[sample_index] = calculate_error_obj(T,X,MEASURED_ARRAY_C3a,MEASURED_ARRAY_C5a)

    @show sample_index
  end

  return error_array
end


text:
module ViewHelper

using Genie, Genie.Helpers

export output_flash

function output_flash(params::Dict{Symbol,Any}) :: String
  ! isempty( flash(params) ) ? """
$(flash(params))
""" : ""
end

end

max_stars_repo_name: Fatalon/Dashboard

max_issues_repo_name: Fatalon/Dashboard

max_forks_repo_name: Fatalon/Dashboard


text:
include("panX.jl")

Export = Command(
   "export",
   "pangraph export <options> [arguments]",
   "exports a pangraph to a chosen file format(s)",
   """zero or one pangraph file (native json)
      if no file given, reads from stdin
      stream can be optionally gzipped.""",
   [
    Arg(
        Int,
        "minimum block length for edge connection",
        (short="-ell", long="--edge-minimum-length"),
        "blocks below this length cutoff will be ignored for edges in graph",
        200,
    ),
    Arg(
        Int,
        "maximum block length for edge connection",
        (short="-elu", long="--edge-maximum-length"),
        "blocks above this length cutoff will be ignored for edges in graph",
        typemax(Int),
    ),
    Arg(
        Int,
        "minimum block depth for edge connection",
        (short="-edl", long="--edge-mininum-depth"),
        "blocks below this depth cutoff will be ignored for edges in graph",
        0,
    ),
    Arg(
        Int,
        "maximum block depth for edge connection",
        (short="-edu", long="--edge-maximum-depth"),
        "blocks above this depth cutoff will be ignored for edges in graph",
        typemax(Int),
    ),
    Arg(
        Int,
        "minimum block length for export",
        (short="-ll", long="--minimum-length"),
        "blocks below this length cutoff will not be exported",
        0,
    ),
    Arg(
        Int,
        "maximum block length for export",
        (short="-lu", long="--maximum-length"),
        "blocks above this length cutoff will not be exported",
        typemax(Int),
    ),
    Arg(
        Int,
        "minimum block depth for export",
        (short="-dl", long="--mininum-depth"),
        "blocks below this depth cutoff will not be exported",
        0,
    ),
    Arg(
        Int,
        "maximum block depth for export",
        (short="-du", long="--maximum-depth"),
        "blocks above this depth cutoff will not be exported",
        typemax(Int),
    ),
    Arg(
        String,
        "output directory",
        (short="-o", long="--output-directory"),
        "relative path to directory where output will be stored",
        "export"
    ),
    Arg(
        String,
        "prefix",
        (short="-p", long="--prefix"),
        "basename of files",
        "pangraph"
    ),
    Arg(
        Bool,
        "no GFA export",
        (short="-ng", long="--no-export-gfa"),
        "do not emit GFA file",
        false,
    ),
    Arg(
        Bool,
        "export panX visualization",
        (short="-pX", long="--export-panX"),
        "emit vis directory to input to panX-visualization",
        false,
    ),
    Arg(
        Bool,
        "no duplications in export",
        (short="-nd", long="--no-duplications"),
        "do not export any block that contains at least one strain more than once",
        false
    ),
   ],

   function(args)
       path = parse(Export, args)
       path = if (path === nothing || length(path) == 0)
           nothing
       elseif length(path) == 1
           path
       else
           return 2
       end

       graph = load(path, Export)

       function cutoffs(prefix)
           (
              depth = (
                   min = arg(Export, "-$(prefix)dl"),
                   max = arg(Export, "-$(prefix)du"),
              ),
              length = (
                   min = arg(Export, "-$(prefix)ll"),
                   max = arg(Export, "-$(prefix)lu"),
              )
           )
       end
       noduplications = arg(Export, "-nd")

       node = cutoffs("")
       edge = cutoffs("e")

       filter = (
            connect = function(node)
                return (!(edge.depth.min  ≤ Blocks.depth(node.block)  ≤ edge.depth.max)
                     || !(edge.length.min ≤ Blocks.length(node.block) ≤ edge.length.max))
            end,
            output = function(segment, isolates)
                return (!(node.depth.min  ≤ segment.depth  ≤ node.depth.max)
                     || !(node.length.min ≤ length(segment.sequence) ≤ node.length.max)
                     ||  (noduplications && (length(isolates) != length(unique(isolates)))))
            end
       )

       # Create directory if it doesn't exist
       directory = arg(Export, "-o")
       if !isdir(directory)
           mkpath(directory)
       end
       prefix = arg(Export, "-p")

       # GFA export (default)
       if !arg(Export, "-ng")
           Base.open("$(directory)/$(prefix).gfa", "w") do io
               marshal(io, graph; fmt=:gfa, opt=filter)
           end
           Base.open("$(directory)/$(prefix).fa", "w") do io
               marshal(io, graph; fmt=:fa)
           end
       end

       # panX export (doesn't fit into marshal paradigm)
       if arg(Export, "-pX")
           if !Shell.havecommand("fasttree")
               panic("external command fasttree not found. please install before attempting to export as panX visualization\n")
           end

           vis = "$(directory)/vis"
           if !isdir(vis)
               mkpath(vis)
           end

           PanX.emit(graph, vis)
       end

        # TODO: others?
   end
)

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
using LinearAlgebraicRepresentation
using ViewerGL
using LinearAlgebra
using SparseArrays
using DataStructures
Lar = LinearAlgebraicRepresentation
GL = ViewerGL

function testarrangement(V,CV,FV,EV)
		cop_EV = Lar.coboundary_0(EV::Lar.Cells);
		cop_FE = Lar.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells);
		W = convert(Lar.Points, V');

		V, copEV, copFE, copCF = Lar.space_arrangement(
				W::Lar.Points, cop_EV::Lar.ChainOp, cop_FE::Lar.ChainOp);

		V = convert(Lar.Points, V');
		V,CVs,FVs,EVs = Lar.pols2tria(V, copEV, copFE, copCF) # whole assembly
		GL.VIEW(GL.GLExplode(V,FVs,1.1,1.1,1.1,99,1));
		GL.VIEW(GL.GLExplode(V,EVs,1.5,1.5,1.5,99,1));
		GL.VIEW(GL.GLExplode(V,CVs,1,1,1,99,1));
end

V = [0.6540617999999998 0.2054992 0.2972308000000001; 0.7142364999999998 0.1455625 0.969203; 0.6872816200343089 0.1724107693563844 0.6681972618528844; 0.5941250999999997 0.8769965 0.36249240000000005; 0.6542998079614951 0.8170598007732168 1.0344646993560203; 0.6929229746417847 0.3843473621810608 0.9924101013070119; 0.6660927570227951 0.3906268097799347 0.66824889023436; 0.6924256885846046 0.3894687221287142 0.9924113129852241; 1.3260341 0.2707609000000001 0.2428771; 1.3862088 0.21082410000000004 0.9148494000000003; 1.2660973 0.9422582000000003 0.3081388; 1.3262719920385087 0.8823213992267827 0.9801110006439793; 0.7946185131077457 0.1533691699791747 0.9627011677876783; 0.795311803386693 0.18207414025178234 0.6679835359595829; 0.6978893364451173 0.3892254686183435 0.9924004714668809; 0.8049461769791598 0.3844586322271989 0.9823884300888783; -0.38740630000000004 0.49022260000000006 0.45363390000000015; 0.3249123 0.707347 0.5231232000000001; -0.17028190000000007 -0.08642419999999995 0.029717700000000125; 0.542036694733776 0.13070013212634024 0.09920695360173326; -0.317917 0.06630639999999999 1.0658723; 0.3944015999999999 0.28343080000000004 1.1353616; -0.10079260000000007 -0.5103404 0.6419561000000001; 0.6115259947337761 -0.2932160678736597 0.7114453536017331; 0.4729938678063347 0.07470259548304928 0.9819171231357612; 0.4673971777131105 0.08956648867512645 0.992844167328483; 0.4666328572520818 0.07498584149882567 0.9828352213039846; 0.4707836778272794 0.18222015069807973 0.6686268661789595; 0.4673140236772549 0.08983923498806379 0.9928443441979011; 0.5039882186318663 0.07330882745780876 0.6685562398308738; 0.7899025999999999 0.060579300000000016 0.6679888999999999; 0.46601 0.0749997 0.6686315999999999; 0.8043229999999999 0.38447249999999994 0.6679745999999999; 0.4804303999999999 0.3988929 0.6686172999999999; 0.7905452000000002 0.060564999999999994 0.9922023000000001; 0.46665270000000003 0.0749854 0.9928450000000001; 0.8049656 0.3844582 0.9921880000000001; 0.4810731000000001 0.3988786 0.9928307000000001]
copEV = SparseArrays.sparse([1, 4, 14, 2, 5, 15, 1, 2, 7, 17, 3, 4, 19, 3, 6, 20, 5, 6, 8, 21, 7, 9, 46, 47, 8, 9, 54, 55, 10, 12, 14, 10, 13, 16, 11, 12, 19, 11, 13, 20, 15, 16, 18, 23, 17, 18, 48, 49, 21, 22, 53, 54, 22, 23, 62, 63, 24, 26, 36, 24, 27, 37, 25, 26, 38, 25, 27, 39, 28, 30, 36, 28, 31, 37, 29, 30, 38, 29, 33, 39, 32, 33, 35, 41, 31, 32, 34, 43, 65, 66, 67, 34, 35, 60, 61, 40, 42, 50, 51, 67, 40, 43, 57, 58, 41, 42, 44, 45, 44, 48, 59, 45, 50, 60, 46, 49, 62, 47, 51, 64, 52, 56, 59, 52, 57, 61, 65, 53, 56, 63, 55, 58, 64, 66], [1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 10, 10, 10, 11, 11, 11, 12, 12, 12, 13, 13, 13, 13, 14, 14, 14, 14, 15, 15, 15, 15, 16, 16, 16, 16, 17, 17, 17, 18, 18, 18, 19, 19, 19, 20, 20, 20, 21, 21, 21, 22, 22, 22, 23, 23, 23, 24, 24, 24, 25, 25, 25, 25, 26, 26, 26, 26, 26, 26, 26, 27, 27, 27, 27, 28, 28, 28, 28, 28, 29, 29, 29, 29, 30, 30, 30, 30, 31, 31, 31, 32, 32, 32, 33, 33, 33, 34, 34, 34, 35, 35, 35, 36, 36, 36, 36, 37, 37, 37, 38, 38, 38, 38], Int8[-1, -1, -1, -1, -1, -1, 1, 1, -1, -1, -1, 1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, -1, 1, 1, -1, -1, -1, -1, 1, 1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1, -1, -1, -1, -1, -1, 1, -1, -1, -1, 1, -1, 1, 1, -1, -1, -1, 1, 1, -1, 1, -1, 1, 1, 1, -1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, -1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, 1, -1, -1, 1, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])
copFE = SparseArrays.sparse([1, 4, 2, 5, 1, 6, 1, 7, 2, 8, 1, 9, 1, 2, 18, 20, 1, 2, 21, 22, 1, 2, 25, 26, 3, 4, 3, 6, 3, 7, 3, 9, 4, 7, 5, 8, 4, 9, 4, 5, 18, 20, 4, 5, 28, 29, 6, 7, 6, 9, 8, 9, 21, 22, 8, 9, 25, 27, 8, 9, 28, 29, 10, 13, 10, 14, 10, 15, 10, 16, 11, 13, 11, 14, 11, 15, 11, 16, 12, 17, 11, 16, 11, 12, 30, 32, 11, 12, 23, 24, 13, 15, 13, 16, 14, 15, 14, 16, 16, 17, 16, 17, 23, 24, 16, 17, 18, 19, 16, 17, 18, 23, 19, 24, 20, 25, 18, 26, 18, 28, 20, 29, 19, 30, 18, 31, 21, 23, 21, 27, 22, 25, 21, 26, 21, 28, 21, 21, 23, 28, 24, 30, 23, 32, 25, 29, 27, 28, 26, 31, 32, 31, 30, 31], [1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15, 16, 16, 17, 17, 17, 17, 18, 18, 18, 18, 19, 19, 20, 20, 21, 21, 21, 21, 22, 22, 22, 22, 23, 23, 23, 23, 24, 24, 25, 25, 26, 26, 27, 27, 28, 28, 29, 29, 30, 30, 31, 31, 32, 32, 33, 33, 34, 34, 34, 34, 35, 35, 35, 35, 36, 36, 37, 37, 38, 38, 39, 39, 40, 40, 41, 41, 41, 41, 42, 42, 42, 42, 43, 43, 44, 44, 45, 45, 46, 46, 47, 47, 48, 48, 49, 49, 50, 50, 51, 51, 52, 52, 53, 53, 54, 54, 55, 55, 56, 56, 57, 58, 59, 59, 60, 60, 61, 61, 62, 62, 63, 63, 64, 64, 65, 66, 67, 67], Int8[1, 1, -1, -1, -1, 1, -1, 1, 1, 1, -1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, -1, 1, -1, 1, -1, -1, -1, -1, -1, 1, -1, -1, -1, 1, -1, -1, 1, 1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, -1, 1, -1, 1, -1, -1, 1, 1, -1, 1, 1, -1, 1, -1, 1, 1, 1, 1, -1, -1, -1, -1, -1, 1, -1, -1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, 1, -1, -1, 1, 1, 1, 1, -1, 1, -1, -1, 1, -1, 1, -1, 1, -1, 1, -1, -1, 1, 1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1, -1, 1, 1, 1, 1, -1, -1, -1, -1, -1, 1, -1, -1, -1, 1, 1, -1, 1, -1, -1, 1, 1, 1, 1, -1, -1, 1])
size(V), size(copEV), size(copFE)

copFV = copFE * copEV
FV = [copFV[f,:].nzind for f=1:size(copFV,1)]
EV = [copEV[e,:].nzind for e=1:size(copEV,1)]
FE = Lar.cop2lar(copFE)
ff = copFV * copFV'
FF = [ff[f,:].nzind for f = 1:size(ff, 2)]
VV = [[v] for v=1:size(V,1)]
Model = (Lar.Points(V'), Lar.Cells[VV,EV,FV])

for k=1:length(FV)
	fVV = [[v] for v in Set(cat([FV[f] for f in FF[k]]))]
	fE = Int64[e for e in  Set(cat([FE[f] for f in FF[k]]))]
	fEV = sort!([EV[e] for e in fE])
	fFV = vcat([FV[f] for f in FF[k]])
	model = (Lar.Points(V'), Lar.Cells[fVV,fEV,fFV]);
	GL.VIEW(GL.numbering(.1)(model, GL.COLORS[1], 0.1));
end


text:
using Flux, Test, CUDA

@info "Testing GPU Support"
CUDA.allowscalar(false)

include("bayesian.jl")
include("layers/mclayers.jl")
include("layers/batchensemble.jl")
include("layers/variational.jl")
include("layers/bayesianBE.jl")
include("layers/spectrallayers.jl")