text:
module Test_iofac

using DescriptorSystems
using LinearAlgebra
using Polynomials
using Test

println("Test_iofac")
@testset "grsfg and grsfg" begin

sys = rss(0,2,3);
g = glinfnorm(sys)[1]*1.001
@time sysf = grsfg(sys,g)
@test iszero(sysf'*sysf+sys'*sys-g*g*I,atol=1.e-7)

@time sysf = glsfg(sys,g)
@test iszero(sysf*sysf'+sys*sys'-g*g*I,atol=1.e-7)

# Error: Norm(sys) < g
@test_throws ErrorException sysf = grsfg(sys,0.5*g)
@test_throws ErrorException sysf = glsfg(sys,0.5*g)

# Eigenvalue(s) on the imaginary axis
s = rtf('s');
sys = dss(1/s);
@test_throws ErrorException sysf = grsfg(sys,1)
@test_throws ErrorException sysf = glsfg(sys,1)

# Improper continuous-system
s = rtf('s');
sys = dss(s);
@test_throws ErrorException sysf = grsfg(sys,1)
@test_throws ErrorException sysf = glsfg(sys,1)

s = rtf('s');
sys = [(s-1)/(s+2) s/(s+2) 1/(s+2);
    0 (s-2)/(s+1)^2 (s-2)/(s+1)^2;
    (s-1)/(s+2) (s^2+2*s-2)/(s+1)/(s+2) (2*s-1)/(s+1)/(s+2)];
sys = gir(dss(sys),atol=1.e-7);
g = glinfnorm(sys)[1]*1.001

sysf = grsfg(sys,g,stabilize=false)
@test iszero(sysf'*sysf+sys'*sys-g*g*I,atol=1.e-7)
@time sysf = glsfg(sys,g,stabilize=false)
@test iszero(sysf*sysf'+sys*sys'-g*g*I,atol=1.e-7)

z = rtf('z');
sys = [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 = gir(dss(sys),atol=1.e-7);
g = glinfnorm(sys)[1]*1.001

sysf = grsfg(sys,g)
@test iszero(sysf'*sysf+sys'*sys-g*g*I,atol=1.e-7)
@time sysf = glsfg(sys,g)
@test iszero(sysf*sysf'+sys*sys'-g*g*I,atol=1.e-7)

sys = rdss(5,4,3,id = ones(Int,2))
g = glinfnorm(sys)[1]*1.001
sysf = grsfg(sys,g,atol=1.e-7)
@test iszero(sysf'*sysf+sys'*sys-g*g*I,atol=1.e-7)
@time sysf = glsfg(sys,g,atol=1.e-7)
@test iszero(sysf*sysf'+sys*sys'-g*g*I,atol=1.e-7)

end # grsfg & glsfg

@testset "giofac and goifac" begin

@testset "giofac" begin

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

@time sysi, syso, info = giofac(sys)
@test gnrank(sys-sysi*syso) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysi'*sysi-I) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysi) & (isproper(sys) ? isproper(syso) : true) &&
      isstable(sysi) && (isstable(sys) ? isstable(syso) : true) &&
      isempty(gzero(gminreal(syso))) &&
      info.nrank == 0 && info.nfuz == 0 && info.niuz == 0

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

@time sysi, syso, info = giofac(sys)
@test gnrank(sys-sysi*syso) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysi'*sysi-I) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysi) & (isproper(sys) ? isproper(syso) : true) &&
      isstable(sysi) && (isstable(sys) ? isstable(syso) : true) &&
      isempty(gzero(gminreal(syso))) &&
      info.nrank == 0 && info.nfuz == 0 && 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 sysi, syso, info = giofac(sys, minphase = false, atol1 = 1.e-7, atol2 = 1.e-7)
@test gnrank(sys-sysi*syso,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysi'*sysi-I) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysi) & isproper(syso) && # checking properness of factors
      isstable(sysi) && (isstable(sys) ? isstable(syso) : true) &&
      count(t -> isinf(t), gzero(gminreal(syso))) == 1 &&
      info.nrank == 1 && ismissing(info.nfuz) && info.niuz == 1

@time sysi, syso, info = giofac(sys, minphase = true, atol1 = 1.e-7, atol2 = 1.e-7)
@test gnrank(sys-sysi*syso,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysi'*sysi-I) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysi) & isproper(syso) && # checking properness of factors
      isstable(sysi) && (isstable(sys) ? isstable(syso) : true) &&
      count(t -> isinf(t), gzero(gminreal(syso))) == 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 sysi, syso, info = giofac(sys, minphase = false, atol1 = 1.e-7, atol2 = 1.e-7)
@test gnrank(sys-sysi*syso,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysi'*sysi-I) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysi) & isproper(syso) && # checking properness of factors
      isstable(sysi) && (isstable(sys) ? isstable(syso) : true) &&
      count(t -> isinf(t), gzero(gminreal(syso))) == 1 &&
      info.nrank == 1 && ismissing(info.nfuz) && info.niuz == 0

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

sys = rss(0,1,0);
@time sysi, syso, info = giofac(sys, minphase = false)
r = size(syso,1);
@test gnrank(sys-sysi[:,1:r]*syso,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysi'*sysi-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysi) & (isproper(sys) ? isproper(syso) : true) &&
      isstable(sysi) && (isstable(sys) ? isstable(syso) : true) &&
      isempty(gzero(gminreal(syso))) &&
      info.nrank == r && ismissing(info.nfuz) && info.niuz == 0

sys = rss(0,0,2);
@time sysi, syso, info = giofac(sys, minphase = false)
r = size(syso,1);
@test gnrank(sys-sysi[:,1:r]*syso,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysi'*sysi-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysi) & (isproper(sys) ? isproper(syso) : true) &&
      isstable(sysi) && (isstable(sys) ? isstable(syso) : true) &&
      isempty(gzero(gminreal(syso))) &&
      info.nrank == r && ismissing(info.nfuz) && info.niuz == 0

sys = dss([1 0]);
@time sysi, syso, info = giofac(sys, minphase = false)
r = size(syso,1);
@test gnrank(sys-sysi[:,1:r]*syso,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysi'*sysi-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysi) & (isproper(sys) ? isproper(syso) : true) &&
      isstable(sysi) && (isstable(sys) ? isstable(syso) : true) &&
      isempty(gzero(gminreal(syso))) &&
      info.nrank == r && ismissing(info.nfuz) && info.niuz == 0

sys = dss([0; 1]);
@time sysi, syso, info = giofac(sys, minphase = false)
r = size(syso,1);
@test gnrank(sys-sysi[:,1:r]*syso,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysi'*sysi-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysi) & (isproper(sys) ? isproper(syso) : true) &&
      isstable(sysi) && (isstable(sys) ? isstable(syso) : true) &&
      isempty(gzero(gminreal(syso))) &&
      info.nrank == r && ismissing(info.nfuz) && info.niuz == 0

sys = rss(0,3,2);
@time sysi, syso, info = giofac(sys, minphase = false)
r = size(syso,1);
@test gnrank(sys-sysi[:,1:r]*syso,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysi'*sysi-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysi) & (isproper(sys) ? isproper(syso) : true) &&
      isstable(sysi) && (isstable(sys) ? isstable(syso) : true) &&
      isempty(gzero(gminreal(syso))) &&
      info.nrank == r && ismissing(info.nfuz) && info.niuz == 0

z = Polynomial([0, 1],'z');     # define the complex variable z

# Gd = [z^2 z/(z-2); 0 1/z]     # define the 2-by-2 improper Gd(z)
Nd = [z^2 z; 0 1]; Dd = [1 z-2; 1 z];

# build LTI minimal descriptor realizations of Gd(z)
sysd = dss(Nd,Dd,minimal = true,Ts = 1);
zeref = gzero(sysd,atol1=1.e-7)

sysi, syso, info = giofac(sysd, atol = 1.e-7)
r = size(syso,1);
zer = gzero(syso,atol1=1.e-7)
@test gnrank(sysd-sysi[:,1:r]*syso,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysi'*sysi-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysi) && (isproper(sysd) ? isproper(syso) : true) && # checking properness of factors
      isstable(sysi) && (isstable(sysd) ? isstable(syso) : true) &&
      count(t -> isinf(t), zer) == info.niuz &&
      count(t -> (abs.(t) .>= 1 && isfinite(t)), zer) == 0 &&
      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)];

sys = dss(gn,gd,minimal = true, atol = 1.e-7);
sysi, syso, info = giofac(sys, atol = 1.e-7)
r = size(syso,1);
zer = gzero(gminreal(syso),atol1=1.e-7);
@test gnrank(sys-sysi[:,1:r]*syso,atol1=1.e-7) == 0   &&   #  G(s) = Gi(s)*Go(s) = 0
      gnrank(sysi'*sysi-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysi) & isproper(syso) && # checking properness of factors
      isstable(sysi) && (isstable(sys) ? isstable(syso) : 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

sysi, syso, info = giofac(sys, minphase = false, atol = 1.e-7)
zeref = gzero(gminreal(sys),atol1=1.e-7)
r = size(syso,1);
zer = gzero(gminreal(syso),atol1=1.e-7);
@test gnrank(sys-sysi[:,1:r]*syso,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysi'*sysi-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysi) && isproper(syso) && # checking properness of factors
      isstable(sysi) && (isstable(sys) ? isstable(syso) : 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);
sysi, syso, info = giofac(sys, atol = 1.e-7)
zeref = gzero(sys,atol1=1.e-7)
r = size(syso,1);
zer = gzero(gminreal(syso),atol1=1.e-7)
@test gnrank(sys-sysi[:,1:r]*syso,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysi'*sysi-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysi) && (isproper(sys) ? isproper(syso) : true) && # checking properness of factors
      isstable(sysi) && (isstable(sys) ? isstable(syso) : 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)
sysi, syso, info = giofac(sys, atol = 1.e-7)
zeref = gzero(sys,atol1=1.e-7)
r = size(syso,1)
zer = gzero(gminreal(syso),atol1=1.e-7)
@test gnrank(sys-sysi[:,1:r]*syso,atol1=1.e-7) == 0   &&   #  G(z) - Gi(z)*Go(z) = 0
      gnrank(sysi'*sysi-I,atol1=1.e-7) == 0  && # conj(Gi(z))*Gi(z)-I = 0
      isproper(sysi) && (isproper(sys) ? isproper(syso) : true) && # checking properness of factors
      isstable(sysi) && (isstable(sys) ? isstable(syso) : 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)
sysi, syso, info = giofac(sys, atol = 1.e-7)
@test gnrank(sys-sysi*syso,atol1=1.e-7) == 0   &&
      gnrank(sysi'*sysi-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];
# error
sys = dss(gn,gd,minimal = true, atol = 1.e-7,Ts = 1);
sysi, syso, info = giofac(sys, atol = 1.e-7)
zeref = gzero(sys,atol1=1.e-7)
r = size(syso,1)
zer = gzero(gminreal(syso),atol1=1.e-7)
@test gnrank(sys-sysi[:,1:r]*syso,atol1=1.e-7) == 0   &&   #  G(z) - Gi(z)*Go(z) = 0
      gnrank(sysi'*sysi-I,atol1=1.e-7) == 0  && # conj(Gi(z))*Gi(z)-I = 0
      isproper(sysi) && (isproper(sys) ? isproper(syso) : true) && # checking properness of factors
      isstable(sysi) && (isstable(sys) ? isstable(syso) : 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);
@time sysi, syso, info = giofac(sys, fast = fast, atol = 1.e-7) ;
zeref = gzero(sys,atol1=1.e-7)
r = size(syso,1);
zer = gzero(gminreal(syso),atol1=1.e-7)
@test gnrank(sys-sysi[:,1:r]*syso,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysi'*sysi-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysi) && (isproper(sys) ? isproper(syso) : true) && # checking properness of factors
      isstable(sysi) && (isstable(sys) ? isstable(syso) : 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, standard
sys = rss(n,p,m,T = Ty,disc=true,stable=true);
@time sysi, syso, info = giofac(sys, fast = fast, atol = 1.e-7) ;
zeref = gzero(sys,atol1=1.e-7)
r = size(syso,1);
zer = gzero(gminreal(syso),atol1=1.e-7)
@test gnrank(sys-sysi[:,1:r]*syso,atol1=1.e-7) == 0   &&   #  G(z) - Gi(z)*Go(z) = 0
      gnrank(sysi'*sysi-I,atol1=1.e-7) == 0  && # conj(Gi(u))*Gi(z)-I = 0
      isproper(sysi) && (isproper(sys) ? isproper(syso) : true) && # checking properness of factors
      isstable(sysi) && (isstable(sys) ? isstable(syso) : 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 sysi, syso, info = giofac(sys, fast = fast, atol = 1.e-7) ;
zeref = gzero(sys,atol1=1.e-7)
r = size(syso,1);
zer = gzero(gminreal(syso),atol1=1.e-7)
@test gnrank(sys-sysi[:,1:r]*syso,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysi'*sysi-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysi) && (isproper(sys) ? isproper(syso) : true) && # checking properness of factors
      isstable(sysi) && (isstable(sys) ? isstable(syso) : 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 sysi, syso, info = giofac(sys, fast = fast, atol = 1.e-7) ;
zeref = gzero(sys,atol1=1.e-7)
r = size(syso,1);
zer = gzero(gminreal(syso),atol1=1.e-7)
@test gnrank(sys-sysi[:,1:r]*syso,atol1=1.e-7) == 0   &&   #  G(z) - Gi(z)*Go(z) = 0
      gnrank(sysi'*sysi-I,atol1=1.e-7) == 0  && # conj(Gi(u))*Gi(z)-I = 0
      isproper(sysi) && (isproper(sys) ? isproper(syso) : true) && # checking properness of factors
      isstable(sysi) && (isstable(sys) ? isstable(syso) : 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 sysi, syso, info = giofac(sys, fast = fast, atol = 1.e-7) ;
zeref = gzero(sys,atol1=1.e-7)
r = size(syso,1);
zer = gzero(gminreal(syso),atol1=1.e-7)
@test gnrank(sys-sysi[:,1:r]*syso,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysi'*sysi-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysi) && (isproper(sys) ? isproper(syso) : true) && # checking properness of factors
      isstable(sysi) && (isstable(sys) ? isstable(syso) : 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 sysi, syso, info = giofac(sys, fast = fast, atol = 1.e-7) ;
zeref = gzero(sys,atol1=1.e-7)
r = size(syso,1);
zer = gzero(gminreal(syso),atol1=1.e-7)
@test gnrank(sys-sysi[:,1:r]*syso,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysi'*sysi-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysi) && (isproper(sys) ? isproper(syso) : true) && # checking properness of factors
      isstable(sysi) && (isstable(sys) ? isstable(syso) : 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 sysi, syso, info = giofac(sys, fast = fast, atol = 1.e-7) ;
zeref = gzero(sys,atol1=1.e-7)
r = size(syso,1);
zer = gzero(gminreal(syso),atol1=1.e-7)
@test gnrank(sys-sysi[:,1:r]*syso,atol1=1.e-7) == 0   &&   #  G(z) - Gi(z)*Go(z) = 0
      gnrank(sysi'*sysi-I,atol1=1.e-7) == 0  && # conj(Gi(u))*Gi(z)-I = 0
      isproper(sysi) && (isproper(sys) ? isproper(syso) : true) && # checking properness of factors
      isstable(sysi) && (isstable(sys) ? isstable(syso) : 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 sysi, syso, info = giofac(sys, fast = fast, atol = 1.e-7) ;
zeref = gzero(sys,atol1=1.e-7)
r = size(syso,1);
zer = gzero(gminreal(syso),atol1=1.e-7)
@test gnrank(sys-sysi[:,1:r]*syso,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysi'*sysi-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysi) && (isproper(sys) ? isproper(syso) : true) && # checking properness of factors
      isstable(sysi) && (isstable(sys) ? isstable(syso) : 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 sysi, syso, info = giofac(sys, fast = fast, atol = 1.e-7, minphase=true) ;
zeref = gzero(sys,atol1=1.e-7)
r = size(syso,1);
zer = gzero(gminreal(syso),atol1=1.e-7)
@test gnrank(sys-sysi[:,1:r]*syso,atol1=1.e-7) == 0   #&&   #  G(z) - Gi(z)*Go(z) = 0
@test       iszero(sysi'*sysi-I,atol1=1.e-7)  #&& # conj(Gi(u))*Gi(z)-I = 0
@test       isproper(sysi) && (isproper(sys) ? isproper(syso) : true) #&& # checking properness of factors
@test       isstable(sysi) && (isstable(sys) ? isstable(syso) : 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 sysi, syso, info = giofac(sys, fast = fast, atol = 1.e-7) ;
zeref = gzero(sys,atol1=1.e-7)
r = size(syso,1);
zer = gzero(gminreal(syso),atol1=1.e-7)
@test gnrank(sys-sysi[:,1:r]*syso,atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysi'*sysi-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysi) && (isproper(sys) ? isproper(syso) : true) && # checking properness of factors
      isstable(sysi) && (isstable(sys) ? isstable(syso) : 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 sysi, syso, info = giofac(sys, fast = fast, atol = 1.e-7) ;
zeref = gzero(gminreal(sys),atol1=1.e-7,atol2=1.e-7)
r = size(syso,1);
zer = gzero(gminreal(syso),atol1=1.e-7,atol2=1.e-7)
@test gnrank(sys-sysi[:,1:r]*syso,atol1=1.e-7) == 0   &&   #  G(z) - Gi(z)*Go(z) = 0
      iszero(sysi'*sysi-I,atol1=1.e-5)  && # conj(Gi(u))*Gi(z)-I = 0
      isproper(sysi) && (isproper(sys) ? isproper(syso) : true) && # checking properness of factors
      isstable(sysi) && (isstable(sys) ? isstable(syso) : 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 # giofac

@testset "goifac" begin

sys = rdss(0,0,0);
@time sysi, syso, info = goifac(sys, minphase = false)
@test gnrank(sys-syso*sysi) == 0   &&   #  G(s) - Go(s)*Gi(s) = 0
       gnrank(sysi'*sysi-I,atol=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
       isproper(sysi) & (isproper(sys) ? isproper(syso) : true) &&
       isstable(sysi) && (isstable(sys) ? isstable(syso) : true) &&
       isempty(gzero(gminreal(syso))) &&
       info.nrank == 0 && ismissing(info.nfuz) && info.niuz == 0

fast = true; Ty = Complex{Float64}; Ty = Float64
n = 5; p = 2; m = 4;
for fast in (true, false)
# random examples
for Ty in (Float64, Complex{Float64})

# continuous, standard
sys = rss(n,p,m,T = Ty,disc=false);
@time sysi, syso, info = goifac(sys, fast = fast, atol = 1.e-7) ;
zeref = gzero(sys,atol1=1.e-7)
r = size(syso,1);
zer = gzero(gminreal(syso),atol1=1.e-7)
@test gnrank(sys-syso*sysi[1:r,:],atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysi'*sysi-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysi) && (isproper(sys) ? isproper(syso) : true) && # checking properness of factors
      isstable(sysi) && (isstable(sys) ? isstable(syso) : 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, standard
sys = rss(n,p,m,T = Ty,disc=true);
@time sysi, syso, info = goifac(sys, fast = fast, atol = 1.e-7) ;
zeref = gzero(sys,atol1=1.e-7)
r = size(syso,1);
zer = gzero(gminreal(syso),atol1=1.e-7)
@test gnrank(sys-syso*sysi[1:r,:],atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysi'*sysi-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysi) && (isproper(sys) ? isproper(syso) : true) && # checking properness of factors
      isstable(sysi) && (isstable(sys) ? isstable(syso) : 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, infinite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=false,id=[3*ones(Int,1);2*ones(Int,1)],randlt=true,randrt=false);
@time sysi, syso, info = goifac(sys, fast = fast, atol = 1.e-7) ; info
zeref = gzero(sys,atol=1.e-7)
r = size(syso,1);
zer = gzero(gminreal(syso),atol=1.e-7)
@test gnrank(sys-syso*sysi[1:r,:],atol1=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
      gnrank(sysi'*sysi-I,atol1=1.e-7) == 0  && # conj(Gi(s))*Gi(s)-I = 0
      isproper(sysi) && (isproper(sys) ? isproper(syso) : true) && # checking properness of factors
      isstable(sysi) && (isstable(sys) ? isstable(syso) : 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 eigenvalues
# sys = rdss(n,p,m,T = Ty,disc=true,id=[ones(Int,3);2*ones(Int,2)],randlt=true,randrt=true);
# @time sysi, syso, info = goifac(sys, fast = fast, atol = 1.e-7) ; info
# zeref = gzero(sys,atol=1.e-7)
# r = size(syso,1);
# zer = gzero(gminreal(syso,atol=1.e-7),atol=1.e-7)
# @test gnrank(sys-syso*sysi[1:r,:],atol=1.e-7) == 0   &&   #  G(s) - Gi(s)*Go(s) = 0
#       iszero(sysi'*sysi-I,atol=1.e-5)  && # conj(Gi(s))*Gi(s)-I = 0
#       isproper(sysi) && (isproper(sys) ? isproper(syso) : true) && # checking properness of factors
#       isstable(sysi) && (isstable(sys) ? isstable(syso) : true) &&
#       count(isinf.(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

end
end

end # goifac

end #test

end #module


text:
using CSV, DataFrames, FeatherFiles

# Data from Gouda, The Netherlands

FILE_NAME = "./data/1_2008Wind_Final_hub(m)_80_lat_long_52.0 4.7.csv"
df1 = CSV.File(FILE_NAME, normalizenames=true) |> DataFrame!

FILE_NAME = "./data/2_2008Wind_Final_hub(m)_80_lat_long_52.0 4.7.csv"
df2 = CSV.File(FILE_NAME, normalizenames=true) |> DataFrame!

FILE_NAME = "./data/3_2008Wind_Final_hub(m)_80_lat_long_52.0 4.7.csv"
df3 = CSV.File(FILE_NAME, normalizenames=true) |> DataFrame!

FILE_NAME = "./data/4_2008Wind_Final_hub(m)_80_lat_long_52.0 4.7.csv"
df4 = CSV.File(FILE_NAME, normalizenames=true) |> DataFrame!

FILE_NAME = "./data/5_2008Wind_Final_hub(m)_80_lat_long_52.0 4.7.csv"
df5 = CSV.File(FILE_NAME, normalizenames=true) |> DataFrame!

FILE_NAME = "./data/6_2008Wind_Final_hub(m)_80_lat_long_52.0 4.7.csv"
df6 = CSV.File(FILE_NAME, normalizenames=true) |> DataFrame!

FILE_NAME = "./data/7_2008Wind_Final_hub(m)_80_lat_long_52.0 4.7.csv"
df7 = CSV.File(FILE_NAME, normalizenames=true) |> DataFrame!

FILE_NAME = "./data/8_2008Wind_Final_hub(m)_80_lat_long_52.0 4.7.csv"
df8 = CSV.File(FILE_NAME, normalizenames=true) |> DataFrame!

FILE_NAME = "./data/9_2008Wind_Final_hub(m)_80_lat_long_52.0 4.7.csv"
df9 = CSV.File(FILE_NAME, normalizenames=true) |> DataFrame!

FILE_NAME = "./data/10_2008Wind_Final_hub(m)_80_lat_long_52.0 4.7.csv"
df10 = CSV.File(FILE_NAME, normalizenames=true) |> DataFrame!

FILE_NAME = "./data/11_2008Wind_Final_hub(m)_80_lat_long_52.0 4.7.csv"
df11 = CSV.File(FILE_NAME, normalizenames=true) |> DataFrame!

FILE_NAME = "./data/12_2008Wind_Final_hub(m)_80_lat_long_52.0 4.7.csv"
df12 = CSV.File(FILE_NAME, normalizenames=true) |> DataFrame!

df = [df1;df2;df3;df4;df5;df6;df7;df8;df9;df10;df11;df12]
save("./data/windyear.feather", df)

ext: jl


text:
__precompile__(false)

module ReplMaker

import REPL
import REPL.LineEdit

export initrepl

"""
```
         initrepl(parser::Function;
                  prompt_text = "myrepl> ",
                  prompt_color = :blue,
                  start_key = ')',
                  repl = Base.active_repl,
                  mode_name = :mylang,
                  valid_input_checker::Function = (s -> true),
                  keymap::Dict = REPL.LineEdit.default_keymap_dict,
                  completion_provider = REPL.REPLCompletionProvider(),
                  sticky = true,
                  startup_text=true
                  )
```
creates a custom repl mode which takes in code and parses it according to whatever parsing function you
provide in the argument `parser`. Choose which key initializes the repl mode with `start_key`, the name of
your repl mode with `mode_name` and optionally provide a function which checks if a given repl input is valid
before parsing with `valid_input_checker`. Autocompletion options are supplied through the argument
`completion_provider` which defaults to the standard julia REPL TAB completions.
"""
function initrepl(parser::Function;
                  prompt_text = "myrepl> ",
                  prompt_color = :blue,
                  start_key = ')',
                  repl = Base.active_repl,
                  mode_name = :mylang,
                  valid_input_checker::Function = (s -> true),
                  keymap::Dict = REPL.LineEdit.default_keymap_dict,
                  completion_provider = REPL.REPLCompletionProvider(),
                  sticky_mode = true,
                  startup_text=true
                  )

    color = Base.text_colors[prompt_color]

    julia_mode = repl.interface.modes[1]
    prefix = repl.hascolor ? color : ""
    suffix = repl.hascolor ? (repl.envcolors ? Base.input_color : repl.input_color()) : ""

    lang_mode = LineEdit.Prompt(prompt_text;
                                prompt_prefix    = prefix,
                                prompt_suffix    = suffix,
                                keymap_dict      = keymap,
                                on_enter         = valid_input_checker,
                                complete         = completion_provider,
                                sticky           = sticky_mode
                                )

    lang_mode.on_done = REPL.respond(parser, repl, lang_mode)

    push!(repl.interface.modes, lang_mode)

    hp                     = julia_mode.hist
    hp.mode_mapping[mode_name |> Symbol] = lang_mode
    lang_mode.hist         = hp

    search_prompt, skeymap = LineEdit.setup_search_keymap(hp)

    mk = REPL.mode_keymap(julia_mode)

    lang_keymap = Dict(
    start_key => (s, args...) ->
      if isempty(s) || position(LineEdit.buffer(s)) == 0
        buf = copy(LineEdit.buffer(s))
        LineEdit.transition(s, lang_mode) do
          LineEdit.state(s, lang_mode).input_buffer = buf
        end
      else
        LineEdit.edit_insert(s, start_key)
      end
    )

    lang_mode.keymap_dict = LineEdit.keymap(Dict[
        skeymap,
        mk,
        LineEdit.history_keymap,
        LineEdit.default_keymap,
        LineEdit.escape_defaults,
    ])
    if start_key in keys(julia_mode.keymap_dict)
        @warn "REPL key '$start_key' overwritten."
    end

    julia_mode.keymap_dict = LineEdit.keymap_merge(julia_mode.keymap_dict, lang_keymap)

    startup_text && println("REPL mode $mode_name initialized. Press $start_key to enter and backspace to exit.")

  nothing
end

end


text:
global const slope_soft  = 1.e-6
global const c_light     = 2.9979e10
global const qe          = 4.8032e-10
global const m_e         = 9.10953e-28
global const nu_c_prefac = 3.0 * qe / (4π * m_e^3 * c_light^5)
global const C_crit      = 3.0 * qe / (4π * m_e^3 * c_light^3)


text:
if VERSION < v"1.2"
    Base.getproperty(x::Tuple, f::Int) = getfield(x, f)
end

if VERSION < v"1.1"
    # Note: these are actually *better* than the ones in julia 1.1, 1.2, 1.3,and 1.4
    # See: https://github.com/JuliaLang/julia/issues/34292
    function fieldtypes(::Type{T}) where {T}
        if @generated
            ntuple(i -> fieldtype(T, i), fieldcount(T))
        else
            ntuple(i -> fieldtype(T, i), fieldcount(T))
        end
    end

    function fieldnames(::Type{T}) where {T}
        if @generated
            ntuple(i -> fieldname(T, i), fieldcount(T))
        else
            ntuple(i -> fieldname(T, i), fieldcount(T))
        end
    end
end


text:
#
# Broadcasting for IndexSets
#

# TODO: delete
## # We're using a specialized type `IndexVector` since `IndexSet` has
## # a compile time check that all indices are unique, but `similar`
## # in general won't make a unique IndexSet. We need a specialized type
## # to dispatch on `copyto!`. This is like
## struct IndexVector{T} <: AbstractVector{T}
##   data::Vector{T}
## end
## size(v::IndexVector) = size(v.data)
##
## struct IndexSetStyle <: Broadcast.BroadcastStyle end
##
## BroadcastStyle(::Type{<: IndexSet}) = IndexSetStyle()
##
## BroadcastStyle(::IndexSetStyle, ::BroadcastStyle) = IndexSetStyle()
##
## broadcastable(is::IndexSet) = is
##
## function _similar(bc::Broadcasted{IndexSetStyle}, ::Type{ElT}) where {ElT}
##   return similar(Array{ElT}, axes(bc))
## end
##
## # In the case when the output type is inferred to be `<: Index`,
## # then output an IndexSet (like `prime.(is)`)
## function similar(bc::Broadcasted{IndexSetStyle}, ::Type{ElT}) where {ElT <: Index}
##   is = _similar(bc, ElT)
##   # We're using a specialized type `IndexVector` since `IndexSet` has
##   # a compile time check that all indices are unique, but `similar`
##   # in general won't make a unique IndexSet. We need a specialized type
##   # to dispatch on `copyto!`
##   return IndexVector(_similar(bc, ElT))
## end
##
## # In general, the output type will be a Vector{ElT} (like `dim.(is)`)
## function similar(bc::Broadcast.Broadcasted{IndexSetStyle}, ::Type{ElT}) where {ElT}
##   return _similar(bc, ElT)
## end
##
## function copyto!(dest::IndexVector, bc::Broadcast.Broadcasted{IndexSetStyle})
##   copyto!(dest.data, bc)
##   return IndexSet(dest.data)
## end

#
# Broadcasting for ITensors
#

#
# ITensorStyle
#

struct ITensorStyle <: BroadcastStyle end

BroadcastStyle(::Type{<:ITensor}) = ITensorStyle()

broadcastable(T::ITensor) = T

function Base.similar(bc::Broadcasted{ITensorStyle}, ::Type{ElT}) where {ElT<:Number}
  A = find_type(ITensor, bc.args)
  return similar(A, ElT)
end

#
# ITensorOpScalarStyle
# Operating with a scalar
#

struct ITensorOpScalarStyle <: BroadcastStyle end

function Base.BroadcastStyle(::ITensorStyle, ::DefaultArrayStyle{0})
  return ITensorOpScalarStyle()
end

Base.BroadcastStyle(::ITensorStyle, ::ITensorOpScalarStyle) = ITensorOpScalarStyle()

instantiate(bc::Broadcasted{ITensorOpScalarStyle}) = bc

function broadcasted(::typeof(Base.literal_pow), ::typeof(^), T::ITensor, x::Val)
  return broadcasted(ITensorOpScalarStyle(), Base.literal_pow, Ref(^), T, Ref(x))
end

function Base.similar(
  bc::Broadcasted{ITensorOpScalarStyle}, ::Type{ElT}
) where {ElT<:Number}
  A = find_type(ITensor, bc.args)
  return similar(A, ElT)
end

#
# For arbitrary function chaining f.(g.(h.(x)))
#

function instantiate(bc::Broadcasted{ITensorStyle,<:Any,<:Function,<:Tuple{Broadcasted}})
  return instantiate(broadcasted(bc.f ∘ bc.args[1].f, bc.args[1].args...))
end

function instantiate(
  bc::Broadcasted{
    ITensorStyle,
    <:Any,
    <:Function,
    <:Tuple{Broadcasted{ITensorStyle,<:Any,<:Function,<:Tuple{<:ITensor}}},
  },
)
  return broadcasted(bc.f ∘ bc.args[1].f, bc.args[1].args...)
end

instantiate(bc::Broadcasted{ITensorStyle}) = bc

#
# Some helper functionality to find certain
# inputs in the argument list
#

"`A = find_type(::Type,As)` returns the first of type Type among the arguments."
function find_type(::Type{T}, args::Tuple) where {T}
  return find_type(T, find_type(T, args[1]), Base.tail(args))
end
find_type(::Type{T}, x) where {T} = x
find_type(::Type{T}, a::T, rest) where {T} = a
find_type(::Type{T}, ::Any, rest) where {T} = find_type(T, rest)
# If not found, return nothing
find_type(::Type{T}, ::Tuple{}) where {T} = nothing

#
# Generic fallback
#

function Base.copyto!(T::ITensor, bc::Broadcasted)
  @show typeof(bc)
  return error(
    "The broadcasting operation you are attempting is not yet implemented for ITensors, please raise an issue if you would like it to be supported.",
  )
end

#
# B .= α .* A
# A .*= α
#

function Base.copyto!(
  T::ITensor,
  bc::Broadcasted{
    ITensorOpScalarStyle,
    <:Any,
    typeof(*),
    <:Tuple{<:Union{<:Number,<:ITensor},<:Union{<:Number,<:ITensor}},
  },
)
  α = find_type(Number, bc.args)
  A = find_type(ITensor, bc.args)
  map!((t, a) -> α * a, T, T, A)
  return T
end

#
# B .= A ./ α
# A ./= α
#

function Base.copyto!(
  T::ITensor,
  bc::Broadcasted{ITensorOpScalarStyle,<:Any,typeof(/),<:Tuple{<:ITensor,<:Number}},
)
  α = find_type(Number, bc.args)
  A = find_type(ITensor, bc.args)
  map!((t, a) -> bc.f(a, α), T, T, A)
  return T
end

#
# C .= A ./ B
#

function Base.copyto!(
  R::ITensor, bc::Broadcasted{ITensorStyle,<:Any,typeof(/),<:Tuple{<:ITensor,<:ITensor}}
)
  T1, T2 = bc.args
  if R === T1
    map!((t1, t2) -> bc.f(t1, t2), R, T1, T2)
  elseif R === T2
    map!((t1, t2) -> bc.f(t2, t1), R, T2, T1)
  else
    error("When dividing two ITensors in-place, one must be the same as the output ITensor")
  end
  return R
end

#
# C .= A .⊙ B
#

function Base.copyto!(
  R::ITensor, bc::Broadcasted{ITensorStyle,<:Any,typeof(⊙),<:Tuple{<:ITensor,<:ITensor}}
)
  T1, T2 = bc.args
  if R === T1
    map!((t1, t2) -> *(t1, t2), R, T1, T2)
  elseif R === T2
    map!((t1, t2) -> *(t2, t1), R, T2, T1)
  else
    error(
      "When Hadamard producting two ITensors in-place, one must be the same as the output ITensor",
    )
  end
  return R
end

#
# B .= α ./ A
#

function Base.copyto!(
  T::ITensor,
  bc::Broadcasted{ITensorOpScalarStyle,<:Any,typeof(/),<:Tuple{<:Number,<:ITensor}},
)
  α = find_type(Number, bc.args)
  A = find_type(ITensor, bc.args)
  map!((t, a) -> bc.f(α, a), T, T, A)
  return T
end

#
# For B .= A .^ 2.5
#

function Base.copyto!(R::ITensor, bc::Broadcasted{ITensorOpScalarStyle,<:Any,typeof(^)})
  α = find_type(Number, bc.args)
  T = find_type(ITensor, bc.args)
  map!((r, t) -> t^α, R, R, T)
  return R
end

#
# For B .= A .^ 2
#

function Base.copyto!(
  R::ITensor, bc::Broadcasted{ITensorOpScalarStyle,<:Any,typeof(Base.literal_pow)}
)
  α = find_type(Base.RefValue{<:Val}, bc.args).x
  powf = find_type(Base.RefValue{<:Function}, bc.args).x
  @assert !isnothing(powf)
  T = find_type(ITensor, bc.args)
  map!((r, t) -> bc.f(^, t, α), R, R, T)
  return R
end

#
# For A .= α
#

function Base.copyto!(
  T::ITensor, bc::Broadcasted{DefaultArrayStyle{0},<:Any,typeof(identity),<:Tuple{<:Number}}
)
  fill!(T, bc.args[1])
  return T
end

#
# For B .= A
#

function Base.copyto!(
  T::ITensor, bc::Broadcasted{ITensorStyle,<:Any,typeof(identity),<:Tuple{<:ITensor}}
)
  A = bc.args[1]
  map!((r, t) -> t, T, T, A)
  return T
end

#
# B .+= A
#

function fmap(bc::Broadcasted{ITensorStyle,<:Any,typeof(+),<:Tuple{Vararg{<:ITensor}}})
  return (r, t) -> bc.f(r, t)
end

function fmap(bc::Broadcasted{ITensorStyle,<:Any,typeof(-),<:Tuple{Vararg{<:ITensor}}})
  return (r, t) -> bc.f(r, t)
end

function Base.copyto!(
  T::ITensor, bc::Broadcasted{ITensorStyle,<:Any,typeof(+),<:Tuple{Vararg{<:ITensor}}}
)
  if T === bc.args[1]
    A = bc.args[2]
  elseif T === bc.args[2]
    A = bc.args[1]
  else
    error("When adding two ITensors in-place, one must be the same as the output ITensor")
  end
  map!(fmap(bc), T, T, A)
  return T
end

#
# B .-= A
#

function Base.copyto!(
  T::ITensor, bc::Broadcasted{ITensorStyle,<:Any,typeof(-),<:Tuple{Vararg{<:ITensor}}}
)
  if T === bc.args[1]
    A = bc.args[2]
  elseif T === bc.args[2]
    A = bc.args[1]
  else
    error(
      "When subtracting two ITensors in-place, one must be the same as the output ITensor"
    )
  end
  map!(fmap(bc), T, T, A)
  return T
end

#
# C .+= α .* A
# C .-= α .* A
#
# C .+= α .* A .* B
# C .-= α .* A .* B
#
# B .+= A .^ 2.5
# B .-= A .^ 2.5
#
# B .+= A .^ 2
# B .-= A .^ 2
#

#

function Base.copyto!(
  T::ITensor, bc::Broadcasted{ITensorOpScalarStyle,<:Any,<:Union{typeof(+),typeof(-)}}
)
  C = find_type(ITensor, bc.args)
  bc_bc = find_type(Broadcasted, bc.args)
  if T === C
    A = find_type(ITensor, bc_bc.args)
    α = find_type(Number, bc_bc.args)

    # Check if it is the case .^(::Int)
    γ = find_type(Base.RefValue{<:Val}, bc_bc.args)
    powf = find_type(Base.RefValue{<:Function}, bc_bc.args)

    if !isnothing(α) && !isnothing(A)
      map!((r, t) -> bc.f(r, bc_bc.f(t, α)), T, T, A)
    elseif !isnothing(γ) && !isnothing(A) && !isnothing(powf)
      map!((r, t) -> bc.f(r, bc_bc.f(powf[], t, γ[])), T, T, A)
    else
      # In-place contraction:
      # C .+= α .* A .* B
      bc_bc_bc = find_type(Broadcasted, bc_bc.args)
      if isnothing(α)
        β = find_type(Number, bc_bc_bc.args)
        B = find_type(ITensor, bc_bc_bc.args)
      else
        A, B = bc_bc_bc.args
      end
      mul!(T, A, B, β, bc.f(1))
    end
  else
    error("When adding two ITensors in-place, one must be the same as the output ITensor")
  end
  return T
end

#
# C .= β .* C .+ α .* A
# C .= β .* C .+ α .* A .* B
#

function Base.copyto!(
  T::ITensor,
  bc::Broadcasted{ITensorOpScalarStyle,<:Any,typeof(+),<:Tuple{Vararg{<:Broadcasted}}},
)
  bc_α = bc.args[1]
  bc_β = bc.args[2]
  α = find_type(Number, bc_α.args)
  A = find_type(ITensor, bc_α.args)
  β = find_type(Number, bc_β.args)
  C = find_type(ITensor, bc_β.args)
  (T !== A && T !== C) &&
    error("When adding two ITensors in-place, one must be the same as the output ITensor")
  if T === A
    bc_α, bc_β = bc_β, bc_α
    α, β = β, α
    A, C = C, A
  end
  if !isnothing(A) && !isnothing(C) && !isnothing(α) && !isnothing(β)
    map!((r, t) -> β * r + α * t, T, T, A)
  else
    bc_bc_α = find_type(Broadcasted, bc_α.args)
    if isnothing(α)
      α = find_type(Number, bc_bc_α.args)
      B = find_type(ITensor, bc_bc_α.args)
    else
      A, B = bc_bc_α.args
    end
    mul!(T, A, B, α, β)
  end
  return T
end

#
# For A .+= α
#

function Base.copyto!(
  T::ITensor,
  bc::Broadcasted{
    ITensorOpScalarStyle,<:Any,typeof(+),<:Tuple{Vararg{<:Union{<:ITensor,<:Number}}}
  },
)
  α = find_type(Number, bc.args)
  A = find_type(ITensor, bc.args)
  if A === T
    tensor(T) .= tensor(A) .+ α
  else
    error(
      "Currently, we don't support `B .= A .+ α` if `B !== A` (i.e. only `A .+= α` is supported",
    )
  end
  return T
end

#
# For C .= A .* B
#

function Base.copyto!(
  T::ITensor, bc::Broadcasted{ITensorStyle,<:Any,typeof(*),<:Tuple{<:ITensor,<:ITensor}}
)
  mul!(T, bc.args[1], bc.args[2])
  return T
end

#
# For C .= α .* A .* B
#

function Base.copyto!(T::ITensor, bc::Broadcasted{ITensorOpScalarStyle,<:Any,typeof(*)})
  A = find_type(Union{<:Number,<:ITensor}, bc.args)
  bc_bc = find_type(Broadcasted, bc.args)
  if A isa Number
    mul!(T, bc_bc.args[1], bc_bc.args[2], A)
  else
    mul!(T, A, find_type(ITensor, bc_bc.args), find_type(Number, bc_bc.args))
  end
  return T
end

#
# For B .= f.(A)
#

function Base.copyto!(
  R::ITensor, bc::Broadcasted{ITensorStyle,<:Any,<:Function,<:Tuple{<:ITensor}}
)
  f = bc.f
  T = bc.args[1]
  map!((r, t) -> f(t), R, R, T)
  return R
end

#
# For B .+= f.(A)
#

function Base.copyto!(
  R::ITensor,
  bc::Broadcasted{
    ITensorStyle,<:Any,typeof(+),<:Tuple{Vararg{Union{<:ITensor,<:Broadcasted}}}
  },
)
  R̃ = find_type(ITensor, bc.args)
  bc2 = find_type(Broadcasted, bc.args)
  if R === R̃
    map!((r, t) -> r + bc2.f(t), R, R, bc2.args[1])
  else
    error("In C .= B .+ f.(A), C and B must be the same ITensor")
  end
  return R
end

#
# For B .= f.(B) + g.(A)
#

function Base.copyto!(
  R::ITensor, bc::Broadcasted{ITensorStyle,<:Any,typeof(+),<:Tuple{Vararg{<:Broadcasted}}}
)
  bc1 = bc.args[1]
  bc2 = bc.args[2]
  T1 = bc1.args[1]
  f1 = bc1.f
  T2 = bc2.args[1]
  f2 = bc2.f
  if R === T1
    map!((r, t) -> f1(r) + f2(t), R, R, T2)
  elseif R === T2
    map!((r, t) -> f2(r) + f1(t), R, R, T1)
  else
    error("In C .= f.(B) .+ g.(A), C and B or A must be the same ITensor")
  end
  return R
end


text:
mutable struct GFProjectState
    constraints::Any
    trace::GFTrace
    score::Float64
    visitor::AddressVisitor
    params::Dict{Symbol,Any}
    discard::DynamicAssignment
end

function GFProjectState(constraints, params)
    visitor = AddressVisitor()
    discard = DynamicAssignment()
    GFProjectState(constraints, GFTrace(), 0., visitor, params, discard)
end

function addr(state::GFProjectState, dist::Distribution{T}, args, addr) where {T}
    visit!(state.visitor, addr)
    retval::T = get_leaf_node(state.constraints, addr)
    score = logpdf(dist, retval, args...)
    call = CallRecord{T}(score, retval, args)
    state.trace = assoc_primitive_call(state.trace, addr, call)
    state.score += score
    retval
end

function addr(state::GFProjectState, gen::Generator{T,U}, args, addr, args_change) where {T,U}
    visit!(state.visitor, addr)
    if has_internal_node(state.constraints, addr)
        constraints = get_internal_node(state.constraints, addr)
    elseif has_leaf_node(state.constraints, addr)
        lightweight_got_leaf_node_err(addr)
    else
        constraints = EmptyAssignment()
    end
    (trace::U, discard) = project(gen, args, constraints)
    call::CallRecord = get_call_record(trace)
    state.trace = assoc_subtrace(state.trace, addr, trace)
    set_internal_node!(state.discard, addr, discard)
    state.score += call.score
    call.retval::T
end

function splice(state::GFProjectState, gen::GenFunction, args::Tuple)
    exec(gf, state, args)
end

function project(gf::GenFunction, args, constraints)
    state = GFProjectState(constraints, gf.params)
    retval = exec(gf, state, args)
    new_call = CallRecord(state.score, retval, args)
    state.trace.call = new_call
    unvisited_constraints = get_unvisited(state.visitor, state.constraints)
    merge!(state.discard, unvisited_constraints)
    (state.trace, state.discard)
end

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
using ComplexOperations
using LinearAlgebra

# testing complex_matrix_inversion
A = rand(3, 3); B = rand(3, 3)
@time for i = 1:1e6; C, D = complex_matrix_inversion(A, B); end
@time for i = 1:1e6; Z = inv(A + im * B); end
C, D = complex_matrix_inversion(A, B)
Z = inv(A + im * B)
Z == complex.(C, D) # TODO: why there is a small difference?

# testing complex_vector_dot_product
A_r = [1,2,3]; A_i = [1, 20, 30];
B_r = [10,2,30]; B_i = [1, 2, 3];
println(); println("A and B are complex")
@time AB = complex_vector_dot_product(A_r, A_i, B_r, B_i)
@time AB_ = complex_vector_dot_product(hcat(A_r, A_i), hcat(B_r, B_i))

println(); println("A is real, B is complex")
@time AB__ = complex_vector_dot_product(hcat(A_r, zero(A_r)), hcat(B_r, B_i))
@time AB__ = complex_vector_dot_product(A_r, hcat(B_r, B_i))

println(); println("A is complex, B is real")
@time AB_ = complex_vector_dot_product(hcat(A_r, A_i), hcat(B_r, zero(B_r)))
@time AB_ = complex_vector_dot_product(hcat(A_r, A_i), B_r)

# testing complex_vector_dot_product, returns SMatrix
A_r = [1,2,3]; A_i = [1, 20, 30];
B_r = [10,2,30]; B_i = [1, 2, 3];
println(); println("A and B are complex")
@time AB = complex_vector_dot_product_SMatrix(A_r, A_i, B_r, B_i)
@time AB_ = complex_vector_dot_product_SMatrix(hcat(A_r, A_i), hcat(B_r, B_i))

println(); println("A is real, B is complex")
@time AB__ = complex_vector_dot_product_SMatrix(hcat(A_r, zero(A_r)), hcat(B_r, B_i))
@time AB__ = complex_vector_dot_product_SMatrix(A_r, hcat(B_r, B_i))

println(); println("A is complex, B is real")
@time AB_ = complex_vector_dot_product_SMatrix(hcat(A_r, A_i), hcat(B_r, zero(B_r)))
@time AB_ = complex_vector_dot_product_SMatrix(hcat(A_r, A_i), B_r)

# testing complex_vector_cross_product
A_r = [1,2,3]; A_i = [1, 20, 30];
B_r = [10,2,30]; B_i = [1, 2, 3];
println(); println("A and B are complex")
@time AB = complex_vector_cross_product(A_r, A_i, B_r, B_i)
@time AB_ = complex_vector_cross_product(hcat(A_r, A_i), hcat(B_r, B_i))

println(); println("A is real, B is complex")
@time AB__ = complex_vector_cross_product(hcat(A_r, zero(A_r)), hcat(B_r, B_i))
@time AB__ = complex_vector_cross_product(A_r, hcat(B_r, B_i))

println(); println("A is complex, B is real")
@time AB_ = complex_vector_cross_product(hcat(A_r, A_i), hcat(B_r, zero(B_r)))
@time AB_ = complex_vector_cross_product(hcat(A_r, A_i), B_r)

# testing complex_vector_cross_product, returns SMatrix
A_r = [1,2,3]; A_i = [1, 20, 30];
B_r = [10,2,30]; B_i = [1, 2, 3];
println(); println("A and B are complex")
@time AB = complex_vector_cross_product_SMatrix(A_r, A_i, B_r, B_i)
@time AB_ = complex_vector_cross_product_SMatrix(hcat(A_r, A_i), hcat(B_r, B_i))

println(); println("A is real, B is complex")
@time AB__ = complex_vector_cross_product_SMatrix(hcat(A_r, zero(A_r)), hcat(B_r, B_i))
@time AB__ = complex_vector_cross_product_SMatrix(A_r, hcat(B_r, B_i))

println(); println("A is complex, B is real")
@time AB_ = complex_vector_cross_product_SMatrix(hcat(A_r, A_i), hcat(B_r, zero(B_r)))
@time AB_ = complex_vector_cross_product_SMatrix(hcat(A_r, A_i), B_r)

#################### which approach is best for complex matrix inversion?
println(); println("check which matrix inversion method is the best")
function complex_matrix_inversion_1(A::AbstractMatrix{R}, B::AbstractMatrix{R}) where R <: Real
    @fastmath @inbounds inv_A_B = inv(A) * B
    @fastmath @inbounds C = inv(A + B * inv_A_B)
    @fastmath @inbounds D = -inv_A_B * C
    return C, D
end

function complex_matrix_inversion_2(A::AbstractMatrix{R}, B::AbstractMatrix{R}) where R <: Real
    @fastmath @inbounds inv_A_B = inv(A) * B
    return ((inv(A + B * inv_A_B)), -inv_A_B * (inv(A + B * inv_A_B)))
end

function complex_matrix_inversion_3(A::AbstractMatrix{R}, B::AbstractMatrix{R}) where R <: Real
    return ((inv(A + B * inv(A) * B)), -inv(A) * B * (inv(A + B * inv(A) * B)))
end

dim_ = 5
@time for i = 1:1e4; complex_matrix_inversion_1(rand(dim_, dim_), rand(dim_, dim_)); end
@time for i = 1:1e4; complex_matrix_inversion_2(rand(dim_, dim_), rand(dim_, dim_)); end
@time for i = 1:1e4; complex_matrix_inversion_3(rand(dim_, dim_), rand(dim_, dim_)); end


text:
using Test, ImageFeatures, Images, ColorTypes

@testset "Orientations" begin
    img = zeros(20, 20)
    img[6:14, 6:14] .= 1
    orientations = corner_orientations(img)
    orientations_deg = map(rad2deg, orientations)
    expected = [45.0, -45.0, 135.0, -135.0]
    @test all(isapprox(orientations_deg[i], e, atol = 0.0001) for (i, e) in enumerate(expected))
    orientations = corner_orientations(img, Keypoints(imcorner(img)))
    orientations_deg = map(rad2deg, orientations)
    @test all(isapprox(orientations_deg[i], e, atol = 0.0001) for (i, e) in enumerate(expected))
    kernel = ones(5, 5)
    orientations = corner_orientations(img, kernel)
    orientations_deg = map(rad2deg, orientations)
    @test all(isapprox(orientations_deg[i], e, atol = 0.0001) for (i, e) in enumerate(expected))
    orientations = corner_orientations(img, Keypoints(imcorner(img)), kernel)
    orientations_deg = map(rad2deg, orientations)
    @test all(isapprox(orientations_deg[i], e, atol = 0.0001) for (i, e) in enumerate(expected))
    img = zeros(RGB{Float64}, 20, 20)
    img[6:14, 6:14] .= RGB{Float64}(1.0, 1.0, 1.0)
    orientations = corner_orientations(img, Keypoints(imcorner(img)), kernel)
    orientations_deg = map(rad2deg, orientations)
    @test all(isapprox(orientations_deg[i], e) for (i, e) in enumerate(expected))

    img = zeros(20, 20)
    diamond = [ 0.0  0.0  0.0  1.0  0.0  0.0  0.0
                0.0  0.0  1.0  1.0  1.0  0.0  0.0
                0.0  1.0  1.0  1.0  1.0  1.0  0.0
                1.0  1.0  1.0  1.0  1.0  1.0  1.0
                0.0  1.0  1.0  1.0  1.0  1.0  0.0
                0.0  0.0  1.0  1.0  1.0  0.0  0.0
                0.0  0.0  0.0  1.0  0.0  0.0  0.0 ]
    img[8:14, 8:14] = diamond
    orientations = corner_orientations(img)
    orientations_deg = map(rad2deg, orientations)
    expected = [0.0, 90.0, -90.0, -180.0]
    @test all(isapprox(orientations_deg[i], e, atol = 0.0001) for (i, e) in enumerate(expected))
end


text:
include("1.jl")
function temp(r::Robot)
    a = collect(snake1(r))
    if isborder(r,Ost)
        a += collect(movehole12(r,West))

    else
        a += collect(movehole12(r,inverse(West)))
    end
    return(a[1]/a[2])
end

function snake1(r::Robot)
    i = [0,0]
    while !isborder(r,Nord)

        if isborder(r,Ost)
            i += collect(movehole12(r,West))

        else
            i += collect(movehole12(r,inverse(West)))
        end
        if isborder(r,Ost) || isborder(r,West)
            move!(r,Nord)
        end
    end
    return (i)
end
function movehole12(r::Robot,side::HorizonSide)
    t = 0
    n = 0
    while (!isborder(r,side))
        move!(r,side)
        if (ismarker(r))
            t += temperature(r)
            n += 1
        end
    end
    return(t,n)
end

max_stars_repo_name: bubasooo/proggggggggg

max_issues_repo_name: bubasooo/proggggggggg

max_forks_repo_name: bubasooo/proggggggggg


text:
#### thermo states for subdomains

using KernelAbstractions: @print

export new_thermo_state_up,
    new_thermo_state_en,
    recover_thermo_state_all,
    recover_thermo_state_up,
    recover_thermo_state_en

####
#### Interface
####

"""
    new_thermo_state_up(bl, state, aux)

Updraft thermodynamic states given:
 - `bl`, parent `BalanceLaw`
 - `state`, state variables
 - `aux`, auxiliary variables

!!! note
    This method calls saturation adjustment for
    EquilMoist models.
"""
new_thermo_state_up(
    bl::AtmosModel,
    state::Vars,
    aux::Vars,
    ts::ThermodynamicState = recover_thermo_state(bl, state, aux),
) = new_thermo_state_up(bl, moisture_model(bl), state, aux, ts)

"""
    new_thermo_state_en(bl, state, aux)

Environment thermodynamic state given:
 - `bl`, parent `BalanceLaw`
 - `state`, state variables
 - `aux`, auxiliary variables

!!! note
    This method calls saturation adjustment for
    EquilMoist models.
"""
new_thermo_state_en(
    bl::AtmosModel,
    state::Vars,
    aux::Vars,
    ts::ThermodynamicState = recover_thermo_state(bl, state, aux),
) = new_thermo_state_en(bl, moisture_model(bl), state, aux, ts)

"""
    recover_thermo_state_all(bl, state, aux)

Recover NamedTuple of all thermo states

# TODO: Define/call `recover_thermo_state` when it's safely implemented
  (see https://github.com/CliMA/ClimateMachine.jl/issues/1648)
"""
function recover_thermo_state_all(bl, state, aux)
    ts = new_thermo_state(bl, state, aux)
    return (
        gm = ts,
        en = new_thermo_state_en(bl, moisture_model(bl), state, aux, ts),
        up = new_thermo_state_up(bl, moisture_model(bl), state, aux, ts),
    )
end

"""
    recover_thermo_state_up(bl, state, aux, ts = new_thermo_state(bl, state, aux))

Recover the updraft thermodynamic states given:
 - `bl`, parent `BalanceLaw`
 - `state`, state variables
 - `aux`, auxiliary variables

!!! warn
    Right now we are directly calling new_thermo_state_up to avoid
    inconsistent aux states in kernels where the aux states are
    out of sync with the boundary state.

# TODO: Define/call `recover_thermo_state` when it's safely implemented
  (see https://github.com/CliMA/ClimateMachine.jl/issues/1648)
"""
function recover_thermo_state_up(
    bl,
    state,
    aux,
    ts = new_thermo_state(bl, state, aux),
)
    return new_thermo_state_up(bl, moisture_model(bl), state, aux, ts)
end

"""
    recover_thermo_state_en(bl, state, aux, ts = recover_thermo_state(bl, state, aux))

Recover the environment thermodynamic state given:
 - `bl`, parent `BalanceLaw`
 - `state`, state variables
 - `aux`, auxiliary variables

!!! warn
    Right now we are directly calling new_thermo_state_up to avoid
    inconsistent aux states in kernels where the aux states are
    out of sync with the boundary state.

# TODO: Define/call `recover_thermo_state` when it's safely implemented
  (see https://github.com/CliMA/ClimateMachine.jl/issues/1648)
"""
function recover_thermo_state_en(
    bl,
    state,
    aux,
    ts = new_thermo_state(bl, state, aux),
)
    return new_thermo_state_en(bl, moisture_model(bl), state, aux, ts)
end

####
#### Implementation
####

function new_thermo_state_up(
    m::AtmosModel{FT},
    moist::DryModel,
    state::Vars,
    aux::Vars,
    ts::ThermodynamicState,
) where {FT}
    N_up = n_updrafts(turbconv_model(m))
    up = state.turbconv.updraft
    p = air_pressure(ts)
    param_set = parameter_set(m)
    # compute thermo state for updrafts
    ts_up = vuntuple(N_up) do i
        ρa_up = up[i].ρa
        ρaθ_liq_up = up[i].ρaθ_liq
        θ_liq_up =
            fix_void_up(ρa_up, ρaθ_liq_up / ρa_up, liquid_ice_pottemp(ts))
        PhaseDry_pθ(param_set, p, θ_liq_up)
    end
    return ts_up
end

function new_thermo_state_up(
    m::AtmosModel{FT},
    moist::EquilMoist,
    state::Vars,
    aux::Vars,
    ts::ThermodynamicState,
) where {FT}

    N_up = n_updrafts(turbconv_model(m))
    up = state.turbconv.updraft
    p = air_pressure(ts)
    param_set = parameter_set(m)
    # compute thermo state for updrafts
    ts_up = vuntuple(N_up) do i
        ρa_up = up[i].ρa
        ρaθ_liq_up = up[i].ρaθ_liq
        ρaq_tot_up = up[i].ρaq_tot
        θ_liq_up =
            fix_void_up(ρa_up, ρaθ_liq_up / ρa_up, liquid_ice_pottemp(ts))
        q_tot_up = fix_void_up(
            ρa_up,
            ρaq_tot_up / ρa_up,
            total_specific_humidity(ts),
        )
        PhaseEquil_pθq(param_set, p, θ_liq_up, q_tot_up)
    end
    return ts_up
end

function new_thermo_state_en(
    m::AtmosModel,
    moist::DryModel,
    state::Vars,
    aux::Vars,
    ts::ThermodynamicState,
)
    turbconv = turbconv_model(m)
    N_up = n_updrafts(turbconv)
    up = state.turbconv.updraft

    # diagnose environment thermo state
    ρ_inv = 1 / state.ρ
    p = air_pressure(ts)
    θ_liq = liquid_ice_pottemp(ts)
    a_en = environment_area(state, N_up)
    ρaθ_liq_up = vuntuple(N_up) do i
        fix_void_up(up[i].ρa, up[i].ρaθ_liq)
    end
    θ_liq_en = (θ_liq - sum(vuntuple(j -> ρaθ_liq_up[j] * ρ_inv, N_up))) / a_en
    a_min = turbconv.subdomains.a_min
    a_max = turbconv.subdomains.a_max
    param_set = parameter_set(m)
    if !(0 <= θ_liq_en)
        @print("ρaθ_liq_up = ", ρaθ_liq_up[Val(1)], "\n")
        @print("θ_liq = ", θ_liq, "\n")
        @print("θ_liq_en = ", θ_liq_en, "\n")
        error("Environment θ_liq_en out-of-bounds in new_thermo_state_en")
    end
    ts_en = PhaseDry_pθ(param_set, p, θ_liq_en)
    return ts_en
end

function new_thermo_state_en(
    m::AtmosModel,
    moist::EquilMoist,
    state::Vars,
    aux::Vars,
    ts::ThermodynamicState,
)
    turbconv = turbconv_model(m)
    N_up = n_updrafts(turbconv)
    up = state.turbconv.updraft

    # diagnose environment thermo state
    ρ_inv = 1 / state.ρ
    p = air_pressure(ts)
    θ_liq = liquid_ice_pottemp(ts)
    q_tot = total_specific_humidity(ts)
    a_en = environment_area(state, N_up)
    θ_liq_en = (θ_liq - sum(vuntuple(j -> up[j].ρaθ_liq * ρ_inv, N_up))) / a_en
    q_tot_en = (q_tot - sum(vuntuple(j -> up[j].ρaq_tot * ρ_inv, N_up))) / a_en
    a_min = turbconv.subdomains.a_min
    a_max = turbconv.subdomains.a_max
    param_set = parameter_set(m)
    if !(0 <= θ_liq_en)
        @print("θ_liq_en = ", θ_liq_en, "\n")
        error("Environment θ_liq_en out-of-bounds in new_thermo_state_en")
    end
    if !(0 <= q_tot_en <= 1)
        @print("q_tot_en = ", q_tot_en, "\n")
        error("Environment q_tot_en out-of-bounds in new_thermo_state_en")
    end
    ts_en = PhaseEquil_pθq(param_set, p, θ_liq_en, q_tot_en)
    return ts_en
end


text:
module Jlrs
using Base.StackTraces

struct TracedException
    exc
    stacktrace::StackTrace
end

const wakerust = Ref{Ptr{Cvoid}}(C_NULL)

function runasync(func::Function, wakeptr::Ptr{Cvoid}, args...)::Any
    try
        func(args...)
    finally
        ccall(wakerust[], Cvoid, (Ptr{Cvoid},), wakeptr)
    end
end

function asynccall(func::Function, wakeptr::Ptr{Cvoid}, args...)::Task
    @assert wakerust[] != C_NULL "wakerust is null"
    Base.Threads.@spawn runasync(func, wakeptr, args...)
end

function tracingcall(func::Function)::Function
    function wrapper(args...)
        try
            func(args...)
        catch
            for s in stacktrace(catch_backtrace(), true)
                println(stderr, s)
            end

            rethrow()
        end
    end

    wrapper
end

function attachstacktrace(func::Function)::Function
    function wrapper(args...)
        try
            func(args...)
        catch exc
            st::StackTrace = stacktrace(catch_backtrace(), true)
            rethrow(TracedException(exc, st))
        end
    end

    wrapper
end
end

max_stars_repo_name: nikhilmitrax/jlrs

max_issues_repo_name: nikhilmitrax/jlrs

max_forks_repo_name: nikhilmitrax/jlrs


text:
using .RationalFunctionFields

logger = Logging.SimpleLogger(stdout, Logging.Debug)
Logging.global_logger(logger)

@testset "lol" begin
    FF = Nemo.QQ

    R, (x1, x2) = Nemo.PolynomialRing(FF, ["x1", "x2"])

    set = [
          ( x1^2 + x2^2 ) // R(1),
          ( x1 ) // R(1),
    ]

    FF = RationalFunctionField(set)

    println( RationalFunctionFields.contains(FF, x1 // R(1), p=1) )

    println( RationalFunctionFields.contains(FF, x2^2 // R(1), p=0.99) )

    println( RationalFunctionFields.groebner_ideal(set, check=true) )
end


text:
const TimeNS = typeof(time_ns())

struct Event
    name::Symbol
    data::NamedTuple
    location::UUID
    __source__::Union{LineNumberNode,Nothing}
    __module__::Union{Module,Nothing}
    time_ns::TimeNS
end

const EVENT_BLOCKSIZE = Ref(1024)

thread_local_event_records() = BlockLinkedList(Vector{Event}, EVENT_BLOCKSIZE[])

const EventRecordShard = typeof(thread_local_event_records())

# TODO: thread safety
mutable struct EventRecord
    @atomic shards::Vector{EventRecordShard}
end

new_shards() = EventRecordShard[thread_local_event_records() for _ in 1:Threads.nthreads()]

create_records() = init_records!(EventRecord(Union{}[]))
init_records!(records::EventRecord) = empty!(records)

function Base.empty!(records::EventRecord)
    shards = new_shards()
    @atomic records.shards = shards
    return records
end

function record!(
    records::EventRecord,
    name::Symbol,
    data::NamedTuple;
    location::UUID = UUID(0),
    __source__::Union{LineNumberNode,Nothing} = nothing,
    __module__::Union{Module,Nothing} = nothing,
    time_ns::TimeNS = time_ns(),
    # Ignored:
    log::NamedTuple = NamedTuple(),
)
    event = Event(name, data, location, __source__, __module__, time_ns)
    shards = @atomic records.shards
    push!(shards[Threads.threadid()], event)
    return event
end


text:
###
### Mer specific specializations of src/biosequence/indexing.jl
###

@inline _offset(K, i) = 2(K - i)
@inline offset(::Type{T}, i::Integer) where {A,K,T<:AbstractMer{A,K}} = 2(K - i)
@inline offset(x::AbstractMer, i::Integer) = offset(typeof(x), i)

@inline function inbounds_getindex(x::AbstractMer, i::Integer)
    return reinterpret(eltype(x), 0x01 << ((encoded_data(x) >> offset(x, i)) & 0b11))
end

function Base.getindex(mer::T, r::OrdinalRange{<:Integer, <:Integer}) where {A,K,T<:AbstractMer{A,K}}
    parameterless_mer_type = Base.typename(T).wrapper
    new_K = length(r)
    new_T = parameterless_mer_type{A, new_K}
    return new_T(mer[i] for i in r)
end


text:
module DandelionSlack

include("WebApi/json_deserialize.jl")
include("WebApi/helpers.jl")
include("WebApi/newtype.jl")
include("WebApi/auxiliary_types.jl")
include("WebApi/types.jl")
include("WebApi/methods.jl")
include("WebApi/rtm.jl")
include("WebApi/channels.jl")
include("WebApi/requests.jl")
include("WebApi/users.jl")

include("RTM/errors.jl")
include("RTM/events.jl")
include("RTM/register.jl")
include("RTM/misc_events.jl")
include("RTM/messages.jl")
include("RTM/team.jl")
include("RTM/throttling.jl")
include("RTM/rtm.jl")

include("Util/util.jl")

end # module


text:
module CVCalib3d

using LibOpenCV
using CVCore
using Cxx

libopencv_calib3d = LibOpenCV.find_library_e("libopencv_calib3d")
try
    Libdl.dlopen(libopencv_calib3d, Libdl.RTLD_GLOBAL)
catch e
    warn("You might need to set DYLD_LIBRARY_PATH to load dependencies proeprty.")
    rethrow(e)
end

cxx"""
#include <opencv2/calib3d.hpp>
"""

end # module


text:
using QuBase
using QuDynamics

# Quantum Harmonic Oscillator using solvers from QuDynamics

# Constructor for constructing the Hamiltonian, Initial state.
function initialize(n::Int, λ)
    a = lowerop(n)
    hamiltonian = a'*a + λ*(a + a')
    init_state = statevec(1,FiniteBasis(n))
    return hamiltonian, init_state
end

# Initializing the system by selecting the size of the basis and parameter lambda
init = initialize(4, 0.5)

# Time step array
t = 0.:0.1:2*π

# Construct the QuPropagator `P`, by passing the Hamiltonian `H`, initial state `psi`, time step array `t` and Method type
# `Method type` can be QuEuler, QuCrankNicolson, QuKrylov, QuODE45, QuODE78, QuODE23s
Η = init[1]
ψ = init[2]
Ρ = QuPropagator(Η, ψ , t, QuKrylov())

# Evolve the ground state, the next state at the next time step 0.1
η = next(Ρ, start(Ρ))
println(η)

# Alternatively using the iterator methods calculate the evolved states at every point in time step array
for (t, β) in Ρ
    println("$t $β")
end


text:
module NetParam

export NetParameterType
export PairFormat
export NetworkData
export readTouchStone

"""
    NetParameterType

    enum serving as a type specifier for what the data in a NetworkData struct
    represents
"""
@enum NetParameterType begin
    Scattering = 1
    Admittance = 2
    Impedance  = 3
    HybridH    = 4
    HybridG    = 5
end

"""
    PairFormat

    enum serving as descriptive format specifier
"""
@enum PairFormat begin
    MA = 1
    DB = 2
    RI = 3
end

"""
    NetworkData

    Struct containing information describing a microwave network
"""
struct NetworkData
    data::AbstractArray{ComplexF64,3}
    freq::AbstractArray{Float64,1}
    paramType::NetParameterType
    ref_imp::Float64
    noiseData::AbstractArray{Float64,2}
end

"""
    Options
"""
struct DataOptions
    freqmult::Int64
    paramtype::NetParameterType
    format::PairFormat
    impedance::Float64
end

# GENERAL TOUCHSTONE SYNTAX RULES
#  1. Touchstone files are case-insensitive
#  2. Only ASCII chars are allowed, specifically ANSI X3.4-1986 compliant
#     Anything greater than 0x07E or less than 0x20 is not allowed with the
#     Exception of tabs, carriage return, linefeed/cr combo.
#  3. Comments are preceded by '!', can appear on individual line or at end
#     of valid line. No multiline comments
#  4. File extension of '.snp', where 'n' is number of network ports.
#  5. Conventionally, angles are measured in degrees
#
#
#  2-port devices can contain noise parameter data (potentially meaningless on other
#   devices without clarification) -- don't parse for non-two-port devices
#
#  OPTION LINE - always starts with '#', always terminated with a newline
#  # <freq unit> <parameter> <format> R <n>
#  * '#' - marks beginning of option line
#  * frequency unit - unit of frequency, acceptable are (GHz, MHz, KHz, Hz)
#                     again, case insensitive - default is GHz
#  * parameter - type of network parameters
#                   S - scattering - default if otherwise unspecified
#                   Y - admittance
#                   Z - impedance
#                   H - hybrid-h
#                   G - hybrid-g
#  * format - specified format of network param pairs
#               DB - for dB-angle (dB = 20*log10|magnitude|)
#               MA - for magnitude-angle - default if otherwise unspecified
#               RI - for real-imaginary
#             Angles are always in degrees
#             Format does not apply to noise parameters
#  * R n - specifies the reference resistance in ohms, where n is a positive number
#          of ohms - defaults to 50
#
#  DATA LINES - contain network data, only one line for 1-port, 2-port,
#               multi line for all else.
#               Always follows the option line (no need to hunt for option)
#               No more than four pairs of network data per line
#               Individual entries always separated by whitespace
#               Always terminated by newline
#               Always increasing order of frequency
#               one-port - <freq> <N11>
#               two-port - <freq> <N11> <N21> <N12> <N22> in that order, no commas
#               three-port - <freq> <N11> <N12> <N13>
#                            <N21> <N22> <N23>
#                            <N31> <N32> <N33>
#               four-port - <freq> <N11> <N12> <N13> <N14>
#                           <N21> <N22> <N23> <N24>
#                           <N31> <N32> <N33> <N34>
#                           <N41> <N42> <N43> <N44>
#               five+ port - <freq> <N11> <N12> <N13> <N14> ! wrap to next line after four
#                            <N15> <N16>
#                            <N21> <N22> <N23> <N24>
#                            <N25> <N26>
#                            ... and so on
#  NOISE DATA - always follows the network parameter data
#               can be any frequency, as long as lowest noise param is lower than
#               highest network param, for parsing this only applies to two port networks
#               source reflection coefficient and effective noise resistance normalized
#               to same as for network parameters always has format of:
#               <x1> <x2> <x3> <x4> <x5>
#               <x1> - frequency in already specified units
#               <x2> - minimum noise figure in dB
#               <x3> - source reflection coefficient to realize minimimum noise figure (magnitude in MA pair)
#               <x4> - phase in degrees of reflection coefficient (angle in MA pair)
#               <x5> - normalized effective noise resistance
#

"""
    parse_options
        (
        optionstr::String
        )

    Parse out file format properties from the touchstone option line
"""
function parse_options( optionstr::String )::DataOptions

    freqmult = 1
    paramtype = Scattering
    format = MA
    ref_imp = 50

    options = split( optionstr )

    for prop in options
        if "ghz" == prop
            freqmult = 1e9
        elseif "mhz" == prop
            freqmult = 1e6
        elseif "khz" == prop
            freqmult = 1e3
        elseif "hz" == prop
            freqmult = 1
        elseif "s" == prop
            paramtype = Scattering
        elseif "y" == prop
            paramtype = Admittance
        elseif "z" == prop
            paramtype = Impedance
        elseif "h" == prop
            paramtype = HybridH
        elseif "g" == prop
            paramtype = HybridG
        elseif "ma" == prop
            format = MA
        elseif "db" == prop
            format = DB
        elseif "ri" == prop
            format = RI
        elseif "r" == prop
            # do nothing
        elseif tryparse( Float64, prop ) !== nothing
            # is number, assume it's the reference impedance
            ref_imp = parse( Float64, prop )
        end
    end

    return DataOptions( freqmult, paramtype, format, ref_imp )
end

"""
    pair_to_complex
        (
        a::Float64,
        b::Float64,
        format::PairFormat
        )

    Convert a formatted pair into a complex representation
"""
function pair_to_complex( a::Float64, b::Float64, format::PairFormat )
    if MA == format
        r = a * cos( deg2rad( b ) )
        i = a * sin( deg2rad( b ) )
        return Complex{Float64}( r, i )
    elseif DB == format
        # dB = 20*log10|magnitude|)
        mag = 10^( a / 20 )
        r = mag * cos( deg2rad( b ) )
        i = mag * sin( deg2rad( b ) )
        return Complex{Float64}( r, i )
    elseif RI == format
        return Complex{Float64}( a, b )
    end
end

"""
    parse_contained_lines
        (
        line::String,
        portcount::Int32,
        format::PairFormat,
        freqmult::Int64
        )

    for devices with < 3 ports the parameters are all on a single line
    break it down and return the corresponding frequency and matrix
"""
function parse_contained_line( line::String, portcount::Int32,
                               format::PairFormat, freqmult::Int64 )
    contents = split( line )
    freq = parse( Float64, contents[1] ) * freqmult

    pair1v1 = parse( Float64, contents[2] )
    pair1v2 = parse( Float64, contents[3] )

    s11 = pair_to_complex( pair1v1, pair1v2, format )

    if 2 == portcount
        pair2v1 = parse( Float64, contents[4] )
        pair2v2 = parse( Float64, contents[5] )
        pair3v1 = parse( Float64, contents[6] )
        pair3v2 = parse( Float64, contents[7] )
        pair4v1 = parse( Float64, contents[8] )
        pair4v2 = parse( Float64, contents[9] )

        s21 = pair_to_complex( pair2v1, pair2v2, format )
        s12 = pair_to_complex( pair3v1, pair3v2, format )
        s22 = pair_to_complex( pair4v1, pair4v2, format )

        return ( freq, [ s11 s12; s21 s22 ] )
    else
        return ( freq, s11 )
    end
end

"""
    parse_extended_line
        (
        line::String,
        portcount::Int32,
        format::PairFormat,
        freqmult::Float64
        )
"""
function parse_extended_line( line::String, portcount::Int32,
                              format::PairFormat, freqmult::Int64 )
    contents = split( line )
    freq = parse( Float64, contents[1] ) * freqmult

    paircount = portcount * portcount
    row = 1
    col = 1

    mat = zeros( Complex{Float64}, portcount, portcount )

    for idx = 1:2:paircount
        v1 = parse( Float64, contents[idx] )
        v2 = parse( Float64, contents[idx + 1] )

        parameter = pair_to_complex( v1, v2, format )

        mat[row, col] = parameter
        col += 1
        if col > portcount
            col = 1
            row += 1
        end
    end

    return ( freq, mat )
end

"""
    read_normal
        (
        portcount::Int32,
        lines::Vector{String}
        )

    parse the data in lines from a normal touchstone file

    TODO: handle noise data
    TODO: not very efficient due to the means of constructing the matrix
"""
function read_normal( portcount::Int32, lines::Vector{String} )::NetworkData

    # Touchstone defaults
    freqmult = 1
    paramtype = Scattering
    format = MA
    ref_imp = 50

    # empty frequency vector
    frequencies = zeros( 0 )
    # empty port data
    if portcount == 1
        params = Vector{ComplexF64}()
    else
        params = Vector{AbstractMatrix{ComplexF64}}()
    end

    longline = ""
    linecnt = 0
    totallines = Int64( ceil( portcount / 4 ) ) * portcount

    for line in lines
        # eliminate comments
        line = lowercase( split( line, "!" )[begin] )

        if '#' in line
            opts = parse_options( line )
            freqmult = opts.freqmult
            paramtype = opts.paramtype
            format = opts.format
            ref_imp = opts.impedance
        else
            # not an options line, make sure it has content
            if !isempty( line )
                # TODO: optimize this
                if portcount < 3
                    # for one and two, all on one line
                    (freq, mat) = parse_contained_line( line, portcount,
                                                        format, freqmult )

                    push!( frequencies, freq )
                    push!( params, mat )
                elseif portcount >= 3
                    # three and four, they get wrapped
                    # 5 and up, wrapped and wrapped
                    linecnt += 1
                    longline = string( longline, ' ', line )

                    if linecnt == totallines
                        (freq, mat) = parse_extended_line( longline, portcount,
                                                           format, freqmult )
                        push!( frequencies, freq )
                        push!( params, mat )
                        longline = ""
                        linecnt = 0
                    end
                end
            end
        end
    end

    # Convert from vector of 2d matrices to fixed 3d matrix
    # TODO: clean things up for the 1x1 matrix
    plen = length( params )
    parameters = zeros( ComplexF64, portcount, portcount, plen )

    if 1 == portcount
        for idx = 1:plen
            parameters[:,:,idx] = [ params[idx] ]
        end
    else
        for idx = 1:plen
            parameters[:,:, idx] = params[idx]
        end
    end

    return NetworkData( parameters, frequencies, paramtype, ref_imp, zeros( 0, 0 ) )
end

"""
    read_mixer
        (
        portcount::Int32,
        lines::Vector{String}
        )

    parse data data in lines for touchstone mixer file
"""
function read_mixer( portcount::Int32, lines::Vector{String} )::NetworkData
    # TODO: implement this function

    error( "mixer files not yet supported" )
    # # Touchstone defaults
    # freqmult = 1
    # paramtype = Scattering
    # format = MA
    # ref_imp = 50
    # csvmode = false
    #
    # for line in lines
    #
    #     if 'CSV' in line
    #         csvmode = true
    #     end
    #
    #     if csvmode
    #     else
    #         contents = split( line )
    #     end
    #
    #
    # end
end

"""
    read_touchstone
        (
        path::String
        )

    Open and parse a file conforming to the Touchstone file format and
    return a corresponding NetworkData object.

    TODO: Check ANSI character encoding compliance
          Error checking is for chumps
    TODO: don't try to handle mixer data in this function - create separate
          datatype to handle mixer data
"""
function read_touchstone( path::String )::NetworkData
    touchstone_file = open( path, "r" )

    if isfile( touchstone_file )
        extension = split( path, '.' )[end]
        # Offset frequency / mixer touchstones files designated
        # by appended x typically: '.sNpx'
        ismixer = ( 'x' in extension )

        if ismixer
            error( "mixer files not supported, use read_mixer() instead" )
        else
            lines = readlines( touchstone_file )
            portcount = parse( Int32, extension[2:end-1] )
            return read_normal( portcount, lines )
        end
    end
end

"""
    write_touchstone
        (
        path::String,
        network_data::NetworkData
        )
"""
function write_touchstone( path::String, network_data::NetworkData )
    # TODO: implement
    error( "write_touchstone not yet implemented" )
end

end # module


text:
"""
    InvertedIndex(idx)
    Not(idx)

Construct an inverted index, selecting all indices not in the passed `idx`.

Upon indexing into an array, the `InvertedIndex` behaves like a 1-dimensional
collection of the indices of the array that are not in `idx`. Bounds are
checked to ensure that all indices in `idx` are within the bounds of the array
— even though they are skipped. The `InvertedIndex` behaves like a
1-dimensional collection of its inverted indices. If `idx` spans multiple
dimensions (like a multidimensional logical mask or `CartesianIndex`), then the
inverted index will similarly span multiple dimensions.

Copied and trimmed down from InvertedIndices.jl repository
"""
struct InvertedIndex{T}
    skip::T
end
const Not = InvertedIndex

import Base.trues
# This is a little tricky because trues isn't indices-aware, but we also don't
# want to use fill(true) in the 1-indexed case since we want to favor BitArrays
@inline trues(tup::Tuple{Vararg{Base.OneTo}}) = Base.trues(map(Base.unsafe_length, tup))
@inline trues(tup::Tuple{Vararg{Base.AbstractUnitRange}}) = fill(true, tup)

@inline function Base.to_indices(A, inds, I::Tuple{InvertedIndex, Vararg{Any}})
    v = trues(spanned_indices(inds, I))
    v[I[1].skip] = false
    to_indices(A, inds, (v, Base.tail(I)...))
end

# Determining the indices that the InvertedIndex spans is tricky due to partial
# linear indexing. Lean on `Base.uncolon` until the deprecation goes through.
@inline spanned_indices(inds, I::Tuple{InvertedIndex,Vararg{Any}}) = (Base.uncolon(inds, (:, Base.tail(I)...)).indices,)


text:
module Bencode

export bencode, bdecode

function bencode(b::Vector{UInt8})::Vector{UInt8}
    blen = length(b)
    blenbytestring = Vector{UInt8}(string(blen))
    [blenbytestring; UInt8(':'); b]
end

function bencode(s::AbstractString)::Vector{UInt8}
    bencode(Vector{UInt8}(s))
end

function bencode(n::Integer)::Vector{UInt8}
    Vector{UInt8}("i" * string(n) * "e")
end

function bencode(l::AbstractVector)::Vector{UInt8}
    if isempty(l)
        Vector{UInt8}("le")
    else
        listcontent = mapreduce(bencode, vcat, l)
        [UInt8('l'); listcontent; UInt8('e')]
    end
end

function bencode(d::Dict)::Vector{UInt8}
    if isempty(d)
        return Vector{UInt8}("de")
    else
        ks = collect(keys(d))
        keyisless = (a, b) -> isless(string(a), string(b))
        sortedkeys = sort(ks, lt=keyisless)
        encodekeyvalue(key) = [bencode(key); bencode(d[key])]
        dictcontent = mapreduce(encodekeyvalue, vcat, sortedkeys)
        [UInt8('d'); dictcontent; UInt8('e')]
    end
end

function bdecode(s::AbstractString; bytestostr::Bool=true, retnbytesread::Bool=false)
    bdecode(Vector{UInt8}(s); bytestostr=bytestostr, retnbytesread=retnbytesread)
end

function bdecode(data::AbstractVector{UInt8}; bytestostr::Bool=true, retnbytesread::Bool=false)
    val, nread = if isdigit(Char(data[1]))
        bytestostr ? bdecodestr(data) : bdecodebytes(data)
    elseif Char(data[1]) == 'i'
        bdecodeint(data)
    elseif Char(data[1]) == 'l'
        bdecodelist(data, bytestostr)
    elseif Char(data[1]) == 'd'
        bdecodedict(data, bytestostr)
    end
    retnbytesread ? (val, nread) : val
end

function bdecodeint(data::AbstractVector{UInt8})::Tuple{Int, Int}
    indexofend = findfirst(isequal(UInt8('e')), data)
    n = parse(Int, String(data[2:indexofend - 1]))
    n, indexofend
end

function bdecodebytes(data::AbstractVector{UInt8})::Tuple{Vector{UInt8}, Int}
    colonindex = findfirst(isequal(UInt8(':')), data)
    strlen = parse(Int, String(copy(data[1:colonindex - 1])))
    str = data[colonindex + 1:colonindex + strlen]
    str, colonindex + strlen
end

function bdecodestr(data::AbstractVector{UInt8})::Tuple{String, Int}
    b, nread = bdecodebytes(data)
    String(copy(b)), nread
end

function bdecodelist(data::AbstractVector{UInt8}, bytestostr::Bool)::Tuple{Vector, Int}
    list = []
    totalread = 1
    data = data[2:end]
    while Char(data[1]) != 'e'
        element, nread = bdecode(data; bytestostr=bytestostr, retnbytesread=true)
        push!(list, element)
        totalread += nread
        data =  data[nread + 1:end]
    end
    list, totalread + 1
end

function bdecodedict(data::AbstractVector{UInt8}, bytestostr::Bool)::Tuple{Dict, Int}
    dict = Dict()
    totalread = 1
    data = data[2:end]
    while Char(data[1]) != 'e'
        key, nreadkey = bdecode(data;
                bytestostr=true, retnbytesread=true) # bytes are not allowed in keys
        data = data[nreadkey + 1:end]
        val, nreadval = bdecode(data; bytestostr=bytestostr, retnbytesread=true)
        dict[key] = val
        data = data[nreadval + 1:end]
        totalread += nreadkey + nreadval
    end
    dict, totalread + 1
end

end # module


text:
import .ArcStd: State, expand_gold, expand_pred, toheads

type DepParser
end

const num_feats = 1 << 26
const weights = [fill(0.0,num_feats) for i=1:3]

function tohash(data::Vector{Int})
    hash = 1
    for x in data
        hash += x
        hash += (hash << 10)
        hash $= (hash >> 6)
    end
    hash += (hash << 3)
    hash $= (hash >> 11)
    hash += (hash << 15)
    abs(hash) % num_feats + 1
end

function scorefun(s::State, acts::Vector{Int})
    feats = Int[]
    tokens = s.tokens
    s1 = tokens[s.top]
    s2 = s.left == nothing ? nulltoken : tokens[s.left.top]
    b1 = s.right <= length(tokens) ? tokens[s.right] : nulltoken
    b2 = s.right+1 <= length(tokens) ? tokens[s.right+1] : nulltoken
    xs = Any[
    Int[1, s1.formid],
    Int[2, s1.catid],
    Int[3, s1.formid, s2.catid],
    Int[4, s2.formid],
    Int[5, s2.catid],
    Int[6, s2.formid, s2.catid],
    Int[7, b1.formid],
    Int[8, b1.catid],
    Int[9, b1.formid, b1.catid],
    Int[10, s1.formid, s2.formid],
    Int[11, s1.catid, s2.catid],
    Int[12, s1.catid, b1.catid],
    Int[13, s1.formid, s1.catid, s2.catid],
    Int[14, s1.catid, s2.formid, s2.catid],
    Int[15, s1.formid, s2.formid, s2.catid],
    Int[16, s1.formid, s1.catid, s2.formid],
    Int[17, s1.formid, s1.catid, s2.formid, s2.catid],
    Int[18, s1.catid, b1.catid, b2.catid],
    Int[19, s2.catid, s1.catid, b1.catid],
    Int[20, s1.formid, b1.catid, b2.catid],
    Int[21, s2.catid, s1.formid, b1.catid],
    ]
    for x in xs
        push!(feats, tohash(x))
    end
    s.feats = feats

    scores = Float64[]
    for a in acts
        s = mapreduce(f -> weights[a][f], +, feats)
        push!(scores, s)
    end
    scores
end

function train_gold(s::State)
    feats::Vector{Int} = s.prev.feats
    for f in feats
        weights[s.prevact][f] += 1.0
    end
end

function train_pred(s::State)
    feats::Vector{Int} = s.prev.feats
    for f in feats
        weights[s.prevact][f] -= 1.0
    end
end

function decode(p::DepParser, tokens::Vector{Token})
end

function train(p::DepParser, traindata, testdata)
    for epoch = 1:20
        println("epoch: $(epoch)")
        loss = 0.0
        for i in randperm(length(traindata))
            #for x in traindata
            x = traindata[i]
            sy = State(x, scorefun)
            sz = State(x, scorefun)
            y = beamsearch(sy, 1, expand_gold)[end][1]
            z = beamsearch(sz, 1, expand_pred)[end][1]
            loss += z.score - y.score
            max_violation!(y, z, train_gold, train_pred)
        end
        println("loss: $(loss)")
        golds, preds = Int[], Int[]
        for x in testdata
            s = State(x, scorefun)
            z = beamsearch(s, 1, expand_pred)[end][1]
            append!(golds, map(t -> t.headid, x))
            append!(preds, toheads(z))
        end
        acc = accuracy(golds, preds)
        println("test acc: $(acc)")
        println("")
    end
end

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

const CURRENT_SAMPLE = -1
const PROPOSED_SAMPLE = -2
const INVALID_SAMPLE = 0
const ACCEPTED_SAMPLE = 1
const REJECTED_SAMPLE = 2

abstract type SampleID end

struct MCMCSampleID <: SampleID
    chainid::Int32
    chaincycle::Int32
    stepno::Int64
    sampletype::Int32
end

const MCMCSampleIDVector{TV<:AbstractVector{<:Int32},UV<:AbstractVector{<:Int64}} = StructArray{
    MCMCSampleID,
    1,
    NamedTuple{(:chainid, :chaincycle, :stepno, :sampletype), Tuple{TV,TV,UV,UV}},
    Int
}

function MCMCSampleIDVector(contents::Tuple{TV,TV,UV,UV}) where {TV<:AbstractVector{<:Int32},UV<:AbstractVector{<:Int64}}
    StructArray{MCMCSampleID}(contents)::MCMCSampleIDVector{TV,UV}
end

MCMCSampleIDVector(::UndefInitializer, len::Integer) = MCMCSampleIDVector((
    Vector{Int32}(undef, len), Vector{Int32}(undef, len),
    Vector{Int64}(undef, len), Vector{Int64}(undef, len)
))

MCMCSampleIDVector() = MCMCSampleIDVector(undef, 0)

_create_undef_vector(::Type{MCMCSampleID}, len::Integer) = MCMCSampleIDVector(undef, len)

# Specialize comparison, currently StructArray seems fall back to `(==)(A::AbstractArray, B::AbstractArray)`
import Base.==
function(==)(A::MCMCSampleIDVector, B::MCMCSampleIDVector)
    A.chainid == B.chainid &&
    A.chaincycle == B.chaincycle &&
    A.stepno == B.stepno &&
    A.sampletype == B.sampletype
end

function Base.merge!(X::MCMCSampleIDVector, Xs::MCMCSampleIDVector...)
    for Y in Xs
        append!(X, Y)
    end
    X
end

Base.merge(X::MCMCSampleIDVector, Xs::MCMCSampleIDVector...) = merge!(deepcopy(X), Xs...)


text:
module SolarModel

using ..vSmartMOM               # For locating default solar T
using DocStringExtensions       # For simplifying docstring
using DelimitedFiles            # For easily reading in solar spectrum
using Interpolations            # For interpolating solar spectrum
using Pkg.Artifacts             # For default solar spectrum

"""
    $(FUNCTIONNAME)(T::Real, ν_grid::Vector)

Produce the black-body planck spectrum (mW/m²-sr-cm⁻¹), given the temperature (K)
and calculation grid (ν in cm⁻¹)

"""
function planck_spectrum_wn(T::Real, ν_grid::Vector)

    c1 = 1.1910427 * 10^(-5)    # mW/m²-sr-cm⁻¹
    c2 = 1.4387752              # K⋅cm

    # L(ν, T) = c1⋅ν³/(exp(c2⋅ν/T) - 1)
    radiance = c1 .* (ν_grid.^3) ./ (exp.(c2 * ν_grid / T) .- 1)

    return radiance
end

"""
    $(FUNCTIONNAME)(T::Real, λ_grid::Vector)

Produce the black-body planck spectrum (W/m²-sr-μm), given the temperature (K)
and calculation grid (λ in μm)

"""
function planck_spectrum_wl(T::Real, λ_grid::Vector)

    c1 = 1.1910427 * 10^8    # W/m²-sr-μm
    c2 = 1.4387752 * 10^4    # K⋅μm

    # L(ν, T) = c1⋅ν³/(exp(c2⋅ν/T) - 1)
    radiance = c1 ./ (λ_grid.^5 .* (exp.(c2 ./ (λ_grid * T)) .- 1))

    return radiance
end

# W/m²-sr-μm to Ph/s-m²-sr-um
# λ_grid in micron
function watts_to_photons(λ_grid::Vector, radiance::Vector)

    h = 6.62607015e-34 # J⋅Hz−1
    c = 299792458 # m/s

    E_per_λ = h * c ./ (λ_grid / 1e6)
    photons = radiance ./ E_per_λ

    return photons
end

"""
    $(FUNCTIONNAME)(T::Real; stride_length::Integer = 100)

Produce the black-body planck spectrum (mW/m²-sr-cm⁻¹), given the temperature (K).
Use a unit calculation grid and check for convergence every `stride_length` cm⁻¹ until the
spectrum dies off.

"""
function planck_spectrum_wn(T::Real; stride_length::Integer = 100)

    # νs, starting with ν0 = 1.0 cm⁻¹
    νs = [1.0]

    # radiances corresponding with νs
    radiances = planck_spectrum_wn(T, νs)

    # Loop until convergence
    while true

        # Add the next ν
        νs = vcat(νs, collect(νs[end] + 1 : νs[end] + stride_length))

        # Compute the next radiance
        radiances = vcat(radiances, planck_spectrum_wn(T, νs[(end - stride_length + 1) : end]))

        # Exit if spectrum has died off
        (radiances[end] < radiances[1]) && break

    end

    return [νs[1:(end-1)] radiances[1:(end-1)]]
end

"""
    $(FUNCTIONNAME)(solar, ν_grid)

Interpolate a solar linelist to the ν_grid
"""
function itp_solar_to_ν_grid(solar, ν_grid)

    solar_idx_start = maximum((argmin(abs.(solar[:, 1] .- minimum(ν_grid))) - 10, 1))
    solar_idx_end   = minimum((argmin(abs.(solar[:, 1] .- maximum(ν_grid))) + 10, length(solar[:,1])))

    solar_subset = solar[solar_idx_start:solar_idx_end, :]

    itp = LinearInterpolation(solar_subset[:, 1],
                              solar_subset[:, 2])

    return itp.(ν_grid)
end

"""
    $(FUNCTIONNAME)(file_name::String)

Get the solar transmission from the specified file
"""
solar_transmission_from_file(file_name::String) = readdlm(file_name)

"""
    $(FUNCTIONNAME)(file_name::String, ν_grid::Union{AbstractRange{<:Real}, AbstractArray})

Get the solar transmission from the specified file, and interpolate to wavenumber grid
"""
function solar_transmission_from_file(file_name::String,
                                      ν_grid::Union{AbstractRange{<:Real}, AbstractArray})

    solar = solar_transmission_from_file(file_name)
    return itp_solar_to_ν_grid(solar, ν_grid)
end

"""
    $(FUNCTIONNAME)(ν_grid::Union{AbstractRange{<:Real}, AbstractArray} = 600.0:0.01:26316.0)

Get the default solar transmission and interpolate to wavenumber grid (entire grid if not specified)
"""
function default_solar_transmission(ν_grid::Union{AbstractRange{<:Real}, AbstractArray} = 600.0:0.01:26316.0)

    @info "Using line-list from:\nToon, G. C., Solar line list for GGG2014, TCCON data archive, hosted by the Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, Oak Ridge, Tennessee, U.S.A., doi:10. 14291/tccon.ggg2014.solar.R0/1221658, 2014."
    @info "Found at: https://mark4sun.jpl.nasa.gov/toon/solar/solar_spectrum.html"

    filename = joinpath(dirname(pathof(RadiativeTransfer)), "SolarModel", "solar.out")
    !isfile(filename) && download("http://web.gps.caltech.edu/~cfranken/hitran_2016/solar_merged_20160127_600_26316_100.out", filename)

    return hcat(ν_grid, solar_transmission_from_file(filename, ν_grid))
end

"""
    $(FUNCTIONNAME)(ν_grid::Union{AbstractRange{<:Real}, AbstractArray} = 600.0:0.01:26316.0)

Get the default solar spectrum and interpolate to wavenumber grid (entire grid if not specified)
"""
function default_solar_spectrum_at_earth(ν_grid::Union{AbstractRange{<:Real}, AbstractArray} = 600.0:0.01:26316.0)

    T = 5777 # K
    λ_grid = reverse(1e4 ./ ν_grid) # Wavenumber to micron
    black_body = reverse(SolarModel.watts_to_photons(λ_grid, planck_spectrum_wl(T, λ_grid) * 2.1629e-05 * pi))
    solar_transmission = default_solar_transmission(ν_grid)[:,2]

    return hcat(ν_grid, black_body .* solar_transmission)
end

export planck_spectrum_wn, planck_spectrum_wl, solar_transmission_from_file, default_solar_transmission

end