text:
module Test_nullrange

using DescriptorSystems
using MatrixPencils
using LinearAlgebra
using Polynomials
using Test

println("Test_nullrange")
@testset "nullrange" begin

@testset "grange" begin

sys = rdss(0,0,0);
@time sysr, sysx, info = grange(sys, inner = true)
@test gnrank(sys-sysr*sysx) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysr'*sysr-I) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysr) & (isproper(sys) ? isproper(sysx) : true) &&
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      isempty(gzero(gminreal(sysx))) &&
      info.nrank == 0 && ismissing(info.nfuz) && info.niuz == 0

@time sysr, sysx, info = grange(sys)
@test gnrank(sys-sysr*sysx) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      isproper(sysr) & (isproper(sys) ? isproper(sysx) : true) &&
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      isempty(gzero(gminreal(sysx))) &&
      isempty(gzero(gminreal(sysr))) &&
      info.nrank == 0 && ismissing(info.nfuz) && info.niuz == 0

sys = rdss(0,0,0,disc=true);
@time sysr, sysx, info = grange(sys, inner = true)
@test gnrank(sys-sysr*sysx) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysr'*sysr-I) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysr) & (isproper(sys) ? isproper(sysx) : true) &&
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      isempty(gzero(gminreal(sysx))) &&
      info.nrank == 0 && ismissing(info.nfuz) && info.niuz == 0

@time sysr, sysx, info = grange(sys)
@test gnrank(sys-sysr*sysx) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      isproper(sysr) & (isproper(sys) ? isproper(sysx) : true) &&
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      isempty(gzero(gminreal(sysx))) &&
      isempty(gzero(gminreal(sysr))) &&
      info.nrank == 0 && ismissing(info.nfuz) && info.niuz == 0

a = -(1. +2im)*ones(1,1); e = (1. -im)*ones(1,1); b = (2. +im)*ones(1,1);
c = (-1. +im)*ones(1,1); d = complex(0)*ones(1,1);
sys = dss(a,e,b,c,d);
@time sysr, sysx, info = grange(sys, inner = true, atol1 = 1.e-7, atol2 = 1.e-7)
@test gnrank(sys-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysr'*sysr-I) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysr) & isproper(sysx) && # checking properness of factors
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      count(t -> isinf(t), gzero(gminreal(sysx))) == 1 &&
      info.nrank == 1 && ismissing(info.nfuz) && info.niuz == 1

@time sysr, sysx, info = grange(sys, zeros = "unstable", inner = true, atol1 = 1.e-7, atol2 = 1.e-7)
@test gnrank(sys-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysr'*sysr-I) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysr) & isproper(sysx) && # checking properness of factors
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      count(t -> isinf(t), gzero(gminreal(sysx))) == 1 &&
      info.nrank == 1 && info.nfuz == 0 && info.niuz == 1

a = -(1. +2im)*ones(1,1); e = (1. -im)*ones(1,1); b = (2. +im)*ones(1,1);
c = (-1. +im)*ones(1,1); d = complex(0)*ones(1,1);
sys = dss(a,e,b,c,d,Ts = 1);
@time sysr, sysx, info = grange(sys, inner = true, atol1 = 1.e-7, atol2 = 1.e-7)
@test gnrank(sys-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysr'*sysr-I) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysr) & isproper(sysx) && # checking properness of factors
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      count(t -> isinf(t), gzero(gminreal(sysx))) == 1 &&
      info.nrank == 1 && ismissing(info.nfuz) && info.niuz == 0

@time sysr, sysx, info = grange(sys, zeros = "unstable", inner = true, atol1 = 1.e-7, atol2 = 1.e-7)
zer = gzero(gminreal(sysx));
@test gnrank(sys-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysr'*sysr-I,atol=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysr) & isproper(sysx) && # checking properness of factors
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      count(t -> abs.(t) .>= 1, zer) == 0 &&
      info.nrank == 1 && info.nfuz == 0 && info.niuz == 0

sys = rss(0,1,0);
@time sysr, sysx, info = grange(sys, inner = true)
r = size(sysr,2);
@test gnrank(sys-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysr'*sysr-I,atol=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysr) & isproper(sysx) && # checking properness of factors
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      isempty(gzero(gminreal(sysx))) &&
      info.nrank == r && ismissing(info.nfuz) && info.niuz == 0

sys = rss(0,0,2);
@time sysr, sysx, info = grange(sys, inner = true)
r = size(sysr,2);
@test gnrank(sys-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysr'*sysr-I,atol=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysr) & isproper(sysx) && # checking properness of factors
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      isempty(gzero(gminreal(sysx))) &&
      info.nrank == r && ismissing(info.nfuz) && info.niuz == 0

sys = dss([1 0]);
@time sysr, sysx, info = grange(sys, inner = true)
r = size(sysr,2);
@test gnrank(sys-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysr'*sysr-I,atol=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysr) & isproper(sysx) && # checking properness of factors
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      isempty(gzero(gminreal(sysx))) &&
      info.nrank == r && ismissing(info.nfuz) && info.niuz == 0

sys = dss([0; 1]);
@time sysr, sysx, info = grange(sys, inner = true)
r = size(sysr,2);
@test gnrank(sys-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysr'*sysr-I,atol=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysr) & isproper(sysx) && # checking properness of factors
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      isempty(gzero(gminreal(sysx))) &&
      info.nrank == r && ismissing(info.nfuz) && info.niuz == 0

sys = rss(0,3,2);
@time sysr, sysx, info = grange(sys, inner = true)
r = size(sysr,2);
@test gnrank(sys-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysr'*sysr-I,atol=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysr) & isproper(sysx) && # checking properness of factors
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      isempty(gzero(gminreal(sysx))) &&
      info.nrank == r && ismissing(info.nfuz) && info.niuz == 0

# large compilation times
@time z = rtf('z');     # define the variable z as rational function
@time Gd = [z^2 z/(z-2); 0 1/z];     # define the 2-by-2 improper Gd(z)
@time sysd = dss(Gd,minimal = true,Ts = 1);
zeref = gzero(sysd,atol1=1.e-7)

@time sysr, sysx, info = grange(sysd, zeros = "none", inner = false, atol = 1.e-7)
zer = gzero(sysr,atol1=1.e-7)
r = size(sysr,2);
@test gnrank(sysd-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      isproper(sysr) && # checking properness of factors
      length(zer) == 0  && # checking lack of zeros
      info.nrank == r && ismissing(info.nfuz) && info.niuz == 0

@time sysr, sysx, info = gcrange(sysd, zeros = "none", coinner = false, atol = 1.e-7)
zer = gzero(sysr,atol1=1.e-7)
r = size(sysr,2);
@test gnrank(sysd-sysx*sysr,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      isproper(sysr) && # checking properness of factors
      length(zer) == 0  && # checking lack of zeros
      info.nrank == r && ismissing(info.nfuz) && info.niuz == 0

@time sysr, sysx, info = grange(sysd, zeros = "none", inner = true, atol = 1.e-7)
zer = gzero(sysr,atol1=1.e-7)
r = size(sysr,2);
@test gnrank(sysd-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysr'*sysr-I,atol=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysr) && # checking properness of factors
      isstable(sysr) && # checking stability of inner factor
      length(zer) == 0  && # checking number of zeros
      info.nrank == r && ismissing(info.nfuz) && info.niuz == 0

@time sysr, sysx, info = gcrange(sysd, zeros = "none", coinner = true, atol = 1.e-7)
zer = gzero(sysr,atol1=1.e-7)
r = size(sysr,2);
@test gnrank(sysd-sysx*sysr,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysr*sysr'-I,atol=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysr) && # checking properness of factors
      isstable(sysr) && # checking stability of inner factor
      length(zer) == 0  && # checking number of zeros
      info.nrank == r && ismissing(info.nfuz) && info.niuz == 0

@time sysr, sysx, info = grange(sysd, zeros = "unstable", inner = false, atol = 1.e-7)
zer = gzero(sysr,atol1=1.e-7)
r = size(sysr,2);
@test gnrank(sysd-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      isproper(sysr) && # checking properness of factors
      sort(zeref[abs.(zeref) .>= 1]) ≈ sort(zer) && # check zeros
      info.nrank == r && info.nfuz == 0 && info.niuz == 0

@time sysr, sysx, info = gcrange(sysd, zeros = "unstable", coinner = false, atol = 1.e-7)
zer = gzero(sysr,atol1=1.e-7)
r = size(sysr,2);
@test gnrank(sysd-sysx*sysr,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      isproper(sysr) && # checking properness of factors
      sort(zeref[abs.(zeref) .>= 1]) ≈ sort(zer) && # check zeros
      info.nrank == r && info.nfuz == 0 && info.niuz == 0

@time sysr, sysx, info = grange(sysd, zeros = "unstable", inner = true, atol = 1.e-7)
zer = gzero(sysr,atol1=1.e-7)
r = size(sysr,2);
@test gnrank(sysd-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysr'*sysr-I,atol=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysr) && # checking properness of factors
      isstable(sysr) && # checking stability of inner factor
      sort(zeref[abs.(zeref) .>= 1]) ≈ sort(zer) && # check zeros
      info.nrank == r && info.nfuz == 0 && info.niuz == 0

@time sysr, sysx, info = gcrange(sysd, zeros = "unstable", coinner = true, atol = 1.e-7)
zer = gzero(sysr,atol1=1.e-7)
r = size(sysr,2);
@test gnrank(sysd-sysx*sysr,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysr*sysr'-I,atol=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysr) && # checking properness of factors
      isstable(sysr) && # checking stability of inner factor
      sort(zeref[abs.(zeref) .>= 1]) ≈ sort(zer) && # check zeros
      info.nrank == r && info.nfuz == 0 && info.niuz == 0

@time sysr, sysx, info = grange(sysd, zeros = "all", atol = 1.e-7)
zer = gzero(sysr,atol1=1.e-7)
r = size(sysr,2);
@test gnrank(sysd-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      isproper(sysr) && # checking properness of factors
      sort(real(zeref[isfinite.(zeref)])) ≈ sort(real(zer[isfinite.(zer)])) &&
      zeref[isinf.(zeref)] == zer[isinf.(zer)] && norm(imag(zer),Inf) < 1.e-7 && # check zeros
      info.nrank == r && ismissing(info.nfuz) && ismissing(info.niuz)

@time sysr, sysx, info = grange(sysd, zeros = "finite", atol = 1.e-7)
zer = gzero(sysr,atol1=1.e-7)
r = size(sysr,2);
@test gnrank(sysd-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      isproper(sysr) && # checking properness of factors
      sort(zeref[isfinite.(zeref)]) ≈ sort(real(zer)) && norm(imag(zer),Inf) < 1.e-7 && # check zeros
      info.nrank == r && ismissing(info.nfuz) && ismissing(info.niuz)

@time sysr, sysx, info = grange(sysd, zeros = "infinite", atol = 1.e-7)
zer = gzero(sysr,atol1=1.e-7)
r = size(sysr,2);
@test gnrank(sysd-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      isproper(sysr) && # checking properness of factors
      sort(zeref[isinf.(zeref)]) ≈ sort(real(zer)) && norm(imag(zer),Inf) < 1.e-7 && # check zeros
      info.nrank == r && ismissing(info.nfuz) && info.niuz == 0

@time sysr, sysx, info = grange(sysd, zeros = "s-unstable", atol = 1.e-7)
zer = gzero(sysr,atol1=1.e-7)
r = size(sysr,2);
@test gnrank(sysd-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      isproper(sysr) && # checking properness of factors
      sort(zeref[abs.(zeref) .>= 1]) ≈ sort(zer) && # check zeros
      info.nrank == r && info.nfuz == 0 && info.niuz == 0

@time sysr, sysx, info = grange(sysd, zeros = "stable", atol = 1.e-7)
zer = gzero(sysr,atol1=1.e-7)
r = size(sysr,2);
@test gnrank(sysd-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      isproper(sysr) && # checking properness of factors
      norm(sort(zeref[abs.(zeref) .< 1]) - sort(abs.(zer)),Inf) < 1.e-7 && # check zeros
      info.nrank == r && info.nfuz == 0 && info.niuz == 0

s = Polynomial([0, 1],'s');
gn = [(s-1) s 1;
    0 (s-2) (s-2);
    (s-1) (s^2+2*s-2) (2*s-1)];
gd = [(s+2) (s+2) (s+2);
    1 (s+1)^2 (s+1)^2;
    (s+2) (s+1)*(s+2) (s+1)*(s+2)];

@time sys = dss(gn,gd,minimal = true, atol = 1.e-7);
@time sysr, sysx, info = grange(sys, zeros = "unstable", inner = true, atol = 1.e-7)
r = size(sysr,2);
zer = gzero(gminreal(sysx),atol1=1.e-7);
@test gnrank(sys-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) = Gi(s)*Go(s) = 0
      gnrank(sysr'*sysr-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysr) & isproper(sysx) && # checking properness of factors
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      count(t -> isinf(t), zer) == info.niuz &&
      count(t -> (real.(t) .>= 0 && isfinite(t)), zer) == 0 &&
      info.nrank == r && info.nfuz == 0 && info.niuz == 1

@time sysr, sysx, info = grange(sys, inner = true, atol = 1.e-7)
zeref = gzero(gminreal(sys),atol1=1.e-7)
r = size(sysr,2);
zer = gzero(gminreal(sysx),atol1=1.e-7);
@test gnrank(sys-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) = Gi(s)*Go(s) = 0
      gnrank(sysr'*sysr-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysr) & isproper(sysx) && # checking properness of factors
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      count(t -> isinf(t), zer) == info.niuz &&
      count(t -> (real.(t) .>= 0 && isfinite(t)), zer) == count(t -> (real.(t) .>= 0 && isfinite(t)), zeref) &&
      info.nrank == r && ismissing(info.nfuz) && info.niuz == 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];

sys = dss(g,minimal = true, atol = 1.e-7);
@time sysr, sysx, info = grange(sys, zeros = "unstable", inner = true, atol = 1.e-7)
zeref = gzero(sys,atol1=1.e-7)
r = size(sysx,1);
zer = gzero(gminreal(sysx),atol1=1.e-7)
@test gnrank(sys-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) = Gi(s)*Go(s) = 0
      gnrank(sysr'*sysr-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysr) & (isproper(sys) ? isproper(sysx) : true) && # checking properness of factors
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      count(t -> isinf(t), zer) == info.niuz &&
      count(t -> (real.(t) .<= 0 && isfinite(t)), zer) >= count(t -> (isfinite(t)), zeref) &&
      info.nrank == r && info.nfuz == 0 && info.niuz == 0

z = Polynomial([0, 1],'z')
g = [z^2+z+1 4*z^2+3*z+2 2*z^2-2;
    z 4*z-1 2*z-2;
    z^2 4*z^2-z 2*z^2-2*z];

sys = dss(g,minimal = true, atol = 1.e-7,Ts = 1)
@time sysr, sysx, info = grange(sys, zeros = "unstable", inner = true, atol = 1.e-7)
zeref = gzero(sys,atol1=1.e-7)
r = size(sysx,1)
zer = gzero(gminreal(sysx),atol1=1.e-7)
@test gnrank(sys-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(z) - Gi(z)*Go(z) = 0
      gnrank(sysr'*sysr-I,atol1=1.e-7) == 0  && # conj(Gi(z))*Gi(z)-I = 0
      isproper(sysr) && (isproper(sys) ? isproper(sysx) : true) && # checking properness of factors
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      count(t -> isinf(t), zer) == info.niuz &&
      count(t -> (abs.(t) .<= 1), zer) >= count(t -> (abs.(t) .>= 1), zeref) &&
      info.nrank == r && info.nfuz == 1 && info.niuz == 0

z = Polynomial([0, 1],'z')
# Eigenvalue(s) on the unit circle
gn = 1
gd = z-1
sys = dss(gn,gd,Ts = 1)
sysr, sysx, info = grange(sys, zeros = "unstable", inner = true, atol = 1.e-7)
@test gnrank(sys-sysr*sysx,atol1=1.e-7) == 0   &&
      gnrank(sysr'*sysr-I,atol1=1.e-7) == 0

z = Polynomial([0, 1],'z')
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];

sys = dss(gn,gd,minimal = true, atol = 1.e-7,Ts = 1);
sysr, sysx, info = grange(sys, zeros = "unstable", inner = true, atol = 1.e-7)
zeref = gzero(sys,atol1=1.e-7)
r = size(sysx,1)
zer = gzero(gminreal(sysx),atol1=1.e-7)
@test gnrank(sys-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(z) - Gi(z)*Go(z) = 0
      gnrank(sysr'*sysr-I,atol1=1.e-7) == 0  && # conj(Gi(z))*Gi(z)-I = 0
      isproper(sysr) && (isproper(sys) ? isproper(sysx) : true) && # checking properness of factors
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      count(t -> isinf(t), zer) == info.niuz &&
      count(t -> (abs.(t) .<= 1), zer) >= count(t -> (abs.(t) .>= 1), zeref) &&
      info.nrank == r && info.nfuz == 1 && info.niuz == 0

fast = true; Ty = Complex{Float64}; Ty = Float64
n = 5; p = 4; m = 2;

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

# continuous, standard
sys = rss(n,p,m,T = Ty,disc=false,stable=true);sys.D[:,:] = 0*sys.D;
sys = sys*rss(n,m,m,T = Ty,disc=false,stable=true)
@time sysr, sysx, info = grange(sys, zeros = "unstable", inner = true, fast = fast, atol = 1.e-7) ;
zeref = gzero(sys,atol1=1.e-7)
r = size(sysx,1);
zer = gzero(gminreal(sysx),atol1=1.e-7)
@test gnrank(sys-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysr'*sysr-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysr) && (isproper(sys) ? isproper(sysx) : true) && # checking properness of factors
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      count(t -> isinf(t), zer) == info.niuz &&
      count(t -> (real.(t) .<= 0 && isfinite(t)), zer) >= count(t -> (isfinite(t)), zeref) &&
      info.nrank == r && info.nfuz == 0 && info.niuz == m

# discrete, standard
sys = rss(n,p,m,T = Ty,disc=true,stable=true);
@time sysr, sysx, info = grange(sys, zeros = "unstable", inner = true, fast = fast, atol = 1.e-7) ;
zeref = gzero(sys,atol1=1.e-7)
r = size(sysx,1);
zer = gzero(gminreal(sysx),atol1=1.e-7)
@test gnrank(sys-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(z) - Gi(z)*Go(z) = 0
      gnrank(sysr'*sysr-I,atol1=1.e-7) == 0  && # conj(Gi(u))*Gi(z)-I = 0
      isproper(sysr) && (isproper(sys) ? isproper(sysx) : true) && # checking properness of factors
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      count(t -> isinf(t), zer) == info.niuz &&
      count(t -> (abs.(t) .<= 1), zer) >= count(t -> (abs.(t) .>= 1), zeref) &&
      info.nrank == r && info.nfuz == 0 && info.niuz == 0

# continuous, descriptor, no infinite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=false,stable=true);
@time sysr, sysx, info = grange(sys, zeros = "unstable", inner = true, fast = fast, atol = 1.e-7) ;
zeref = gzero(sys,atol1=1.e-7)
r = size(sysx,1);
zer = gzero(gminreal(sysx),atol1=1.e-7)
@test gnrank(sys-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysr'*sysr-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysr) && (isproper(sys) ? isproper(sysx) : true) && # checking properness of factors
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      count(t -> isinf(t), zer) == info.niuz &&
      count(t -> (real.(t) .<= 0 && isfinite(t)), zer) >= count(t -> (isfinite(t)), zeref) &&
      info.nrank == r && info.nfuz == 0 && info.niuz == 0

# discrete, descriptor, no infinite eigenvalues
sys = rdss(n,p,m,T = Ty, disc=true,stable=true);
@time sysr, sysx, info = grange(sys, zeros = "unstable", inner = true, fast = fast, atol = 1.e-7) ;
zeref = gzero(sys,atol1=1.e-7)
r = size(sysx,1);
zer = gzero(gminreal(sysx),atol1=1.e-7)
@test gnrank(sys-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(z) - Gi(z)*Go(z) = 0
      gnrank(sysr'*sysr-I,atol1=1.e-7) == 0  && # conj(Gi(u))*Gi(z)-I = 0
      isproper(sysr) && (isproper(sys) ? isproper(sysx) : true) && # checking properness of factors
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      count(t -> isinf(t), zer) == info.niuz &&
      count(t -> (abs.(t) .<= 1), zer) >= count(t -> (abs.(t) .>= 1), zeref) &&
      info.nrank == r && info.nfuz == 0 && info.niuz == 0

# continuous, descriptor, proper
sys = rdss(n,p,m,T = Ty,disc=false,id=ones(Int,3),stable=true);
@time sysr, sysx, info = grange(sys, zeros = "unstable", inner = true, fast = fast, atol = 1.e-7) ;
zeref = gzero(sys,atol1=1.e-7)
r = size(sysx,1);
zer = gzero(gminreal(sysx),atol1=1.e-7)
@test gnrank(sys-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysr'*sysr-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysr) && (isproper(sys) ? isproper(sysx) : true) && # checking properness of factors
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      count(t -> isinf(t), zer) == info.niuz &&
      count(t -> (real.(t) .<= 0 && isfinite(t)), zer) >= count(t -> (isfinite(t)), zeref) &&
      info.nrank == r && info.nfuz == 0 && info.niuz == 0

# continuous, descriptor, proper
sys = rdss(n,p,m,T = Ty,disc=false,id=ones(Int,3),stable=true);
@time sysr, sysx, info = grange(sys, zeros = "unstable", inner = true, fast = fast, atol = 1.e-7) ;
zeref = gzero(sys,atol1=1.e-7)
r = size(sysx,1);
zer = gzero(gminreal(sysx),atol1=1.e-7)
@test gnrank(sys-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysr'*sysr-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysr) && (isproper(sys) ? isproper(sysx) : true) && # checking properness of factors
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      count(t -> isinf(t), zer) == info.niuz &&
      count(t -> (real.(t) .<= 0 && isfinite(t)), zer) >= count(t -> (isfinite(t)), zeref) &&
      info.nrank == r && info.nfuz == 0 && info.niuz == 0

# discrete, descriptor, proper
sys = rdss(n,p,m,T = Ty, disc=true,id=ones(Int,3));
@time sysr, sysx, info = grange(sys, zeros = "unstable", inner = true, fast = fast, atol = 1.e-7) ;
zeref = gzero(sys,atol1=1.e-7)
r = size(sysx,1);
zer = gzero(gminreal(sysx),atol1=1.e-7)
@test gnrank(sys-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(z) - Gi(z)*Go(z) = 0
      gnrank(sysr'*sysr-I,atol1=1.e-7) == 0  && # conj(Gi(u))*Gi(z)-I = 0
      isproper(sysr) && (isproper(sys) ? isproper(sysx) : true) && # checking properness of factors
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      count(t -> isinf(t), zer) == info.niuz &&
      count(t -> (abs.(t) .<= 1), zer) >= count(t -> (abs.(t) .>= 1), zeref) &&
      info.nrank == r && info.nfuz == 0 && info.niuz == 0

# continuous, descriptor, infinite poles
sys = rdss(n,p,m,T = Ty,disc=false,id=[3*ones(Int,1);2*ones(Int,1)]);
@time sysr, sysx, info = grange(sys, zeros = "unstable", inner = true, fast = fast, atol = 1.e-7) ;
zeref = gzero(sys,atol1=1.e-7)
r = size(sysx,1);
zer = gzero(gminreal(sysx),atol1=1.e-7)
@test gnrank(sys-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysr'*sysr-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysr) && (isproper(sys) ? isproper(sysx) : true) && # checking properness of factors
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      count(t -> isinf(t), zer) == info.niuz &&
      count(t -> (real.(t) .<= 0 && isfinite(t)), zer) >= count(t -> (isfinite(t)), zeref) &&
      info.nrank == r && info.nfuz == 0 && info.niuz == 0

# discrete, descriptor, infinite poles
sys = rdss(n,p,m,T = Ty, disc=true,id=[3*ones(Int,1);2*ones(Int,1)],stable=true);
@time sysr, sysx, info = grange(sys, zeros = "unstable", inner = true, fast = fast, atol = 1.e-7) ;
zeref = gzero(sys,atol1=1.e-7)
r = size(sysx,1);
zer = gzero(gminreal(sysx),atol1=1.e-7)
@test gnrank(sys-sysr*sysx,atol1=1.e-7) == 0   #&&   #  G(z) - Gi(z)*Go(z) = 0
@test       iszero(sysr'*sysr-I,atol1=1.e-7)  #&& # conj(Gi(u))*Gi(z)-I = 0
@test       isproper(sysr) && (isproper(sys) ? isproper(sysx) : true) #&& # checking properness of factors
@test       isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) #&&
@test       count(t -> isinf(t), zer) == info.niuz #&&
@test       count(t -> (abs.(t) .<= 1), zer) >= count(t -> (abs.(t) .>= 1), zeref) #&&
@test       info.nrank == r && info.nfuz == 0 && info.niuz == 0

# continuous, descriptor, proper, uncontrollable infinite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=false,iduc=[3*ones(Int,1);2*ones(Int,1)]);
@time sysr, sysx, info = grange(sys, zeros = "unstable", inner = true, fast = fast, atol = 1.e-7) ;
zeref = gzero(sys,atol1=1.e-7)
r = size(sysx,1);
zer = gzero(gminreal(sysx),atol1=1.e-7)
@test gnrank(sys-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysr'*sysr-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysr) && (isproper(sys) ? isproper(sysx) : true) && # checking properness of factors
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      count(t -> isinf(t), zer) == info.niuz &&
      count(t -> (real.(t) .<= 0 && isfinite(t)), zer) >= count(t -> (isfinite(t)), zeref) &&
      info.nrank == r && info.nfuz == 0 && info.niuz == 0

# discrete, descriptor, proper, uncontrollable infinite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=true,iduc=[3*ones(Int,1);2*ones(Int,1)]);
@time sysr, sysx, info = grange(sys, zeros = "unstable", inner = true, fast = fast, atol = 1.e-7) ;
zeref = gzero(gminreal(sys),atol1=1.e-7,atol2=1.e-7)
r = size(sysx,1);
zer = gzero(gminreal(sysx),atol1=1.e-7,atol2=1.e-7)
@test gnrank(sys-sysr*sysx,atol1=1.e-7) == 0   &&   #  G(z) - Gi(z)*Go(z) = 0
      iszero(sysr'*sysr-I,atol1=1.e-5)  && # conj(Gi(u))*Gi(z)-I = 0
      isproper(sysr) && (isproper(sys) ? isproper(sysx) : true) && # checking properness of factors
      isstable(sysr) && (isstable(sys) ? isstable(sysx) : true) &&
      count(t -> isinf(t), zer) == info.niuz &&
      count(t -> (abs.(t) .<= 1), zer) >= count(t -> (abs.(t) .>= 1), zeref) &&
      info.nrank == r && info.nfuz == 0 && info.niuz == 0

end
end

end # grange

@testset "grnull" begin

p1 = 0; p2 = 0; m = 0;
sys = dss(ones(p1+p2,m));
@time sysnull, info = grnull(sys,p2)
isys1 = 1:p1; isys2 = p1+1:p1+p2;
indf = m+1:m+p2; indnull = 1:m;
@test gnrank(sys[isys1,:]*sysnull[indnull,:]) == 0 &&
      gnrank(sys[isys2,:]*sysnull[indnull,:]-sysnull[indf,:], atol = 1.e-7) == 0 &&
      info.nrank == 0 && info.degs == [] && info.stdim == []

p1 = 0; p2 = 0; m = 1;
sys = dss(ones(p1+p2,m));
@time sysnull, info = grnull(sys,p2)
isys1 = 1:p1; isys2 = p1+1:p1+p2;
indf = m+1:m+p2; indnull = 1:m;
@test gnrank(sys[isys1,:]*sysnull[indnull,:]) == 0 &&
      gnrank(sys[isys2,:]*sysnull[indnull,:]-sysnull[indf,:], atol = 1.e-7) == 0 &&
      info.nrank == 0 && info.degs == [0] && info.stdim == []

p1 = 1; p2 = 1; m = 1;
sys = dss(ones(p1+p2,m));
@time sysnull, info = grnull(sys,p2)
isys1 = 1:p1; isys2 = p1+1:p1+p2;
indf = m+1:m+p2; indnull = 1:m
@test gnrank(sys[isys1,:]*sysnull[indnull,:]) == 0 &&
      gnrank(sys[isys2,:]*sysnull[indnull,:]-sysnull[indf,:], atol = 1.e-7) == 0 &&
      info.nrank == 1 && info.degs == Int[]  && info.stdim == []

p1 = 1; p2 = 1; m = 2;
sys = dss(ones(p1+p2,m));
@time sysnull, info = grnull(sys,p2)
isys1 = 1:p1; isys2 = p1+1:p1+p2;
indf = m+1:m+p2; indnull = 1:m;
@test gnrank(sys[isys1,:]*sysnull[indnull,:], atol = 1.e-7) == 0 &&
      gnrank(sys[isys2,:]*sysnull[indnull,:]-sysnull[indf,:], atol = 1.e-7) == 0 &&
      info.nrank == 1 && info.degs == [0] && info.stdim == Int64[]

p1 = 0; p2 = 1; m = 2;
sys = dss(ones(p1+p2,m));
@time sysnull, info = grnull(sys,p2)
isys1 = 1:p1; isys2 = p1+1:p1+p2;
indf = m+1:m+p2; indnull = 1:m;
@test gnrank(sys[isys1,:]*sysnull[indnull,:]) == 0 &&
      gnrank(sys[isys2,:]*sysnull[indnull,:]-sysnull[indf,:], atol = 1.e-7) == 0 &&
      info.nrank == 0 && info.degs == [0, 0] && info.stdim == []

p1 = 1; p2 = 1; m = 2;
sys = dss([1 0; 0 1]);
@time sysnull, info = grnull(sys,p2)
isys1 = 1:p1; isys2 = p1+1:p1+p2;
indf = m+1:m+p2; indnull = 1:m;
@test gnrank(sys[isys1,:]*sysnull[indnull,:]) == 0 &&
      gnrank(sys[isys2,:]*sysnull[indnull,:]-sysnull[indf,:], atol = 1.e-7) == 0 &&
      info.nrank == 1 && info.degs == [0] && info.stdim == []

# uncontrollable realization
s = Polynomial([0, 1],'s')
m = 2; p1 = 0; p2 = 1;
G = dss(zeros(0,m)); F = dss([1 1], [(s+1) (s+2)]);
sys = [G; F]

@time nr, info = grnull(sys,p2; atol = 1.e-7); info

@test gnrank(sys[1:p1,:]*nr[1:m,:],atol=1.e-7) == 0 &&
      gnrank(sys*nr[1:m,:]-nr[m+1:m+p2,:],atol= 1.e-7) == 0 &&
      info.nrank == 0 && info.degs == [1, 1] && info.stdim == [2]

# uncontrollable realization
s = Polynomial([0, 1],'s')
m = 2; p1 = 0; p2 = 1;
G = dss(zeros(0,m)); F = dss([1 1], [(s+1) (s+2)]);
sys = [G; F]

@time nr, info = grnull(sys,p2; atol = 1.e-7, simple = true); info

@test gnrank(sys[1:p1,:]*nr[1:m,:],atol=1.e-7) == 0 &&
      gnrank(sys*nr[1:m,:]-nr[m+1:m+p2,:],atol= 1.e-7) == 0 &&
      info.nrank == 0 && info.degs == [1, 1] && info.stdim == [1, 1]

# uncontrollable realization
s = Polynomial([0, 1],'s')
m = 2; p1 = 0; p2 = 1;
G = dss(zeros(0,m)); F = dss([1 1], [(s-1) (s+2)]);
sys = [G; F]

@time nr, info = grnull(sys,p2; atol = 1.e-7, inner = true); info

@test gnrank(sys[1:p1,:]*nr[1:m,:],atol=1.e-7) == 0 &&
      gnrank(sys*nr[1:m,:]-nr[m+1:m+p2,:],atol= 1.e-7) == 0 &&
      gnrank(nr[1:m,:]'*nr[1:m,:]-I,atol= 1.e-7) == 0 &&
      info.nrank == 0 && info.degs == [1, 1] && info.stdim == [2]

# Example 3, Zuniga & Henrion AMC (2009)

s = Polynomial([0, 1],'s');
m = 5; p1 = 2; p2 = 0;
G = copy(transpose([ 1 2; -2 s; 2 2*s; -2*s^2 s^3; -2+s^3 s+s^3]));
sys = dss(G[:,1:m]);

# minimal rational basis
@time nr, info = grnull(sys,p2,atol=1.e-7); info
@test gnrank(sys[1:p1,:]*nr[1:m,:],atol=1.e-7) == 0 &&
      gnrank(sys[p1+1:p1+p2,:]*nr[1:m,:]-nr[m+1:m+p2,:],atol= 1.e-7) == 0 &&
      info.nrank == 2 && info.degs == [1, 2, 3] && info.stdim == [3, 2, 1]

# minimal polynomial basis
@time nr, info = grnull(sys,p2;atol=1.e-7,polynomial = true); info
@test gnrank(sys[1:p1,:]*nr[1:m,:],atol=1.e-7) == 0 &&
      gnrank(sys[p1+1:p1+p2,:]*nr[1:m,:]-nr[m+1:m+p2,:],atol= 1.e-7) == 0 &&
      all(isinf.(gpole(nr,atol=1.e-7))) &&
      info.nrank == 2 && info.degs == [1, 2, 3] && info.stdim == [3, 2, 1]

# minimal inner rational basis
@time nr, info = grnull(sys,p2;atol=1.e-7, inner = true); info
@test gnrank(sys[1:p1,:]*nr[1:m,:],atol=1.e-7) == 0 &&
      gnrank(sys[p1+1:p1+p2,:]*nr[1:m,:]-nr[m+1:m+p2,:],atol= 1.e-7) == 0 &&
      gnrank(nr[1:m,:]'*nr[1:m,:]-I,atol= 1.e-7) == 0 &&
      info.nrank == 2 && info.degs == [1, 2, 3] && info.stdim == [3, 2, 1]

z = Polynomial([0, 1],'z');
m = 5; p1 = 2; p2 = 0;
G = copy(transpose([ 1 2; -2 z; 2 2*z; -2*z^2 z^3; -2+z^3 z+z^3]));
sys = dss(G,Ts = 1);

# minimal inner rational basis (discrete-time)
@time nr, info = grnull(sys,atol=1.e-7, inner = true); info
@test gnrank(sys[1:p1,:]*nr[1:m,:],atol=1.e-7) == 0 &&
      gnrank(sys[p1+1:p1+p2,:]*nr[1:m,:]-nr[m+1:m+p2,:],atol= 1.e-7) == 0 &&
      gnrank(nr[1:m,:]'*nr[1:m,:]-I,atol= 1.e-7) == 0 &&
      info.nrank == 2 && info.degs == [1, 2, 3] && info.stdim == [3, 2, 1]

z = Polynomial([0, 1],'z');
m = 5; p1 = 2; p2 = 0;
G = copy(transpose([ 1 2; -2 z; 2 2*z; -2*z^2 z^3; -2+z^3 z+z^3]));
sys = dss(G,Ts = 1);

# minimal rational simple basis (discrete-time)
@time nr, info = grnull(sys,atol=1.e-7, inner = false, simple = true); info
@test gnrank(sys[1:p1,:]*nr[1:m,:],atol=1.e-7) == 0 &&
      gnrank(sys[p1+1:p1+p2,:]*nr[1:m,:]-nr[m+1:m+p2,:],atol= 1.e-7) == 0 &&
      info.nrank == 2 && info.degs == [1, 2, 3] && info.stdim == [1, 2, 3]

z = Polynomial([0, 1],'z');
m = 5; p1 = 2; p2 = 0;
G = copy(transpose([ 1 2; -2 z; 2 2*z; -2*z^2 z^3; -2+z^3 z+z^3]));
sys = dss(G,Ts = 1);

# minimal rational simple basis with innner vectors (discrete-time)
@time nr, info = grnull(sys,atol=1.e-7, inner = true, simple = true); info
@test gnrank(sys[1:p1,:]*nr[1:m,:],atol=1.e-7) == 0 &&
      gnrank(sys[p1+1:p1+p2,:]*nr[1:m,:]-nr[m+1:m+p2,:],atol= 1.e-7) == 0 &&
      gnrank(nr[:,1]'*nr[:,1]-I,atol= 1.e-7) == 0 &&
      gnrank(nr[:,2]'*nr[:,2]-I,atol= 1.e-7) == 0 &&
      gnrank(nr[:,3]'*nr[:,3]-I,atol= 1.e-7) == 0 &&
      info.nrank == 2 && info.degs == [1, 2, 3] && info.stdim == [1, 2, 3]

z = Polynomial([0, 1],'z');
m = 5; p1 = 2; p2 = 0;
G = copy(transpose([ 1 2; -2 z; 2 2*z; -2*z^2 z^3; -2+z^3 z+z^3]));
sys = dss(G,Ts = 1);

# minimal polynomial simple basis (discrete-time)
@time nr, info = grnull(sys,atol=1.e-7, polynomial = true, simple = false); info
@test gnrank(sys[1:p1,:]*nr[1:m,:],atol=1.e-7) == 0 &&
      gnrank(sys[p1+1:p1+p2,:]*nr[1:m,:]-nr[m+1:m+p2,:],atol= 1.e-7) == 0 &&
      all(isinf.(gpole(nr,atol=1.e-7))) &&
      info.nrank == 2 && info.degs == [1, 2, 3] && info.stdim == [3, 2, 1]

# minimal polynomial simple basis (discrete-time)
@time nr, info = grnull(sys,atol=1.e-7, polynomial = true, simple = true); info
@test gnrank(sys[1:p1,:]*nr[1:m,:],atol=1.e-7) == 0 &&
      gnrank(sys[p1+1:p1+p2,:]*nr[1:m,:]-nr[m+1:m+p2,:],atol= 1.e-7) == 0 &&
      all(isinf.(gpole(nr,atol=1.e-7))) &&
      info.nrank == 2 && info.degs == [1, 2, 3] && info.stdim == [2, 3, 4]

z = Polynomial([0, 1],'z');
m = 5; p1 = 2; p2 = 0;
G = copy(transpose([ 1 2; -2 z; 2 2*z; -2*z^2 z^3; -2+z^3 z+z^3]));
sys = dss(G,Ts = 1);

# minimal simple basis with poles assigned to 0 (discrete-time)
@time nr, info = grnull(sys,atol=1.e-7, poles = [0,0,0], simple = true); info
@test gnrank(sys[1:p1,:]*nr[1:m,:],atol=1.e-7) == 0 &&
      gnrank(sys[p1+1:p1+p2,:]*nr[1:m,:]-nr[m+1:m+p2,:],atol= 1.e-7) == 0 &&
      all(abs.(gpole(nr,atol=1.e-7)).< 0.001) &&
      info.nrank == 2 && info.degs == [1, 2, 3] && info.stdim == [1, 2, 3]

z = Polynomial([0, 1],'z');
m = 5; p1 = 2; p2 = 0;
G = copy(transpose([ 1 2; -2 z; 2 2*z; -2*z^2 z^3; -2+z^3 z+z^3]));
sys = dss(G,Ts = 1);

# minimal rational basis with poles assigned to 0 (discrete-time)
@time nr, info = grnull(sys,atol=1.e-7, sdeg = 0, simple = false); info
@test gnrank(sys[1:p1,:]*nr[1:m,:],atol=1.e-7) == 0 &&
      gnrank(sys[p1+1:p1+p2,:]*nr[1:m,:]-nr[m+1:m+p2,:],atol= 1.e-7) == 0 &&
      all(abs.(gpole(nr,atol=1.e-7)).< 0.01) &&
      info.nrank == 2 && info.degs == [1, 2, 3] && info.stdim == [3, 2, 1]

# Kailath 1980, pag. 459, rank 2 matrix with both left and right nullspaces
s = Polynomial([0, 1],'s');
m = 4; p1 = 3; p2 = 0;
gn = [1 0 1 s;
    0 (s+1)^2 (s+1)^2 0;
    -1 (s+1)^2 s^2+2*s -s^2];
gd = [s 1 s 1;
    1 1 1 1;
    1 1 1 1];
sys = dss(gn,gd,minimal = true, atol = 1.e-7);

# minimal rational basis with poles assigned to -1
@time nr, info = grnull(sys,atol=1.e-7,poles = [-1,-1]); info
@test gnrank(sys[1:p1,:]*nr[1:m,:],atol=1.e-7) == 0 &&
      gnrank(sys[p1+1:p1+p2,:]*nr[1:m,:]-nr[m+1:m+p2,:],atol= 1.e-7) == 0 &&
      all(real(gpole(nr).+1) .< 2*sqrt(eps())) && all(abs.(imag(gpole(nr))) .< 0.000001) &&
      info.nrank == 2 && info.degs == [0, 2] && info.stdim == [1, 1]

# minimal simple basis with poles assigned to -1
@time nr, info = grnull(sys,atol=1.e-7,poles = [-1,-1], simple = true); info
@test gnrank(sys[1:p1,:]*nr[1:m,:],atol=1.e-7) == 0 &&
      gnrank(sys[p1+1:p1+p2,:]*nr[1:m,:]-nr[m+1:m+p2,:],atol= 1.e-7) == 0 &&
      all(real(gpole(nr)) .≈ -1) && all(abs.(imag(gpole(nr))) .< 0.000001) &&
      info.nrank == 2 && info.degs == [0, 2] && info.stdim == [0, 2]

# minimal simple polynomial basis
@time nr, info = grnull(sys,atol=1.e-7,polynomial = true, simple = true); info
@test gnrank(sys[1:p1,:]*nr[1:m,:],atol=1.e-7) == 0 &&
      gnrank(sys[p1+1:p1+p2,:]*nr[1:m,:]-nr[m+1:m+p2,:],atol= 1.e-7) == 0 &&
      all(isinf.(gpole(nr,atol=1.e-7))) &&
      info.nrank == 2 && info.degs == [0, 2] && info.stdim == [0, 3]

# random examples
m = 5; p1 = 2; p2 = 2; p = p1+p2; n = 5;

Ty = Float64; fast = true; simple = false; polynomial = false;
for Ty in (Float64, Complex{Float64})

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

for fast in (true, false)

for polynomial in (false, true)

#non-simple basis
# continuous, standard
sys = rss(n,p,m,T = Ty, disc=false);
@time nr, info = grnull(sys,p2, atol=1.e-7,polynomial = polynomial, simple = false); info
@test gnrank(sys[1:p1,:]*nr[1:m,:],atol=1.e-7) == 0 &&
      gnrank(sys[p1+1:p1+p2,:]*nr[1:m,:]-nr[m+1:m+p2,:],atol= 1.e-7) == 0 &&
      info.nrank == p1 && info.degs == [1, 2, 2] && info.stdim == [3, 2]

# discrete, descriptor
sys = rdss(n,p,m,T = Ty, disc=true);
@time nr, info = grnull(sys,p2, atol=1.e-7,polynomial = polynomial, simple = false); info
@test gnrank(sys[1:p1,:]*nr[1:m,:],atol=1.e-7) == 0 &&
      gnrank(sys[p1+1:p1+p2,:]*nr[1:m,:]-nr[m+1:m+p2,:],atol= 1.e-7) == 0 &&
      info.nrank == p1 && info.degs == [1, 2, 2] && info.stdim == [3, 2]

# discrete, descriptor
sys = rdss(n,p,m,T = Ty, disc=true);
@time nr, info = grnull(sys,p2, atol=1.e-7,polynomial = polynomial, inner = true, simple = false); info
@test gnrank(sys[1:p1,:]*nr[1:m,:],atol=1.e-7) == 0 &&
      gnrank(sys[p1+1:p1+p2,:]*nr[1:m,:]-nr[m+1:m+p2,:],atol= 1.e-7) == 0 &&
      (polynomial || gnrank(nr[1:m,:]'*nr[1:m,:]-I,atol= 1.e-7) == 0) &&
      info.nrank == p1 && info.degs == [1, 2, 2] && info.stdim == [3, 2]

# discrete, polynomial
sys = rdss(5,p,m,T = Ty, disc=true, id=[3*ones(Int,1);2*ones(Int,1)]);
@time nr, info = grnull(sys,p2, atol=1.e-7,polynomial = polynomial, simple = false); info
@test gnrank(sys[1:p1,:]*nr[1:m,:],atol=1.e-7) == 0 &&
      gnrank(sys[p1+1:p1+p2,:]*nr[1:m,:]-nr[m+1:m+p2,:],atol= 1.e-7) == 0 &&
      info.nrank == p1 && info.degs == [2, 3, 3] && info.stdim == [3, 3, 2]

# simple basis
# continuous, standard
sys = rss(n,p,m,T = Ty, disc=false);
@time nr, info = grnull(sys,p2, atol=1.e-7,polynomial = polynomial, simple = true); info
@test gnrank(sys[1:p1,:]*nr[1:m,:],atol=1.e-7) == 0 &&
      gnrank(sys[p1+1:p1+p2,:]*nr[1:m,:]-nr[m+1:m+p2,:],atol= 1.e-7) == 0 &&
      info.nrank == p1 && info.degs == [1, 2, 2] && info.stdim == (polynomial ? [2, 3, 3] : [1, 2, 2])

# discrete, descriptor
sys = rdss(n,p,m,T = Ty, disc=true);
@time nr, info = grnull(sys,p2, atol=1.e-7,polynomial = polynomial, simple = true); info
@test gnrank(sys[1:p1,:]*nr[1:m,:],atol=1.e-7) == 0 &&
      gnrank(sys[p1+1:p1+p2,:]*nr[1:m,:]-nr[m+1:m+p2,:],atol= 1.e-7) == 0 &&
      info.nrank == p1 && info.degs == [1, 2, 2] && info.stdim == (polynomial ? [2, 3, 3] : [1, 2, 2])

# discrete, descriptor
sys = rdss(n,p,m,T = Ty, disc=true);
@time nr, info = grnull(sys,p2, atol=1.e-7,polynomial = polynomial, inner = true, simple = true); info
@test gnrank(sys[1:p1,:]*nr[1:m,:],atol=1.e-7) == 0 &&
      gnrank(sys[p1+1:p1+p2,:]*nr[1:m,:]-nr[m+1:m+p2,:],atol= 1.e-7) == 0 &&
      (polynomial || gnrank(nr[1:m,1]'*nr[1:m,1]-I,atol= 1.e-7) == 0) &&
      (polynomial || gnrank(nr[1:m,2]'*nr[1:m,2]-I,atol= 1.e-7) == 0) &&
      info.nrank == p1 && info.degs == [1, 2, 2] && info.stdim == (polynomial ? [2, 3, 3] : [1, 2, 2])

# discrete, polynomial
sys = rdss(5,p,m,T = Ty, disc=true, id=[3*ones(Int,1);2*ones(Int,1)]);
@time nr, info = grnull(sys,p2, atol=1.e-7,polynomial = polynomial, simple = true); info
@test gnrank(sys[1:p1,:]*nr[1:m,:],atol=1.e-7) == 0 &&
      gnrank(sys[p1+1:p1+p2,:]*nr[1:m,:]-nr[m+1:m+p2,:],atol= 1.e-7) == 0 &&
       info.nrank == p1 && info.degs == [2, 3, 3] && info.stdim == (polynomial ? [3, 4, 4] : [2, 3, 3] )

end # polynomial

end # fast

end # Ty

end # grnull

@testset "glnull" begin

m1 = 0; m2 = 0; p = 0;
sys = dss(ones(p,m1+m2));
@time sysnull, info = glnull(sys,m2)
isys1 = 1:m1; isys2 = m1+1:m1+m2;
indf = p+1:p+m2; indnull = 1:p;
@test gnrank(sysnull[:,indnull]*sys[:,isys1]) == 0 &&
      gnrank(sysnull[:,indnull]*sys[:,isys2]-sysnull[:,indf], atol = 1.e-7) == 0 &&
      info.nrank == 0 && info.degs == [] && info.stdim == []

m1 = 0; m2 = 0; p = 1;
sys = dss(ones(p,m1+m2));
@time sysnull, info = glnull(sys,m2)
isys1 = 1:m1; isys2 = m1+1:m1+m2;
indf = p+1:p+m2; indnull = 1:p;
@test gnrank(sysnull[:,indnull]*sys[:,isys1]) == 0 &&
      gnrank(sysnull[:,indnull]*sys[:,isys2]-sysnull[:,indf], atol = 1.e-7) == 0 &&
      info.nrank == 0 && info.degs == [0] && info.stdim == []

m1 = 1; m2 = 1; p = 1;
sys = dss(ones(p,m1+m2));
@time sysnull, info = glnull(sys,m2)
isys1 = 1:m1; isys2 = m1+1:m1+m2;
indf = p+1:p+m2; indnull = 1:p;
@test gnrank(sysnull[:,indnull]*sys[:,isys1]) == 0 &&
      gnrank(sysnull[:,indnull]*sys[:,isys2]-sysnull[:,indf], atol = 1.e-7) == 0 &&
      info.nrank == 1 && info.degs == Int[]  && info.stdim == []

m1 = 1; m2 = 1; p = 2;
sys = dss(ones(p,m1+m2));
@time sysnull, info = glnull(sys,m2)
isys1 = 1:m1; isys2 = m1+1:m1+m2;
indf = p+1:p+m2; indnull = 1:p;
@test gnrank(sysnull[:,indnull]*sys[:,isys1], atol = 1.e-7) == 0 &&
      gnrank(sysnull[:,indnull]*sys[:,isys2]-sysnull[:,indf], atol = 1.e-7) == 0 &&
      info.nrank == 1 && info.degs == [0] && info.stdim == Int[]

m1 = 0; m2 = 1; p = 2;
sys = dss(ones(p,m1+m2));
@time sysnull, info = glnull(sys,m2)
isys1 = 1:m1; isys2 = m1+1:m1+m2;
indf = p+1:p+m2; indnull = 1:p;
@test gnrank(sysnull[:,indnull]*sys[:,isys1]) == 0 &&
      gnrank(sysnull[:,indnull]*sys[:,isys2]-sysnull[:,indf], atol = 1.e-7) == 0 &&
      info.nrank == 0 && info.degs == [0, 0] && info.stdim == []

m1 = 1; m2 = 1; p = 2;
sys = dss([1 0; 0 1]);
@time sysnull, info = glnull(sys,m2)
isys1 = 1:m1; isys2 = m1+1:m1+m2;
indf = p+1:p+m2; indnull = 1:p;
@test gnrank(sysnull[:,indnull]*sys[:,isys1]) == 0 &&
      gnrank(sysnull[:,indnull]*sys[:,isys2]-sysnull[:,indf], atol = 1.e-7) == 0 &&
      info.nrank == 1 && info.degs == [0] && info.stdim == []

# unobservable realization
s = Polynomial([0, 1],'s')
p = 2; m1 = 0; m2 = 1;
G = dss(zeros(p,0)); F = dss([1; 1], [(s+1); (s+2)]);
sys = [G F]

@time nl, info = glnull(sys,m2; atol = 1.e-7); info

@test gnrank(nl[:,1:p]*sys[:,1:m1],atol=1.e-7) == 0 &&
      gnrank(nl[:,1:p]*sys-nl[:,p+1:p+m2],atol= 1.e-7) == 0 &&
      info.nrank == 0 && info.degs == [1, 1] && info.stdim == [2]

# unobservable realization
s = Polynomial([0, 1],'s')
p = 2; m1 = 0; m2 = 1;
G = dss(zeros(p,0)); F = dss([1; 1], [(s+1); (s+2)]);
sys = [G F]

@time nl, info = glnull(sys,m2; atol = 1.e-7, simple = true); info

@test gnrank(nl[:,1:p]*sys[:,1:m1],atol=1.e-7) == 0 &&
      gnrank(nl[:,1:p]*sys-nl[:,p+1:p+m2],atol= 1.e-7) == 0 &&
      info.nrank == 0 && info.degs == [1, 1] && info.stdim == [1, 1]

# unobservable realization
s = Polynomial([0, 1],'s')
p = 2; m1 = 0; m2 = 1;
G = dss(zeros(p,0)); F = dss([1; 1], [(s-1); (s+2)]);
sys = [G F]

@time nl, info = glnull(sys,m2; atol = 1.e-7, coinner = true); info

@test gnrank(nl[:,1:p]*sys[:,1:m1],atol=1.e-7) == 0 &&
      gnrank(nl[:,1:p]*sys-nl[:,p+1:p+m2],atol= 1.e-7) == 0 &&
      gnrank(nl[:,1:p]*nl[:,1:p]'-I,atol= 1.e-7) == 0 &&
      info.nrank == 0 && info.degs == [1, 1] && info.stdim == [2]

# Example 3, Zuniga & Henrion AMC (2009)

s = Polynomial([0, 1],'s');
p = 5; m1 = 2; m2 = 0;
G = [ 1 2; -2 s; 2 2*s; -2*s^2 s^3; -2+s^3 s+s^3];
sys = dss(G[1:p,:]);

# minimal rational basis
@time nl, info = glnull(sys,m2,atol=1.e-7); info
@test gnrank(nl[:,1:p]*sys[:,1:m1],atol=1.e-7) == 0 &&
      gnrank(nl[:,1:p]*sys[:,m1+1:m1+m2]-nl[:,p+1:p+m2],atol= 1.e-7) == 0 &&
      info.nrank == 2 && info.degs == [1, 2, 3] && info.stdim == [1, 2, 3]

# minimal polynomial basis
@time nl, info = glnull(sys,m2,atol=1.e-7,polynomial = true); info
@test gnrank(nl[:,1:p]*sys[:,1:m1],atol=1.e-7) == 0 &&
      gnrank(nl[:,1:p]*sys[:,m1+1:m1+m2]-nl[:,p+1:p+m2],atol= 1.e-7) == 0 &&
      all(isinf.(gpole(nl,atol=1.e-7))) &&
      info.nrank == 2 && info.degs == [1, 2, 3] && info.stdim == [1, 2, 3]

# minimal inner rational basis
@time nl, info = glnull(sys,m2,atol=1.e-7,coinner = true); info
@test gnrank(nl[:,1:p]*sys[:,1:m1],atol=1.e-7) == 0 &&
      gnrank(nl[:,1:p]*sys[:,m1+1:m1+m2]-nl[:,p+1:p+m2],atol= 1.e-7) == 0 &&
      gnrank(nl[:,1:p]*nl[:,1:p]'-I,atol= 1.e-7) == 0 &&
      info.nrank == 2 && info.degs == [1, 2, 3] && info.stdim == [1, 2, 3]

z = Polynomial([0, 1],'z');
p = 5; m1 = 2; m2 = 0;
G = [ 1 2; -2 z; 2 2*z; -2*z^2 z^3; -2+z^3 z+z^3];
sys = dss(G,Ts = 1);

# minimal inner rational basis (discrete-time)
@time nl, info = glnull(sys,m2,atol=1.e-7,coinner = true); info
@test gnrank(nl[:,1:p]*sys[:,1:m1],atol=1.e-7) == 0 &&
      gnrank(nl[:,1:p]*sys[:,m1+1:m1+m2]-nl[:,p+1:p+m2],atol= 1.e-7) == 0 &&
      gnrank(nl[:,1:p]*nl[:,1:p]'-I,atol= 1.e-7) == 0 &&
      info.nrank == 2 && info.degs == [1, 2, 3] && info.stdim == [1, 2, 3]

z = Polynomial([0, 1],'z');
p = 5; m1 = 2; m2 = 0;
G = [ 1 2; -2 z; 2 2*z; -2*z^2 z^3; -2+z^3 z+z^3];
sys = dss(G,Ts = 1);

# minimal rational simple basis (discrete-time)
@time nl, info = glnull(sys,m2,atol=1.e-7,simple = true); info
@test gnrank(nl[:,1:p]*sys[:,1:m1],atol=1.e-7) == 0 &&
      gnrank(nl[:,1:p]*sys[:,m1+1:m1+m2]-nl[:,p+1:p+m2],atol= 1.e-7) == 0 &&
      info.nrank == 2 && info.degs == [1, 2, 3] && info.stdim == [3, 2, 1]

z = Polynomial([0, 1],'z');
p = 5; m1 = 2; m2 = 0;
G = [ 1 2; -2 z; 2 2*z; -2*z^2 z^3; -2+z^3 z+z^3];
sys = dss(G,Ts = 1);

# minimal rational simple basis with innner vectors (discrete-time)
@time nl, info = glnull(sys,atol=1.e-7, coinner = true, simple = true); info
@test gnrank(nl[:,1:p]*sys[:,1:m1],atol=1.e-7) == 0 &&
      gnrank(nl[:,1:p]*sys[:,m1+1:m1+m2]-nl[:,p+1:p+m2],atol= 1.e-7) == 0 &&
      gnrank(nl[1,:]*nl[1,:]'-I,atol= 1.e-7) == 0 &&
      gnrank(nl[1,:]*nl[1,:]'-I,atol= 1.e-7) == 0 &&
      gnrank(nl[1,:]*nl[1,:]'-I,atol= 1.e-7) == 0 &&
      info.nrank == 2 && info.degs == [1, 2, 3] && info.stdim == [3, 2, 1]

z = Polynomial([0, 1],'z');
p = 5; m1 = 2; m2 = 0;
G = [ 1 2; -2 z; 2 2*z; -2*z^2 z^3; -2+z^3 z+z^3];
sys = dss(G,Ts = 1);

# minimal polynomial simple basis (discrete-time)
@time nl, info = glnull(sys,atol=1.e-7, polynomial = true, simple = false); info
@test gnrank(nl[:,1:p]*sys[:,1:m1],atol=1.e-7) == 0 &&
      gnrank(nl[:,1:p]*sys[:,m1+1:m1+m2]-nl[:,p+1:p+m2],atol= 1.e-7) == 0 &&
      all(isinf.(gpole(nl,atol=1.e-7))) &&
      info.nrank == 2 && info.degs == [1, 2, 3] && info.stdim == [1, 2, 3]

# minimal polynomial simple basis (discrete-time)
@time nl, info = glnull(sys,atol=1.e-7, polynomial = true, simple = true); info
@test gnrank(nl[:,1:p]*sys[:,1:m1],atol=1.e-7) == 0 &&
      gnrank(nl[:,1:p]*sys[:,m1+1:m1+m2]-nl[:,p+1:p+m2],atol= 1.e-7) == 0 &&
      all(isinf.(gpole(nl,atol=1.e-7))) &&
      info.nrank == 2 && info.degs == [1, 2, 3] && info.stdim == [4, 3, 2]

z = Polynomial([0, 1],'z');
p = 5; m1 = 2; m2 = 0;
G = [ 1 2; -2 z; 2 2*z; -2*z^2 z^3; -2+z^3 z+z^3];
sys = dss(G,Ts = 1);

# minimal simple basis with poles assigned to 0 (discrete-time)
@time nl, info = glnull(sys,atol=1.e-7, poles = [0,0,0], simple = true); info
@test gnrank(nl[:,1:p]*sys[:,1:m1],atol=1.e-7) == 0 &&
      gnrank(nl[:,1:p]*sys[:,m1+1:m1+m2]-nl[:,p+1:p+m2],atol= 1.e-7) == 0 &&
      all(abs.(gpole(nl,atol=1.e-7)).< 0.001) &&
      info.nrank == 2 && info.degs == [1, 2, 3] && info.stdim == [3, 2, 1]

z = Polynomial([0, 1],'z');
p = 5; m1 = 2; m2 = 0;
G = [ 1 2; -2 z; 2 2*z; -2*z^2 z^3; -2+z^3 z+z^3];
sys = dss(G,Ts = 1);

# minimal rational basis with poles assigned to 0 (discrete-time)
@time nl, info = glnull(sys,atol=1.e-7, sdeg = 0, simple = false); info
@test gnrank(nl[:,1:p]*sys[:,1:m1],atol=1.e-7) == 0 &&
      gnrank(nl[:,1:p]*sys[:,m1+1:m1+m2]-nl[:,p+1:p+m2],atol= 1.e-7) == 0 &&
      all(abs.(gpole(nl,atol=1.e-7)).< 0.01) &&
      info.nrank == 2 && info.degs == [1, 2, 3] && info.stdim == [1, 2, 3]

# Kailath 1980, pag. 459, rank 2 matrix with both left and right nullspaces
s = Polynomial([0, 1],'s');
p = 3; m1 = 4; m2 = 0;
gn = [1 0 1 s;
    0 (s+1)^2 (s+1)^2 0;
    -1 (s+1)^2 s^2+2*s -s^2];
gd = [s 1 s 1;
    1 1 1 1;
    1 1 1 1];
sys = dss(gn,gd,minimal = true, atol = 1.e-7);

# minimal rational basis with poles assigned to -1
@time nl, info = glnull(sys,atol=1.e-7,poles = [-1]); info
@test gnrank(nl[:,1:p]*sys[:,1:m1],atol=1.e-7) == 0 &&
      gnrank(nl[:,1:p]*sys[:,m1+1:m1+m2]-nl[:,p+1:p+m2],atol= 1.e-7) == 0 &&
      all(real(gpole(nl)) .≈ -1) && all(abs.(imag(gpole(nl))) .< 0.000001) &&
      info.nrank == 2 && info.degs == [1] && info.stdim == [1]

# minimal simple basis with poles assigned to -1
@time nl, info = glnull(sys,atol=1.e-7,poles = [-1], simple = true); info
@test gnrank(nl[:,1:p]*sys[:,1:m1],atol=1.e-7) == 0 &&
      gnrank(nl[:,1:p]*sys[:,m1+1:m1+m2]-nl[:,p+1:p+m2],atol= 1.e-7) == 0 &&
      all(real(gpole(nl)) .≈ -1) && all(abs.(imag(gpole(nl))) .< 0.000001) &&
      info.nrank == 2 && info.degs == [1] && info.stdim == [1]

# minimal simple polynomial basis
@time nl, info = glnull(sys,atol=1.e-7,polynomial = true, simple = true); info
@test gnrank(nl[:,1:p]*sys[:,1:m1],atol=1.e-7) == 0 &&
      gnrank(nl[:,1:p]*sys[:,m1+1:m1+m2]-nl[:,p+1:p+m2],atol= 1.e-7) == 0 &&
      all(isinf.(gpole(nl,atol=1.e-7))) &&
      info.nrank == 2 && info.degs == [1] && info.stdim == [2]

# random examples
p = 5; m1 = 2; m2 = 2; m = m1+m2; n = 5;

Ty = Float64; fast = true; simple = false; polynomial = true;
for Ty in (Float64, Complex{Float64})

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

for fast in (true, false)

for polynomial in (false, true)

#non-simple basis
# continuous, standard
sys = rss(n,p,m,T = Ty, disc=false);
@time nl, info = glnull(sys,m2, atol=1.e-7,polynomial = polynomial, simple = false); info
@test gnrank(nl[:,1:p]*sys[:,1:m1],atol=1.e-7) == 0 &&
      gnrank(nl[:,1:p]*sys[:,m1+1:m1+m2]-nl[:,p+1:p+m2],atol= 1.e-7) == 0 &&
      (polynomial ? all(isinf.(ordeigvals(nl.A,nl.E)[1])) : true) &&
      info.nrank == m1 && info.degs == [1, 2, 2] && info.stdim == [2, 3]

# discrete, descriptor
sys = rdss(n,p,m,T = Ty, disc=true);
@time nl, info = glnull(sys,m2, atol=1.e-7,polynomial = polynomial, simple = false); info
@test gnrank(nl[:,1:p]*sys[:,1:m1],atol=1.e-7) == 0 &&
      gnrank(nl[:,1:p]*sys[:,m1+1:m1+m2]-nl[:,p+1:p+m2],atol= 1.e-7) == 0 &&
      (polynomial ? all(isinf.(ordeigvals(nl.A,nl.E)[1])) : true) &&
      info.nrank == m1 && info.degs == [1, 2, 2] && info.stdim == [2, 3]

# discrete, descriptor
sys = rdss(n,p,m,T = Ty, disc=true);
@time nl, info = glnull(sys,m2, atol=1.e-7,polynomial = polynomial, coinner = true, simple = false); info
@test gnrank(nl[:,1:p]*sys[:,1:m1],atol=1.e-7) == 0 &&
      gnrank(nl[:,1:p]*sys[:,m1+1:m1+m2]-nl[:,p+1:p+m2],atol= 1.e-7) == 0 &&
      (polynomial || gnrank(nl[:,1:p]*nl[:,1:p]'-I,atol= 1.e-7) == 0) &&
      (polynomial ? all(isinf.(ordeigvals(nl.A,nl.E)[1])) : true) &&
      info.nrank == m1 && info.degs == [1, 2, 2] && info.stdim == [2, 3]

# discrete, polynomial
sys = rdss(n,p,m,T = Ty, disc=true, id=[3*ones(Int,1);2*ones(Int,1)]);
@time nl, info = glnull(sys, m2, atol=1.e-7,polynomial = polynomial, simple = false); info
@test gnrank(nl[:,1:p]*sys[:,1:m1],atol=1.e-7) == 0 &&
      gnrank(nl[:,1:p]*sys[:,m1+1:m1+m2]-nl[:,p+1:p+m2],atol= 1.e-7) == 0 &&
      (polynomial ? all(isinf.(ordeigvals(nl.A,nl.E)[1])) : true) &&
      info.nrank == m1 && info.degs == [2, 3, 3] && info.stdim == [2, 3, 3]

# simple basis
# continuous, standard
sys = rss(n,p,m,T = Ty, disc=false);
@time nl, info = glnull(sys,m2, atol=1.e-7,polynomial = polynomial, simple = true); info
@test gnrank(nl[:,1:p]*sys[:,1:m1],atol=1.e-7) == 0 &&
      gnrank(nl[:,1:p]*sys[:,m1+1:m1+m2]-nl[:,p+1:p+m2],atol= 1.e-7) == 0 &&
      (polynomial ? all(isinf.(ordeigvals(nl.A,nl.E)[1])) : true) &&
      info.nrank == m1 && info.degs == [1, 2, 2] && info.stdim == (polynomial ? [3, 3, 2] : [2, 2, 1])

# discrete, descriptor
sys = rdss(n,p,m,T = Ty, disc=true);
@time nl, info = glnull(sys,m2, atol=1.e-7,polynomial = polynomial, simple = true); info
@test gnrank(nl[:,1:p]*sys[:,1:m1],atol=1.e-7) == 0 &&
      gnrank(nl[:,1:p]*sys[:,m1+1:m1+m2]-nl[:,p+1:p+m2],atol= 1.e-7) == 0 &&
      (polynomial ? all(isinf.(ordeigvals(nl.A,nl.E)[1])) : true) &&
      info.nrank == m1 && info.degs == [1, 2, 2] && info.stdim == (polynomial ? [3, 3, 2] : [2, 2, 1])

# discrete, descriptor
sys = rdss(n,p,m,T = Ty, disc=true);
@time nl, info = glnull(sys,m2, atol=1.e-7,polynomial = polynomial, coinner = true, simple = true); info
@test gnrank(nl[:,1:p]*sys[:,1:m1],atol=1.e-7) == 0 &&
      gnrank(nl[:,1:p]*sys[:,m1+1:m1+m2]-nl[:,p+1:p+m2],atol= 1.e-7) == 0 &&
      (polynomial || gnrank(nl[1,1:p]*nl[1,1:p]'-I,atol= 1.e-7) == 0) &&
      (polynomial ? all(isinf.(ordeigvals(nl.A,nl.E)[1])) : true) &&
      info.nrank == m1 && info.degs == [1, 2, 2] && info.stdim == (polynomial ? [3, 3, 2] : [2, 2, 1])

# discrete, polynomial
sys = rdss(5,p,m,T = Ty, disc=true, id=[3*ones(Int,1);2*ones(Int,1)]);
@time nl, info = glnull(sys,m2, atol=1.e-7,polynomial = polynomial, simple = true); info
@test gnrank(nl[:,1:p]*sys[:,1:m1],atol=1.e-7) == 0 &&
      gnrank(nl[:,1:p]*sys[:,m1+1:m1+m2]-nl[:,p+1:p+m2],atol= 1.e-7) == 0 &&
      (polynomial ? all(isinf.(ordeigvals(nl.A,nl.E)[1])) : true) &&
      info.nrank == m1 && info.degs == [2, 3, 3] && info.stdim == (polynomial ? [4, 4, 3] : [3, 3, 2] )

end # polynomial

end # fast

end # Ty

end # glnull

end # nullrange
end # module


text:
using OpenType
using OpenType: curves
using GeometryExperiments

function intensity(curve_points, pixel_per_em)
    @assert length(curve_points) == 3
    res = 0.
    for coord in 1:2
        (x̄₁, x̄₂, x̄₃) = getindex.(curve_points, 3 - coord)
        if maximum(getindex.(curve_points, coord)) * pixel_per_em ≤ -0.5
            continue
        end
        rshift = sum(((i, x̄),) -> x̄ > 0 ? (1 << i) : 0, enumerate((x̄₁, x̄₂, x̄₃)))
        code = (0x2e74 >> rshift) & 0x0003
        if code ≠ 0
            a = x̄₁ - 2x̄₂ + x̄₃
            b = x̄₁ - x̄₂
            c = x̄₁
            if isapprox(a, 0, atol=1e-7)
                t₁ = t₂ = c / 2b
            else
                Δ = b ^ 2 - a * c
                if Δ < 0
                    # in classes C and F, only x̄₂ is of the opposite sign
                    # and there may be no real roots.
                    continue
                end
                δ = sqrt(Δ)
                t₁ = (b - δ) / a
                t₂ = (b + δ) / a
            end
            bezier = BezierCurve()
            if code & 0x0001 == 0x0001
                val = clamp(pixel_per_em * bezier(t₁, curve_points)[coord] + 0.5, 0, 1)
            end
            if code > 0x0001
                val = -clamp(pixel_per_em * bezier(t₂, curve_points)[coord] + 0.5, 0, 1)
            end
            res += val * (coord == 1 ? 1 : -1)
        end
    end
    res
end

function intensity(point, glyph::OpenType.Glyph, units_per_em; font_size=12)
    res = sum(curves(glyph)) do p
        poffset = map(Translation(-point), p)
        intensity(poffset, font_size)
    end
    sqrt(abs(res))
end

function plot_outline(glyph)
    cs = curves(glyph)
    p = plot()
    for (i, curve) in enumerate(cs)
        for (i, point) in enumerate(curve)
            color = i == 1 ? :blue : i == 2 ? :cyan : :green
            scatter!(p, [point[1]], [point[2]], legend=false, color=color)
        end
        points = BezierCurve().(0:0.1:1, Ref(curve))
        curve_color = UInt8[255 - Int(floor(i / length(cs) * 255)), 40, 40]
        plot!(p, first.(points), last.(points), color=string('#', bytes2hex(curve_color)))
    end
    p
end

function render_glyph(font, glyph, font_size)
    step = 0.01
    n = Int(inv(step))
    xs = 0:step:1
    ys = 0:step:1

    grid = map(xs) do x
        map(ys) do y
            Point(x, y)
        end
    end

    grid = hcat(grid...)

    is = map(grid) do p
        try
            intensity(p, glyph, font.head.units_per_em; font_size)
        catch e
            if e isa DomainError
                NaN
            else
                rethrow(e)
            end
        end
    end
    @assert !all(iszero, is)

    p = heatmap(is)
    xticks!(p, 1:n ÷ 10:n, string.(xs[1:n ÷ 10:n]))
    yticks!(p, 1:n ÷ 10:n, string.(ys[1:n ÷ 10:n]))
end

render_glyph(font, char::Char, font_size) = render_glyph(font, font[char], font_size)

using Plots

const BezierCurve = GeometryExperiments.BezierCurve

f = OpenTypeFont(joinpath(dirname(@__DIR__), "assets", "fonts", "juliamono-regular.ttf"))

glyph = f.glyphs[64]
plot_outline(glyph)
render_glyph(f, glyph, 12)

glyph = f.glyphs[75]
plot_outline(glyph)
render_glyph(f, glyph, 12)

glyph = f.glyphs[13]
plot_outline(glyph)
render_glyph(f, glyph, 12)

render_glyph(f, '€', 12)

render_glyph(f, 'A', 12)

ext: jl


text:
# These should be moved to GeometryTypes, once figured out the future of Meshes
function uvmesh(prim::GeometryPrimitive, resolution = (2, 2))
    uv = decompose(UV{Float32}, prim, resolution)
    positions = decompose(Point2f0, prim, resolution)
    faces = decompose(GLTriangle, prim, resolution)
    vertices = map(identity, zip(positions, uv))
    view(vertices, faces)
end
function normalmesh(prim)
    V = GeometryTypes.vertices(prim)
    N = GeometryTypes.normals(prim)
    F = decompose(GLTriangle, prim)
    verts = map((a, b)-> (a, b), V, N)
    Base.view(verts, F)
end


text:
function asinh(x::DoubleFloat{T}) where {T<:AbstractFloat}
    result = abs(x)
    result = result + sqrt(square(result) + 1.0)
    result = log(result)
    return copysign(result, x)
end

function acosh(x::DoubleFloat{T}) where {T<:AbstractFloat}
    x < 1.0 && throw(DomainError("$x"))
    !isfinite(x) && return nan(typeof(x))
    result = x + sqrt(square(x) - 1.0)
    result = log(result)
    return result
end

function atanh(x::DoubleFloat{T}) where {T<:AbstractFloat}
    abs(x) > 1.0 && throw(DomainError("$x"))
    !isfinite(x) && return nan(typeof(x))
    twox = DoubleFloat{T}(x.hi+x.hi, x.lo+x.lo)
    result = 1.0 + twox / (1.0 - x)
    result = log(result)
    result = DoubleFloat{T}(result.hi*0.5, result.lo*0.5)
    return result
end

function acsch(x::DoubleFloat{T}) where {T<:AbstractFloat}
    iszero(x) && throw(DomainError("$x"))
    !isfinite(x) && return nan(typeof(x))
    result = inv(x)
    result = sinh(result)
    return result
end

function asech(x::DoubleFloat{T}) where {T<:AbstractFloat}
    iszero(x) && return inf(DoubleFloat{T})
    isone(x) && return zero(DoubleFloat{T})
    !isfinite(x) && return nan(typeof(x))
    (x < 0.0 || x > 1.0) && throw(DomainError("$x"))
    result = inv(x)
    result = cosh(result)
    return result
end

function acoth(x::DoubleFloat{T}) where {T<:AbstractFloat}
    isone(x) && return inf(DoubleFloat{T})
    isone(-x) && return -inf(DoubleFloat{T})
    (-1.0 < x < 1.0) && throw(DomainError("$x"))
    !isfinite(x) && return nan(typeof(x))
    result = inv(x)
    result = tanh(result)
    return result
end

#=
function asinh(x::DoubleFloat{T}) where {T<:AbstractFloat}
    signbit(x) && return -asinh(abs(x))

    result = x + sqrt(square(x) + 1.0)
    result = log(result)
    return result
end

function acosh(x::DoubleFloat{T}) where {T<:AbstractFloat}
    x < 1.0 && throw(DomainError("$x"))
    !isfinite(x) && return nan(typeof(x))

    result = x + sqrt(square(x) - 1.0)
    result = log(result)
    return result
end

function atanh(x::DoubleFloat{T}) where {T<:AbstractFloat}
    abs(x) > 1.0 && throw(DomainError("$x"))
    !isfinite(x) && return nan(typeof(x))

    if abs(x.hi) == one(T)
        if x.hi > 0
            z = posinf(typeof(x))
        else
            z = neginf(typeof(x))
        end
    else
        oneplus  = 1.0 + x
        oneminus = 1.0 - x
        z = oneplus / oneminus
        z = log(z)
        z = DoubleFloat{T}(z.hi/2, z.lo/2)
    end

    return z
end

function acsch(x::DoubleFloat{T}) where {T<:AbstractFloat}
    iszero(x) && throw(DomainError("$x"))
    !isfinite(x) && return nan(typeof(x))
    result = inv(x)
    result = sinh(result)
    return result
end

function asech(x::DoubleFloat{T}) where {T<:AbstractFloat}
    iszero(x) && return inf(DoubleFloat{T})
    isone(x) && return zero(DoubleFloat{T})
    !isfinite(x) && return nan(typeof(x))
    (x < 0.0 || x > 1.0) && throw(DomainError("$x"))
    result = inv(x)
    result = cosh(result)
    return result
end

function acoth(x::DoubleFloat{T}) where {T<:AbstractFloat}
    isone(x) && return inf(DoubleFloat{T})
    isone(-x) && return -inf(DoubleFloat{T})
    (-1.0 < x < 1.0) && throw(DomainError("$x"))
    !isfinite(x) && return nan(typeof(x))
    result = inv(x)
    result = tanh(result)
    return result
end
=#


text:
abstract type AbstractNextBasisSolver end

include("nextbasis/levmarq.jl")
include("nextbasis/corrrelu.jl")
include("nextbasis/projgraddesc.jl")

################################################################################
type SolverProxLiner <: AbstractNextBasisSolver
end

################################################################################
type SolverStageWiseGaussNewton <: AbstractNextBasisSolver
end

################################################################################
type SolverGaussNewton <: AbstractNextBasisSolver
end

################################################################################
type SolverStocGradientDescent <: AbstractNextBasisSolver
end

################################################################################
type SolverNLOpt <: AbstractNextBasisSolver
end

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


text:
#=

  Set covering problem in Julia ConstraintSolver.jl

  Example 9.1-2, page 354ff, from Taha "Operations Research - An Introduction"
  Minimize the number of security telephones in street corners on a campus.

  AMPL model: http://taha.ineg.uark.edu/setcover.txt

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function set_covering2(problem,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   # "all_solutions"=> all_solutions,
                                                            "all_optimal_solutions"=>true,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    n = problem[:n]
    corners = problem[:corners]
    num_rows,_ = size(corners)

    # Where to place the telephones
    @variable(model, x[1:n], Bin)
    @variable(model, 1 <= z <= n, Int)

    # All streets must be covered
    for i in 1:num_rows
        @constraint(model,x[corners[i,1]] + x[corners[i,2]] >= 1)
    end
    @constraint(model, sum(x) == z)
    @objective(model,Min,z)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                z_val = convert.(Integer,JuMP.value.(z; result=sol))
                println("z:$z_val")
                println("x:$x_val")

            end
        end
    else
        println("status:$status")
    end

    return status
end

problem = Dict(
   :n => 8,
   :corners => resize_matrix([[1,2],
                              [2,3],
                              [4,5],
                              [7,8],
                              [6,7],
                              [2,6],
                              [1,6],
                              [4,7],
                              [2,4],
                              [5,8],
                              [3,5]]),

)
@time set_covering2(problem)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
using BenchmarkTools
using InteractiveUtils
using Printf

if length(ARGS) != 2
    println("Usage: $PROGRAM_FILE <filename> <count>")
	exit()
end

versioninfo()
println()

load_time = @elapsed using SASLib
@printf "Loaded library in %.3f seconds\n" load_time

b = @benchmark readsas($ARGS[1], verbose_level=0) samples=parse(Int, ARGS[2])
display(b)

println()
println()


text:
using Revise
using JuMP, EAGO

m = Model(optimizer_with_attributes(EAGO.Optimizer, "verbosity" => 1,
                                                    "output_iterations" => 10000,
                                                    "cp_repetitions" => -1,
                                                    "iteration_limit" => 1000000,
                                                    "cut_min_iterations" => 1,
                                                    "cut_max_iterations" => 1,
                                                    "objective_cut_on" => true,
                                                    "subgrad_tighten" => true,
                                                    "obbt_depth" => -1,
                                                    "fbbt_lp_depth" => -1))

x_Idx = Any[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24]

@variable(m, 0 <= x[x_Idx] <= 100)

# ----- Constraints ----- #
@constraint(m, e2, x[1]+x[2]+x[3]+x[4] == 8.0)
@constraint(m, e3, x[5]+x[6]+x[7]+x[8] == 24.0)
@constraint(m, e4, x[9]+x[10]+x[11]+x[12] == 20.0)
@constraint(m, e5, x[13]+x[14]+x[15]+x[16] == 24.0)
@constraint(m, e6, x[17]+x[18]+x[19]+x[20] == 16.0)
@constraint(m, e7, x[21]+x[22]+x[23]+x[24] == 12.0)
@constraint(m, e8, x[1]+x[5]+x[9]+x[13]+x[17]+x[21] == 29.0)
@constraint(m, e9, x[2]+x[6]+x[10]+x[14]+x[18]+x[22] == 41.0)
@constraint(m, e10, x[3]+x[7]+x[11]+x[15]+x[19]+x[23] == 13.0)
@constraint(m, e11, x[4]+x[8]+x[12]+x[16]+x[20]+x[24] == 21.0)

@NLobjective(m, Min, (300*x[1]-7*x[1]*x[1]-4*x[2]*x[2]+270*x[2]-6*x[3]*x[3]+460*x[3]-8*x[4]*x[4]+800*x[4]-12*x[5]*x[5]+740*x[5]-9*x[6]*x[6]+600*x[6]-14*x[7]*x[7]+540*x[7]-7*x[8]*x[8]+380*x[8]-13*x[9]*x[9]+300*x[9]-12*x[10]*x[10]+490*x[10]-8*x[11]*x[11]+380*x[11]-4*x[12]*x[12]+760*x[12]-7*x[13]*x[13]+430*x[13]-9*x[14]*x[14]+250*x[14]-16*x[15]*x[15]+390*x[15]-8*x[16]*x[16]+600*x[16]-4*x[17]*x[17]+210*x[17]-10*x[18]*x[18]+830*x[18]-21*x[19]*x[19]+470*x[19]-13*x[20]*x[20]+680*x[20]-17*x[21]*x[21]+360*x[21]-9*x[22]*x[22]+290*x[22]-8*x[23]*x[23]+400*x[23]-4*x[24]*x[24]+310*x[24]))

JuMP.optimize!(m)
run_time = backend(m).optimizer.model.optimizer._run_time
#objective_value =
#objective_bound =
#println("objective between: ")
println("run time: $run_time")


text:
######################################################################
# tsp.jl
#
# Some code for the traveling salesman problem
######################################################################

# uses nearest neighbors to heuristically solve tsp
# tsp is traveling salesman problem
function tsp_nn!(xd::VVF, points::VVF)
	xc = xd[1]
	next_p = xd[1]
	n = 1
	while length(points) > 0
		best_d = Inf
		best_ind = 1
		for (p_ind,p) in enumerate(points)
			r = xc - p
			d = dot(r, r)
			if d < best_d
				best_d = d
				best_ind = p_ind
			end
		end
		xd[n+1] = points[best_ind]
		deleteat!(points, best_ind)
		n += 1
	end
end

# simply takes the points as they are
# this creates a shitty trajectory
function tsp_rand!(xd::VVF, points::VVF)
	n = 1
	for p in points
		xd[n+1] = p
		n += 1
	end
end


text:
## subarrays ##

typealias RangeIndex Union(Int, Range{Int}, UnitRange{Int})

type SubArray{T,N,A<:AbstractArray,I<:(RangeIndex...,)} <: AbstractArray{T,N}
    parent::A
    indexes::I
    dims::NTuple{N,Int}
    strides::Array{Int,1}  # for accessing parent with linear indexes
    first_index::Int

    # Note: no bounds-checking on construction. See issue #4044
    #linear indexing constructor (scalar)
    if N == 0 && length(I) == 1 && A <: Array
        function SubArray(p::A, i::(Int,))
            new(p, i, (), Int[], i[1])
        end
    #linear indexing constructor (ranges)
    elseif N == 1 && length(I) == 1 && A <: Array
        function SubArray(p::A, i::(UnitRange{Int},))
            new(p, i, (length(i[1]),), [1], first(i[1]))
        end
        function SubArray(p::A, i::(Range{Int},))
            new(p, i, (length(i[1]),), [step(i[1])], first(i[1]))
        end
    else
        function SubArray(p::A, i::I)
            newdims = Array(Int, 0)
            newstrides = Array(Int, 0)
            newfirst = 1
            pstride = 1
            for j = 1:length(i)
                if isa(i[j], Int)
                    newfirst += (i[j]-1)*pstride
                else
                    push!(newdims, length(i[j]))
                    #may want to return error if step(i[j]) <= 0
                    push!(newstrides, isa(i[j],UnitRange) ? pstride :
                         pstride * step(i[j]))
                    newfirst += (first(i[j])-1)*pstride
                end
                pstride *= size(p,j)
            end
            new(p, i, tuple(newdims...), newstrides, newfirst)
        end
    end
end

#linear indexing sub (may want to rename as slice)
function sub{T,N}(A::Array{T,N}, i::(Union(Range{Int}, UnitRange{Int}),))
    SubArray{T,1,typeof(A),typeof(i)}(A, i)
end

# if `I` were a vector, index_ranges would do the following:
# j = length(I)
# while j > 0 && isa(I[j], Int)
#     j -= 1
# end
# for i = 1:j
#     if isa(I[i], Int)
#         I[i] = I[i]:I[i]
#     end
# end
to_range(j::Int) = j:j
to_range(j::RangeIndex) = j
index_ranges(I::Int...) = I
index_ranges(i, I...) = tuple(to_range(i), index_ranges(I...)...)

function sub_internal{T,N,L}(A::AbstractArray{T,N}, i::NTuple{N,RangeIndex}, ::NTuple{L,Int})
    SubArray{T,L,typeof(A),typeof(i)}(A, i)
end

function sub{T,N}(A::AbstractArray{T,N}, i::NTuple{N,RangeIndex})
    sub_internal(A, index_ranges(i...), index_shape(i...))
end

sub{N}(A::SubArray, i::NTuple{N,RangeIndex}) = sub(A, i...)

sub(A::AbstractArray, i::RangeIndex...) = sub(A, i)

function sub(A::SubArray, i::RangeIndex...)
    L = length(i)
    while L > 0 && isa(i[L], Int); L-=1; end
    j = 1
    newindexes = Array(RangeIndex,length(A.indexes))
    for k = 1:length(A.indexes)
        if isa(A.indexes[k], Int)
            newindexes[k] = A.indexes[k]
        else
            r = A.indexes[k]
            ri = (isa(i[j],Int) && j<=L) ? (i[j]:i[j]) : i[j]
            newindexes[k] = step(r) == 1 ? (first(r)-1) + ri : first(r) + (ri-1)*step(r)
            j += 1
        end
    end
    ni = tuple(newindexes...)
    SubArray{eltype(A),L,typeof(A.parent),typeof(ni)}(A.parent, ni)
end

# Drops all Ints from a tuple of RangeIndexes
ranges_only(I::Int...) = ()
ranges_only(i::Int, I...) = ranges_only(I...)
ranges_only(i::Union(Range{Int}, UnitRange{Int}), I...) = tuple(i, ranges_only(I...)...)

function slice_internal{T,N,L}(A::AbstractArray{T,N}, i::NTuple{N,RangeIndex}, ::NTuple{L,RangeIndex})
    SubArray{T,L,typeof(A),typeof(i)}(A, i)
end
slice{T,N}(A::AbstractArray{T,N}, i::NTuple{N,RangeIndex}) = slice_internal(A, i, ranges_only(i...))

# Throw error on slice dimension mismatch
slice{T,N,M}(A::AbstractArray{T,N}, i::NTuple{M,RangeIndex}) = throw(BoundsError())

slice(A::AbstractArray, i::RangeIndex...) = slice(A, i)

function slice(A::SubArray, i::RangeIndex...)
    j = 1
    newindexes = Array(RangeIndex,length(A.indexes))
    for k = 1:length(A.indexes)
        if isa(A.indexes[k], Int)
            newindexes[k] = A.indexes[k]
        else
            r = A.indexes[k]
            newindexes[k] = step(r) == 1 ? (first(r)-1) + i[j] : first(r) + (i[j]-1)*step(r)
            j += 1
        end
    end
    slice(A.parent, tuple(newindexes...))
end

# Colon translation
sub(A::AbstractArray, I::Union(RangeIndex, Colon)...) = sub(A, ntuple(length(I), i-> isa(I[i], Colon) ? (1:size(A,i)) : I[i])...)
slice(A::AbstractArray, I::Union(RangeIndex, Colon)...) = slice(A, ntuple(length(I), i-> isa(I[i], Colon) ? (1:size(A,i)) : I[i])...)

### rename the old slice function ###
##squeeze all dimensions of length 1
#slice{T,N}(a::AbstractArray{T,N}) = sub(a, map(i-> i == 1 ? 1 : (1:i), size(a)))
#slice{T,N}(s::SubArray{T,N}) =
#    sub(s.parent, map(i->!isa(i, Int) && length(i)==1 ?i[1] : i, s.indexes))
#
##slice dimensions listed, error if any have length > 1
##silently ignores dimensions that are greater than N
#function slice{T,N}(a::AbstractArray{T,N}, sdims::Integer...)
#    newdims = ()
#    for i = 1:N
#        next = 1:size(a, i)
#        for j in sdims
#            if i == j
#                if size(a, i) != 1
#                    error("dimension ", i, " has length greater than 1")
#                end
#                next = 1
#                break
#            end
#        end
#        newdims = tuple(newdims..., next)
#    end
#    sub(a, newdims)
#end
#function slice{T,N}(s::SubArray{T,N}, sdims::Integer...)
#    newdims = ()
#    for i = 1:length(s.indexes)
#        next = s.indexes[i]
#        for j in sdims
#            if i == j
#                if length(next) != 1
#                    error("dimension ", i," has length greater than 1")
#                end
#                next = isa(next, Int) ? next : first(next)
#                break
#            end
#        end
#        newdims = tuple(newdims..., next)
#    end
#    sub(s.parent, newdims)
#end
### end commented code ###

size(s::SubArray) = s.dims
ndims{T,N}(s::SubArray{T,N}) = N

parent(s::SubArray) = s.parent
parentindexes(s::SubArray) = s.indexes

parent(a::AbstractArray) = a
parentindexes(a::AbstractArray) = ntuple(ndims(a), i->1:size(a,i))

copy(s::SubArray) = copy!(similar(s.parent, size(s)), s)
similar(s::SubArray, T, dims::Dims) = similar(s.parent, T, dims)

getindex{T}(s::SubArray{T,0}) = s.parent[s.first_index]

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

getindex(s::SubArray, i::Real) = getindex(s, to_index(i))
getindex(s::SubArray, i0::Real, i1::Real) =
    getindex(s, to_index(i0), to_index(i1))
getindex(s::SubArray, i0::Real, i1::Real, i2::Real) =
    getindex(s, to_index(i0), to_index(i1), to_index(i2))
getindex(s::SubArray, i0::Real, i1::Real, i2::Real, i3::Real) =
    getindex(s, to_index(i0), to_index(i1), to_index(i2), to_index(i3))
getindex(s::SubArray, i0::Real, i1::Real, i2::Real, i3::Real, i4::Real) =
    getindex(s, to_index(i0), to_index(i1), to_index(i2), to_index(i3), to_index(i4))
getindex(s::SubArray, i0::Real, i1::Real, i2::Real, i3::Real, i4::Real, i5::Real) =
    getindex(s, to_index(i0), to_index(i1), to_index(i2), to_index(i3), to_index(i4), to_index(i5))
getindex(s::SubArray, i0::Real, i1::Real, i2::Real, i3::Real, i4::Real, i5::Real, is::Real...) =
    getindex(s, to_index(i0), to_index(i1), to_index(i2), to_index(i3), to_index(i4), to_index(i5), to_index(is)...)

getindex(s::SubArray, i::Integer) = s[ind2sub(size(s), i)...]

function getindex(s::SubArray, is::Integer...)
    index = s.first_index
    for i = 1:length(is)
        isi = is[i]
        if isi != 1
            index += (is[i]-1)*s.strides[i]
        end
    end
    s.parent[index]
end

function getindex_bool_1d(S::SubArray, I::AbstractArray{Bool})
    n = sum(I)
    out = similar(S, n)
    c = 1
    for i = 1:length(I)
        if I[i]
            out[c] = S[i]
            c += 1
        end
    end
    out
end

getindex{T}(S::SubArray{T,1}, I::AbstractArray{Bool,1}) = getindex_bool_1d(S, I)
getindex{T}(S::SubArray{T,2}, I::AbstractArray{Bool,2}) = getindex_bool_1d(S, I)
getindex{T}(S::SubArray{T,3}, I::AbstractArray{Bool,3}) = getindex_bool_1d(S, I)
getindex{T}(S::SubArray{T,4}, I::AbstractArray{Bool,4}) = getindex_bool_1d(S, I)
getindex{T}(S::SubArray{T,5}, I::AbstractArray{Bool,5}) = getindex_bool_1d(S, I)

getindex{T}(s::SubArray{T,1}, I::UnitRange{Int}) =
    getindex(s.parent, (s.first_index+(first(I)-1)*s.strides[1]):s.strides[1]:(s.first_index+(last(I)-1)*s.strides[1]))

getindex{T}(s::SubArray{T,1}, I::Range{Int}) =
    getindex(s.parent, (s.first_index+(first(I)-1)*s.strides[1]):(s.strides[1]*step(I)):(s.first_index+(last(I)-1)*s.strides[1]))

function getindex{T,S<:Integer}(s::SubArray{T,1}, I::AbstractVector{S})
    t = Array(Int, length(I))
    for i = 1:length(I)
        t[i] = s.first_index + (I[i]-1)*s.strides[1]
    end
    getindex(s.parent, t)
end

function translate_indexes(s::SubArray, I::Union(Real,AbstractArray)...)
    n = length(I)
    newindexes = Any[s.indexes...]
    pdims = parentdims(s)
    havelinear = n < ndims(s)
    for i = 1:n-havelinear
        newindexes[pdims[i]] = s.indexes[pdims[i]][I[i]]
    end
    if havelinear
        newindexes = newindexes[1:pdims[n]]
        newindexes[pdims[n]] = translate_linear_indexes(s, n, I[end], pdims)
    end
    newindexes
end

# translate a linear index vector I for dim n to a linear index vector for
# the parent array
function translate_linear_indexes(s, n, I, pdims)
    idx = Array(Int, length(I))
    ssztail = size(s)[n:end]
    indexestail = s.indexes[pdims[n:end]]
    # The next gets the strides of dimensions listed in pdims[n:end], relative to the stride of pdims[n]
    pstrd = [1]
    j = n+1
    strd = 1
    for i = pdims[n]+1:ndims(s.parent)
        strd *= size(s.parent, i-1)
        if j <= length(pdims) && i == pdims[j]
            push!(pstrd, strd)
            j += 1
        end
    end
    # Compute the offset from any omitted dimensions
    taildimsoffset = 0
    for i = pdims[n]+1:ndims(s.parent)
        thisI = s.indexes[i]
        if isa(thisI, Integer)
            taildimsoffset += (thisI-1)*stride(s.parent, i)
        end
    end
    nd = length(pstrd)
    for j=1:length(I)
        su = ind2sub(ssztail,I[j])  # convert to particular location within indexes
        K = taildimsoffset + 1
        for k = 1:nd
            K += pstrd[k]*(indexestail[k][su[k]]-1)   # convert to particular location in parent
        end
        idx[j] = K
    end
    idx
end

function parentdims(s::SubArray)
    nd = ndims(s)
    dimindex = Array(Int, nd)
    sp = strides(s.parent)
    j = 1
    for i = 1:ndims(s.parent)
        r = s.indexes[i]
        if j <= nd && (isa(r,Range) ? sp[i]*step(r) : sp[i]) == s.strides[j]
            dimindex[j] = i
            j += 1
        end
    end
    dimindex
end

function getindex(s::SubArray, I::Union(Real,AbstractVector)...)
    newindexes = translate_indexes(s, I...)

    rs = index_shape(I...)
    result = getindex(s.parent, newindexes...)
    if isequal(rs, size(result))
        return result
    else
        return reshape(result, rs)
    end
end

setindex!(s::SubArray, v, i::Integer) = setindex!(s, v, ind2sub(size(s), i)...)

function setindex!(s::SubArray, v, is::Integer...)
    index = s.first_index
    for i = 1:length(is)
        index += (is[i]-1)*s.strides[i]
    end
    s.parent[index] = v
    return s
end

setindex!{T}(s::SubArray{T,0}, v) = setindex!(s.parent, v, s.first_index)

setindex!{T}(s::SubArray{T,1}, v, i::Integer) =
    setindex!(s.parent, v, s.first_index + (i-1)*s.strides[1])

setindex!{T}(s::SubArray{T,2}, v, i::Integer, j::Integer) =
    setindex!(s.parent, v, s.first_index +(i-1)*s.strides[1]+(j-1)*s.strides[2])

setindex!{T}(s::SubArray{T,3}, v, i::Integer, j::Integer, k::Integer) =
    setindex!(s.parent, v, s.first_index +(i-1)*s.strides[1]+(j-1)*s.strides[2]+(k-1)*s.strides[3])

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

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

setindex!{T}(s::SubArray{T,1}, v, I::UnitRange{Int}) =
    setindex!(s.parent, v, (s.first_index+(first(I)-1)*s.strides[1]):s.strides[1]:(s.first_index+(last(I)-1)*s.strides[1]))

setindex!{T}(s::SubArray{T,1}, v, I::Range{Int}) =
    setindex!(s.parent, v, (s.first_index+(first(I)-1)*s.strides[1]):(s.strides[1]*step(I)):(s.first_index+(last(I)-1)*s.strides[1]))

function setindex!{T,S<:Integer}(s::SubArray{T,1}, v, I::AbstractVector{S})
    t = Array(Int, length(I))
    for i = 1:length(I)
        t[i] = s.first_index + (I[i]-1)*s.strides[1]
    end
    setindex!(s.parent, v, t)
end

# to avoid ambiguity warning
function setindex!(s::SubArray, v, I::Real)
    newindexes = translate_indexes(s, (to_index(I),))
    setindex!(s.parent, v, newindexes...)
end
function setindex!(s::SubArray, v, I::Union(Real,AbstractArray)...)
    newindexes = translate_indexes(s, to_index(I)...)
    setindex!(s.parent, v, newindexes...)
end

stride(s::SubArray, i::Integer) = i <= length(s.strides) ? s.strides[i] : s.strides[end]*s.dims[end]

convert{T}(::Type{Ptr{T}}, x::SubArray{T}) =
    pointer(x.parent) + (x.first_index-1)*sizeof(T)
convert{T,S,N}(::Type{Array{T,N}},A::SubArray{S,N}) = copy!(Array(T,size(A)), A)

pointer(s::SubArray, i::Int) = pointer(s, ind2sub(size(s), i))

function pointer(s::SubArray, is::(Int...))
    index = s.first_index
    for n = 1:length(is)
        index += (is[n]-1)*s.strides[n]
    end
    return pointer(s.parent, index)
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
__precompile__()

module DiscreteLog

################################################################################
#
#  Import/export
#
################################################################################

using LinearAlgebra, Distributed, SparseArrays, Serialization, Random, Pkg, Test

import AbstractAlgebra

import SparseArrays

import Serialization: serialize, deserialize

###############################################################################
#
#   Library initialization
#
###############################################################################

const pkgdir = joinpath(dirname(pathof(DiscreteLog)), "..")

function __init__()

    if myid() == 1
        println("")
        print("Welcome to DiscreteLog.jl")
        println()
    end
end

################################################################################
#
#  Submodules
#
################################################################################

include("../src/ElimProcedures/Deg32Elim.jl")
include("../src/ElimProcedures/Deg43Elim.jl")

include("../src/EllCrv/EllCrv.jl")
include("../src/EllCrv/EllCrvDiv.jl")
include("../src/EllCrv/EllCrvFin.jl")
include("../src/EllCrv/EllCrvMap.jl")

include("../src/Misc/BigIntInterface.jl")
include("../src/Misc/FieldsRings.jl")
include("../src/Misc/PolyFac.jl")

end


text:
module atn_module

using Statistics

function infected_c(
    AgentState,
    AgentAct,
    m,
    mu,
    lambda,
    etaS,
    eta,
    InfectedOverATrial,
    tmax,
    TWindow,
)
    N = length(AgentAct)
    AgentStateNext = zeros(Bool, N)
    # AgentState = zeros(Bool, N)

    # AgentState = zeros(Bool, N)
    # AgentState[1:Ni] .= 1
    for k = 1:N
        AgentStateNext[k] = AgentState[k]
    end
    for i = 1:N
        if AgentState[i] == 1
            if rand() < mu
                AgentStateNext[i] = 0
            end
        end
        # Activity =
        #     (1 - AgentState[i]) * (etaS * AgentAct[i]) +
        #     AgentState[i] * (eta * AgentAct[i])
        Activity =
            (1 - AgentState[i]) * (eta * AgentAct[i]) +
            AgentState[i] * (eta * AgentAct[i])

        if Activity > rand()
            for j = 1:m
                rInt = rand(1:N)
                if AgentState[i] == 1
                    r = rand()
                    if ((AgentState[rInt] == 0) && (rand() < lambda))
                        AgentStateNext[rInt] = 1
                    end
                else
                    if ((AgentState[rInt] == 1) && (rand() < lambda))
                        AgentStateNext[i] = 1
                    end
                end
            end
        end
    end
    for k = 1:N
        AgentState[k] = AgentStateNext[k]
        # if AgentState[k] == 1
        #     NumInf = NumInf + 1
        # end
    end

    return AgentState

end

function avg_infected(
    N,
    NTrials,
    m,
    mu,
    eta,
    lambda,
    etaS,
    AgentAct,
    tmax,
    TWindow,
)
    InfectedOverTrials = 0.0
    for nt = 1:NTrials
        AgentState = zeros(Bool, N)
        AgentState[1:Int(0.01 * N)] .= 1
        InfectedOverATrial = 0.0
        for i_time = 1:tmax
            NumInf, AgentState = infected_c(
                AgentState,
                AgentAct,
                m,
                mu,
                lambda,
                etaS,
                eta,
                InfectedOverATrial,
                tmax,
                TWindow,
            )
            if i_time > (tmax - TWindow)
                NumInf = sum(AgentState)
                InfectedOverATrial = InfectedOverATrial + NumInf
            end
        end
        InfectedOverATrial = InfectedOverATrial / TWindow
        InfectedOverTrials = InfectedOverTrials + InfectedOverATrial
    end
    return (InfectedOverTrials / NTrials) / N
end

function avg_infected_par(
    N,
    NTrials,
    m,
    mu,
    grid,
    etaS,
    AgentAct,
    tmax,
    TWindow,
)
    eta = grid[1]
    lambda = grid[2]
    InfectedOverTrials = 0.0
    for nt = 1:NTrials
        AgentState = zeros(Bool, N)
        for k = 1:100
            AgentState[k] = 1
        end
        InfectedOverATrial = 0.0
        for i_time = 1:tmax
            AgentState = infected_c(
                AgentState,
                AgentAct,
                m,
                mu,
                lambda,
                etaS,
                eta,
                InfectedOverATrial,
                tmax,
                TWindow,
            )
            if i_time > (tmax - TWindow)
                NumInf = sum(AgentState)
                InfectedOverATrial = InfectedOverATrial + NumInf
            end
        end
        InfectedOverATrial = InfectedOverATrial / TWindow
        InfectedOverTrials = InfectedOverTrials + InfectedOverATrial
    end
    return (InfectedOverTrials / NTrials) / N
end

end


text:
function update_tau!(state::State, data::Data, prior::Prior)
  a, b = params(prior.tau)
  new_shape = a + prior.K * prior.a_omega
  new_rate = b + sum(1 ./ state.omega)
  state.tau = rand(Gamma(new_shape, 1 / new_rate))
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
module ISICArchive

import Requests: get
using JSON
using Parameters
using Images
using ImageMagick

export

    DatasetListRequest,
    DatasetMetadataRequest,

    ImageListRequest,
    ImageMetadataRequest,
    ImageDownloadRequest,
    ImageThumbnailRequest

include("common.jl")
include("dataset.jl")
include("image.jl")

end # module


text:
### Zero variance (ZV) Monte Carlo estimators
### Reference:
### Mira A, Solgi R and Imparato D
### Zero Variance Markov Chain Monte Carlo for Bayesian Estimators
### Statistics and Computing, 2013, 23 (5), pp 653-662

## Zero variance (ZV) Monte Carlo estimators linear polynomials

function lzv(chain::Vector{N}, grad::Vector{N}) where N<:Real
  z = -0.5*grad
  augmentedcov = cov([z chain])
  a = -augmentedcov[1, 2]/augmentedcov[1, 1]
  return chain+z*a, a
end

function lzv(chain::Matrix{N}, grad::Matrix{N}) where N<:Real
  npars = size(chain, 2)

  augmentedcov = Array(N, npars+1, npars+1, npars)
  precision = Array(N, npars, npars, npars)
  sigma = Array(N, npars, npars)
  a = Array(N, npars, npars)

  z = -0.5*grad

  for i in 1:npars
     augmentedcov[:, :, i] = cov([z chain[:, i]])
     precision[:, :, i] = inv(augmentedcov[1:npars, 1:npars, i])
     sigma[:, i] = augmentedcov[1:npars, npars+1, i]
     a[:, i] = -precision[:, :, i]*sigma[:, i]
  end

  return chain+z*a, a
end

lzv(s::ParameterNState{Continuous, Univariate}) = lzv(s.value, s.gradlogtarget)

lzv(s::ParameterNState{Continuous, Multivariate}) = lzv(transpose(s.value), transpose(s.gradlogtarget))

## Zero variance (ZV) Monte Carlo estimators quadratic polynomials

function qzv(chain::Vector{N}, grad::Vector{N}) where N<:Real
  z = -0.5*grad
  qz = [z 2*z.*chain-1]
  augmentedcov = cov([qz chain])
  a = -inv(augmentedcov[1:2, 1:2])*augmentedcov[1:2, 3]
  return chain+qz*a, a
end

function qzv(chain::Matrix{N}, grad::Matrix{N}) where N<:Real
  nsamples, npars = size(chain)
  k = convert(Integer, npars*(npars+3)/2)
  l = 2*npars+1

  qz = Array(N, nsamples, k)
  augmentedcov = Array(N, k+1, k+1, npars)
  precision = Array(N, k, k, npars)
  sigma = Array(N, k, npars)
  a = Array(N, k, npars)

  z = -grad/2

  qz[:, 1:npars] = z
  qz[:, (npars+1):(2*npars)] = 2*z.*chain-1
  for i in 1:(npars-1)
    for j = (i+1):npars
      qz[:, l] = chain[:, i].*z[:, j]+chain[:, j].*z[:, i]
      l += 1
    end
  end

  for i in 1:npars
    augmentedcov[:, :, i] = cov([qz chain[:, i]])
    precision[:, :, i] = inv(augmentedcov[1:k, 1:k, i])
    sigma[:, i] = augmentedcov[1:k, k+1, i]
    a[:, i] = -precision[:, :, i]*sigma[:, i]
  end

  return chain+qz*a, a
end

qzv(s::ParameterNState{Continuous, Univariate}) = qzv(s.value, s.gradlogtarget)

qzv(s::ParameterNState{Continuous, Multivariate}) = qzv(transpose(s.value), transpose(s.gradlogtarget))


text:
# V-measure of contingency table
function _vmeasure(A::AbstractMatrix{<:Integer}; β::Real)
    (β >= 0) || throw(ArgumentError("β should be nonnegative"))

    N = sum(A)
    (N == 0.0) && return 0.0

    entA = entropy(A)
    entArows = entropy(sum(A, dims=2))
    entAcols = entropy(sum(A, dims=1))

    hck = (entA - entAcols)/N
    hkc = (entA - entArows)/N
    hc = entArows/N + log(N)
    hk = entAcols/N + log(N)

    # Homogeneity
    h = hc == 0.0 ? 1.0 : 1.0 - hck/hc
    # Completeness
    c = hk == 0.0 ? 1.0 : 1.0 - hkc/hk

    # V-measure
    V_β = (1 + β)*h*c/(β*h + c)
    return V_β
end

"""
    vmeasure(assign1, assign2; β = 1.0)

V-measure between two clustering assignments.

`assign1` and `assign2` can be either `ClusteringResult` objects or
assignments vectors (`AbstractVector{<:Integer}`).

The `β` parameter defines trade-off between _homogeneity_ and _completeness_:
 * if `β` is greater than 1, _completeness_ is weighted more strongly,
 * if `β` is less than 1, _homogeneity_ is weighted more strongly.

*Ref:* Andrew Rosenberg and Julia Hirschberg, 2007. "V-Measure: A conditional entropy-based external cluster evaluation measure"
"""
function vmeasure(assign1::Union{AbstractVector{<:Integer}, ClusteringResult},
                  assign2::Union{AbstractVector{<:Integer}, ClusteringResult};
                  β::Real = 1.0)
    _assign1 = isa(assign1, AbstractVector) ? assign1 : assignments(assign1)
    _assign2 = isa(assign2, AbstractVector) ? assign2 : assignments(assign2)
    return _vmeasure(counts(_assign1, _assign2,
                            (1:maximum(_assign1), 1:maximum(_assign2))), β=β)
end