text:
module Test_model_matching

using DescriptorSystems
using LinearAlgebra
using Polynomials
using Random
using Test

println("Test_model_matching")

Random.seed!(2123)

@testset "Model_matching " begin

@testset "grasol" begin

# continuous-time, both G and F unstable
n1 = 3; n2 = 2; m = 2; p = 3; mf = 1;
G = rss(n1,p,m); F = rss(n2,p,mf);
@time X, info = grasol(G, F; reltol=0.0001, atol = 1.e-7); info
@test abs(glinfnorm(G*X-F,offset=1.e-13)[1] - info.mindist) < 0.01 &&  (!info.nonstandard && isstable(X,offset=1.e-13))

@time X, info = grasol(G, F; nehari = true, atol = 1.e-7); info
@test abs(glinfnorm(G*X-F)[1] - info.mindist) < 0.01 &&  (!info.nonstandard && isstable(X))

# discrete-time, both G and F unstable
n1 = 3; n2 = 2; m = 2; p = 3; mf = 1;
G = rss(n1,p,m,disc=true); F = rss(n2,p,mf,disc=true);
@time X, info = grasol(G, F; reltol=0.0001, atol = 1.e-9); info
@test abs(glinfnorm(G*X-F,offset=1.e-13)[1] - info.mindist) < 0.01 &&  (!info.nonstandard && isstable(X,offset=1.e-13))

@time X, info = grasol(G, F; nehari = true, atol = 1.e-7); info
@test abs(glinfnorm(G*X-F)[1] - info.mindist) < 0.01 &&  (!info.nonstandard && isstable(X))

# continuous-time, both G and F unstable and complex
n1 = 3; n2 = 2; m = 2; p = 3; mf = 1;
G = rss(n1,p,m,T = Complex{Float64}); F = rss(n2,p,mf,T = Complex{Float64});
@time X, info = grasol(G, F; reltol=0.0001, atol = 1.e-9); info
@test abs(glinfnorm(G*X-F,offset=1.e-13)[1] - info.mindist) < 0.01 &&  (!info.nonstandard && isstable(X,offset=1.e-13))

@time X, info = grasol(G, F; nehari = true, atol = 1.e-7); info
@test abs(glinfnorm(G*X-F)[1] - info.mindist) < 0.01 &&  (!info.nonstandard && isstable(X))

# discrete-time, both G and F unstable and complex
n1 = 3; n2 = 2; m = 2; p = 3; mf = 1;
G = rdss(n1,p,m,T = Complex{Float64},disc=true); F = rss(n2,p,mf,T = Complex{Float64},disc=true);
@time X, info = grasol(G, F; reltol=0.001, offset=1.e-13, atol = 1.e-13); info
@test abs(glinfnorm(G*X-F,offset=1.e-13,atol=1.e-10)[1] - info.mindist) < 0.01 &&  (!info.nonstandard && isstable(X,offset=1.e-13))

@time X, info = grasol(G, F; nehari = true, atol = 1.e-7); info
@test abs(glinfnorm(G*X-F)[1] - info.mindist) < 0.01 &&  (!info.nonstandard && isstable(X))

# both G and F stable
n1 = 3; n2 = 2; m = 2; p = 3; mf = 1;
G = rss(n1,p,m,stable=true); F = rss(n2,p,mf,stable=true);
@time X, info = grasol(G, F; atol = 1.e-8); info
@test abs(glinfnorm(G*X-F,offset=1.e-13)[1] - info.mindist) < 0.01 &&  (!info.nonstandard && isstable(X,offset=1.e-13))

@time X, info = grasol(G, F; nehari = true, atol = 1.e-7); info
@test abs(glinfnorm(G*X-F)[1] - info.mindist) < 0.01 &&  (!info.nonstandard && isstable(X))

# exact stable solution exists, no free poles
n1 = 3; n2 = 2; m = 2; p = 3; mf = 1;
#G = rss(n1,p,m,stable=true); X0 = rss(n2,pf,p,stable=true);  F = X0*G;
G = rss(n1,p,m); X0 = rss(n2,m,mf,stable=true);  F = G*X0;
@time X, info = grasol(G, F, offset = 1.e-13, atol = 1.e-7); info
@test abs(glinfnorm(G*X-F)[1] - info.mindist) < 1.e-5  &&  (!info.nonstandard && isstable(X))

# exact solution exists, free poles exist
n1 = 3; n2 = 2; m = 2; p = 3; mf = 1;
G = rss(n1,p,m,stable=true); X0 = rss(n2,m,mf,stable=true);  F = G*X0;
@time X, info = grasol(G, F, atol = 1.e-7, poles = [-2, -3], sdeg = -1, mindeg = true); info
@test abs(glinfnorm(G*X-F)[1] - info.mindist) < 1.e-5 &&  (!info.nonstandard && isstable(X)) &&
      order(X) == order(X0)

## Glover & Packard (SCL,2017)
s = rtf('s');
epsi = -1;
g = (s-1)*(s-epsi)/(s+1)/(s+2); f = 1/(s+1);
sysg = gir(dss(g)); sysf = dss(f);
@time sysx, info = grasol(sysg, sysf, mindeg = true); info
@test glinfnorm(sysx*sysg-sysf)[1] ≈ info.mindist && info.nrank == 1 && info.nonstandard == false &&
      gpole(sysx) ≈ [-1]

s = rtf('s');
epsi = 1;
g = (s-1)*(s-epsi)/(s+1)/(s+2); f = 1/(s+1);
sysg = gir(dss(g)); sysf = dss(f);
@time sysx, info = grasol(sysg, sysf, atol=1.e-7); info
# check solution
γ0 = 0.25*(1+sqrt(1+8/(1+epsi)))
α = 1 + 2(1 + epsi)*γ0
@test glinfnorm(sysx*sysg-sysf)[1] ≈ info.mindist && info.mindist ≈ γ0 &&
       gpole(sysx) ≈ [-α]

s = rtf('s');
epsi = 1;
gf = [(s-1)*(s-epsi)/(s+1)/(s+2) 1/(s+1)];
sysgf = gir(dss(gf));
@time sysx, info = grasol(sysgf, 1, atol=1.e-7); info
# check solution
γ0 = 0.25*(1+sqrt(1+8/(1+epsi)))
α = 1 + 2(1 + epsi)*γ0
@test glinfnorm(sysx*sysgf[:,1]-sysgf[:,2])[1] ≈ info.mindist && info.mindist ≈ γ0 &&
       gpole(sysx) ≈ [-α]

# improper L∞ solution for a nonstandard problem
s = rtf('s');
sysg = dss([1/(s+1); 1/(s+2)]); sysf = dss([ 1/(s+1); 0]);
@time sysx, info = grasol(sysg, sysf, atol=1.e-7); info
@test abs(glinfnorm(sysg*sysx-sysf)[1] - info.mindist) < 0.0001 && info.nonstandard ≈ true &&
       gpole(sysx) ≈ [-1.5811388300841827, Inf]

# proper L2 solution for a nonstandard problem
@time sysx, info = grasol(sysg, sysf, atol=1.e-7, L2sol = true); info
@test gl2norm(sysg*sysx-sysf) ≈ info.mindist && info.nonstandard ≈ true &&
       gpole(sysx) ≈ [-1.5811388300841827]

G = rss(0,3,2);
@time X, info = grasol(G, G, atol = 1.e-7); info
@test gl2norm(G*X-G,atolinf=1.e-7) < 1.e-7 && X.D ≈ eye(2)

G = rss(0,3,2);
@time X, info = grasol(G, G, atol = 1.e-7, L2sol = true); info
@test gl2norm(G*X-G,atolinf=1.e-7) < 1.e-7 && X.D ≈ eye(2)

G = rss(3,3,2,stable=true);
@time X, info = grasol(G, G, atol = 1.e-7); info
@test glinfnorm(G*X-G)[1] < 1.e-7 && X.D ≈ eye(2)

G = rss(3,3,2,stable=true);
@time X, info = grasol(G, G, atol = 1.e-7, L2sol = true); info
@test gl2norm(G*X-G,atolinf=1.e-7) < 1.e-7 && X.D ≈ eye(2)

##  Francis (1987) Example 1, p. 112
s = rtf('s'); # define the complex variable s
# enter W(s), G(s) and F(s)
W = (s+1)/(10*s+1); # weighting function
#G = W.*[ -(s-1)/(s^2+s+1); (s^2-2*s)/(s^2+s+1)];  # erroneous
G = [ -(s-1)/(s^2+s+1)*W; (s^2-2*s)/(s^2+s+1)*W];
F = [ W; 0 ];
# solve G(s)*X(s) = F(s) for the least order solution
sysg = dss(G); sysf = dss(F);
X, info = grasol(sysg,sysf,atol=1.e-7,mindeg = true); info
@test glinfnorm(sysg*X-sysf)[1] ≈ info.mindist
# compute the suboptimal solution for gamma = 0.2729
Xsub, info = grasol(sysg,sysf,0.2729,atol=1.e-7,mindeg = true); info
@test glinfnorm(sysg*Xsub-sysf)[1] ≈ info.mindist

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

@time sysx, info = grasol(sysg, sysf); info
@test abs(glinfnorm(sysg*sysx-sysf)[1] ≈ info.mindist) < 1.e-7 && isstable(sysx)

@time sysx, info = grasol(sysg, sysf, mindeg = true); info
@test abs(glinfnorm(sysg*sysx-sysf)[1] ≈ info.mindist) < 1.e-7 && isstable(sysx)

evals=[-1, -2, -3];
@time sysx, info = grasol(sysg, sysf, poles = evals); info
@test abs(glinfnorm(sysg*sysx-sysf)[1] ≈ info.mindist) < 1.e-7 && sort(gpole(sysx)) ≈ sort(evals)

sdeg = -1;
@time sysx, info = grasol(sysg, sysf, sdeg = sdeg); info
@test abs(glinfnorm(sysg*sysx-sysf)[1] ≈ info.mindist) < 1.e-7 &&
      all(real(gpole(sysx)) .< sdeg*(1-eps()^(1/info.nr)))

@time sysx, info = glasol(gdual(sysg), sysf); info
@test abs(glinfnorm(sysx*gdual(sysg)-sysf,atol=1.e-7)[1] - info.mindist) < 1.e-7 && isstable(sysx)

end #grasol

@testset "glasol" begin

# continuous-time, both G and F unstable
n1 = 3; n2 = 2; m = 3; p = 2; pf = 1;
G = rss(n1,p,m); F = rss(n2,pf,m);
@time X, info = glasol(G, F; reltol=0.001, atol = 1.e-7); info
@test abs(glinfnorm(X*G-F,offset=1.e-13)[1] - info.mindist) < 0.01 &&  (!info.nonstandard && isstable(X,offset=1.e-13))

@time X, info = glasol(G, F; nehari = true, atol = 1.e-7); info
@test abs(glinfnorm(X*G-F)[1] - info.mindist) < 0.01 &&  (!info.nonstandard && isstable(X))

# discrete-time, both G and F unstable #fails
n1 = 3; n2 = 2; m = 3; p = 2; pf = 1;
G = rss(n1,p,m,disc=true); F = rss(n2,pf,m,disc=true);
@time X, info = glasol(G, F; reltol=0.0001, atol = 1.e-9); info
@test abs(glinfnorm(X*G-F,offset=1.e-13)[1] - info.mindist) < 0.01 &&  (!info.nonstandard && isstable(X,offset=1.e-13))

@time X, info = glasol(G, F; nehari = true, atol = 1.e-7); info
@test abs(glinfnorm(X*G-F)[1] - info.mindist) < 0.01 &&  (!info.nonstandard && isstable(X))

# continuous-time, both G and F unstable and complex
n1 = 3; n2 = 2; m = 3; p = 2; pf = 1;
G = rss(n1,p,m,T = Complex{Float64}); F = rss(n2,pf,m,T = Complex{Float64});
@time X, info = glasol(G, F; reltol=0.0001, atol = 1.e-9); info
@test abs(glinfnorm(X*G-F,offset=1.e-13)[1] - info.mindist) < 0.01 &&  (!info.nonstandard && isstable(X,offset=1.e-13))

@time X, info = glasol(G, F; nehari = true, atol = 1.e-7); info
@test abs(glinfnorm(X*G-F)[1] - info.mindist) < 0.01 &&  (!info.nonstandard && isstable(X))

# # discrete-time, both G and F unstable and complex
# n1 = 3; n2 = 2; m = 3; p = 2; pf = 1;
# G = rdss(n1,p,m,T = Complex{Float64},disc=true); F = rss(n2,pf,m,T = Complex{Float64},disc=true);
# @time X, info = glasol(G, F; reltol=0.0001, atol = 1.e-9); info
# @test abs(glinfnorm(X*G-F,offset=1.e-10)[1] - info.mindist) < 0.01 &&  (!info.nonstandard && isstable(X,offset=1.e-13))

# @time X, info = glasol(G, F; nehari = true, atol = 1.e-7); info
# @test abs(glinfnorm(X*G-F)[1] - info.mindist) < 0.01 &&  (!info.nonstandard && isstable(X))

# both G and F stable
n1 = 3; n2 = 2; m = 3; p = 2; pf = 1;
G = rss(n1,p,m,stable=true); F = rss(n2,pf,m,stable=true);
@time X, info = glasol(G, F; atol = 1.e-8); info
@test abs(glinfnorm(X*G-F,offset=1.e-13)[1] - info.mindist) < 0.01 &&  (!info.nonstandard && isstable(X,offset=1.e-13))

@time X, info = glasol(G, F; nehari = true, atol = 1.e-7); info
@test abs(glinfnorm(X*G-F)[1] - info.mindist) < 0.01 &&  (!info.nonstandard && isstable(X))

# exact solution exists, no free poles
n1 = 3; n2 = 2; m = 3; p = 2; pf = 1;
G = rss(n1,p,m); X0 = rss(n2,pf,p,stable=true);  F = X0*G;
@time X, info = glasol(G, F, offset = 1.e-13, atol = 1.e-7); info
@test abs(glinfnorm(X*G-F)[1] - info.mindist) < 1.e-5 &&  (!info.nonstandard && isstable(X))

# exact solution exists, free poles exist
n1 = 3; n2 = 2; m = 2; p = 3; pf = 1;
G = rss(n1,p,m,stable=true); X0 = rss(n2,pf,p,stable=true);  F = X0*G;
@time X, info = glasol(G, F, atol = 1.e-7, poles = [-2, -3], sdeg = -1); info
@test abs(glinfnorm(X*G-F)[1] - info.mindist) < 1.e-7 &&  (!info.nonstandard && isstable(X))

n1 = 3; n2 = 2; m = 2; p = 3; pf = 1;
G = rss(n1,p,m,stable=true); X0 = rss(n2,pf,p,stable=true);  F = X0*G;
@time X, info = glasol(G, F, atol = 1.e-7, poles = [-2, -3], sdeg = -1, mindeg = true); info
@test abs(glinfnorm(X*G-F)[1] - info.mindist) < 1.e-7 &&  (!info.nonstandard && isstable(X)) &&
      order(X) == order(X0)

# n = 3;
# G = rss(n,2,3,stable=true); F = G+0.001*rss(1,2,3,stable=true);
# @time X, info = glasol(G, F, atol = 1.e-7); info
# abs(glinfnorm(X*G-F)[1] - info.mindist)
# @test abs(glinfnorm(X*G-F)[1] -info.mindist) < 0.1 &&  (!info.nonstandard && isstable(X))
# gpole(X)

# F = rss(3,2,3,stable=true); G = F+0.001*rss(3,2,3,stable=true);
# @time X, info = glasol(G, F, atol = 1.e-13); info
# @test abs(glinfnorm(X*G-F)[1] -info.mindist) < 0.1 &&  (!info.nonstandard && isstable(X))
# gpole(X,atol=1.e-10)

## Glover & Packard (SCL,2017)
s = rtf('s');
epsi = -1;
g = (s-1)*(s-epsi)/(s+1)/(s+2); f = 1/(s+1);
sysg = gir(dss(g)); sysf = dss(f);
@time sysx, info = glasol(sysg, sysf, mindeg = true); info
@test glinfnorm(sysx*sysg-sysf)[1] ≈ info.mindist && info.nrank == 1 && info.nonstandard == false &&
      gpole(sysx) ≈ [-1]

s = rtf('s');
epsi = 1;
g = (s-1)*(s-epsi)/(s+1)/(s+2); f = 1/(s+1);
sysg = gir(dss(g)); sysf = dss(f);
@time sysx, info = glasol(sysg, sysf, mindeg = true); info
# check solution
γ0 = 0.25*(1+sqrt(1+8/(1+epsi)))
α = 1 + 2(1 + epsi)*γ0
@test glinfnorm(sysx*sysg-sysf)[1] ≈ info.mindist && info.mindist ≈ γ0 &&
       gpole(sysx) ≈ [-α]

s = rtf('s');
epsi = 1;
gf = [(s-1)*(s-epsi)/(s+1)/(s+2); 1/(s+1)];
sysgf = gir(dss(gf));
@time sysx, info = glasol(sysgf, 1, mindeg = true); info
# check solution
γ0 = 0.25*(1+sqrt(1+8/(1+epsi)))
α = 1 + 2(1 + epsi)*γ0
@test glinfnorm(sysx*sysgf[1,:]-sysgf[2,:])[1] ≈ info.mindist && info.mindist ≈ γ0 &&
       gpole(sysx) ≈ [-α]

# improper L∞ solution for a nonstandard problem
s = rtf('s');
sysg = dss([1/(s+1) 1/(s+2)]); sysf = dss([ 1/(s+1) 0]);
@time sysx, info = glasol(sysg, sysf, atol=1.e-7); info
@test abs(glinfnorm(sysx*sysg-sysf)[1] - info.mindist) < 0.0001 && info.nonstandard ≈ true &&
       gpole(sysx) ≈ [-1.5811388300841827, Inf]

# proper L2 solution for a nonstandard problem
@time sysx, info = glasol(sysg, sysf, atol=1.e-7, L2sol = true); info
@test gl2norm(sysx*sysg-sysf) ≈ info.mindist && info.nonstandard ≈ true &&
       gpole(sysx) ≈ [-1.5811388300841827]

G = rss(0,2,3);
@time X, info = glasol(G, G, atol = 1.e-7); info
@test gl2norm(X*G-G,atolinf=1.e-7) < 1.e-7 && X.D ≈ eye(2)

G = rss(0,2,3);
@time X, info = glasol(G, G, atol = 1.e-7, L2sol = true); info
@test gl2norm(X*G-G,atolinf=1.e-7) < 1.e-7 && X.D ≈ eye(2)

G = rss(3,2,3,stable=true);
@time X, info = glasol(G, G, atol = 1.e-7); info
@test glinfnorm(X*G-G)[1] < 1.e-7 && X.D ≈ eye(2)

G = rss(3,2,3,stable=true);
@time X, info = glasol(G, G, atol = 1.e-7, L2sol = true); info
@test gl2norm(X*G-G,atolinf=1.e-7) < 1.e-7 && X.D ≈ eye(2)

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

@time sysx, info = glasol(sysg, sysf); info
@test abs(glinfnorm(sysx*sysg-sysf,atol=1.e-7)[1] - info.mindist) < 1.e-7 && isstable(sysx)

@time sysx, info = glasol(sysg, sysf, mindeg = true); info
@test abs(glinfnorm(sysx*sysg-sysf,atol=1.e-7)[1] - info.mindist) < 1.e-7 && isstable(sysx)

evals=[-1, -2, -3];
@time sysx, info = glasol(sysg, sysf, poles = evals); info
@test abs(glinfnorm(sysx*sysg-sysf)[1] ≈ info.mindist) < 1.e-7 && sort(gpole(sysx)) ≈ sort(evals)

sdeg = -1;
@time sysx, info = glasol(sysg, sysf, sdeg = sdeg); info
@test abs(glinfnorm(sysx*sysg-sysf)[1] ≈ info.mindist) < 1.e-7 &&
      all(real(gpole(sysx)) .< sdeg*(1-eps()^(1/info.nl)))

end #glasol

@testset "glinfldp" begin

## Parrott's problem
sys1 = rss(0,3,4); sys2 = rss(0,3,2);
@time sysx, mindist = glinfldp(sys1,sys2)
@test glinfnorm([sys1-sysx sys2])[1] ≈ mindist

# Example from [1]
# [1] C.-C. Chu, J. C. Doyle, and E. B. Lee
#     "The general distance problem in H-infinity optimal control theory",
#     Int. J. Control, vol 44, pp. 565-596, 1986.

s = rtf('s');
a = 1; f1 = 1/(s-1); f2 = 1/(s-a); Q0 = -1;
sys1 = dss(f1); sys2 = dss(f2);
@time sysx, mindist = glinfldp([sys1 sys2], 1, atol = 1.e-7, reltol=1.e-9);
@test glinfnorm([sys1-Q0 sys2])[1] ≈  glinfnorm([sys1-sysx sys2])[1] ≈ mindist ≈ 1/a

γ = 1.1;
@time sysx, mindist = glinfldp([sys1 sys2], 1, γ, reltol=1.e-9);
@test glinfnorm([sys1-sysx sys2])[1] ≈ mindist < γ

s = rtf('s');
a = 0.5; f1 = 1/(s-1); f2 = 1/(s-a); Q0 = -((1+a)*s+2*a)/2/(s+a);
sys1 = dss(f1); sys2 = dss(f2); sys0 = dss(Q0);
@time sysx, mindist = glinfldp([sys1 sys2],1, reltol=1.e-9);
@test glinfnorm([sys1-sys0 sys2])[1] ≈  glinfnorm([sys1-sysx sys2])[1] ≈ mindist ≈ 1/a

γ = 1.1/a;
@time sysx, mindist = glinfldp([sys1 sys2], 1, γ, reltol=1.e-9);
@test glinfnorm([sys1-sysx sys2])[1] ≈ mindist < γ

s = rtf('s');
a = 2; f1 = 1/(s-1); f2 = 1/(s-a);
g0 = 1/(2*(a-1))*(-1 + sqrt(a^2+ 4*(a-1)/(a+1))); Q0 = -(g0*s+g0+(a-1)/2)/(s+a);
sys1 = dss(f1); sys2 = dss(f2); sys0 = dss(Q0);
@time sysx, mindist = glinfldp([sys1 sys2],1, reltol=1.e-9);
@test glinfnorm([sys1-sys0 sys2])[1] ≈  glinfnorm([sys1-sysx sys2])[1] ≈ mindist ≈ g0

γ = 1.1*g0;
@time sysx, mindist = glinfldp([sys1 sys2], 1, γ, reltol=1.e-9);
@test glinfnorm([sys1-sysx sys2])[1] ≈ mindist < γ

s = rtf('s'); t = rtfbilin("tustin",Ts=0.5)[1];
a = 0.5; f1 = confmap(1/(s-1),t); f2 = confmap(1/(s-a),t);
Q0 = confmap(-((1+a)*s+2*a)/2/(s+a),t);
sys1 = dss(f1); sys2 = dss(f2); sys0 = dss(Q0);
@time sysx, mindist = glinfldp(gir([sys1 sys2],atol=1.e-7),1, offset=1.e-13, atol=1.e-7, reltol=1.e-9); #fails
@test glinfnorm(gir([sys1-sys0 sys2],atol=1.e-7))[1] ≈  glinfnorm(gir([sys1-sysx sys2],atol=1.e-7))[1] ≈ mindist ≈ 1/a

γ = 1.1/a;
@time sysx, mindist = glinfldp(gir([sys1 sys2],atol=1.e-7), 1, γ, reltol=1.e-9);
@test glinfnorm(gir([sys1-sysx sys2],atol=1.e-7))[1] ≈ mindist < γ

s = rtf('s'); t = rtfbilin("tustin",Ts=0.5)[1];
a = 2; f1 = confmap(1/(s-1),t); f2 = confmap(1/(s-a),t);
g0 = 1/(2*(a-1))*(-1 + sqrt(a^2+ 4*(a-1)/(a+1))); Q0 = -(g0*s+g0+(a-1)/2)/(s+a);
Q0 = confmap(-(g0*s+g0+(a-1)/2)/(s+a),t);
sys1 = dss(f1); sys2 = dss(f2); sys0 = dss(Q0);
@time sysx, mindist = glinfldp(gir([sys1 sys2],atol=1.e-7),1, reltol=1.e-9);
@test glinfnorm(gir([sys1-sys0 sys2],atol=1.e-7))[1] ≈  glinfnorm(gir([sys1-sysx sys2],atol=1.e-7))[1] ≈ mindist ≈ g0

z = rtf('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 = gir(dss(G),atol=1.e-7);
@time sysx, s2 = glinfldp(sys; reltol=1.e-7);
@test glinfnorm(sysx-sys)[1] ≈ s2 && s2 ≈ 8.662176191833835

# poles on the imaginary axis
s = rtf('s');
sys = dss(1/s);
@test_throws ErrorException sysn, s1 = glinfldp(sys)

# Improper continuous-time system sys
s = rtf('s');
sys = dss(s);
@test_throws ErrorException sysn, s1 = glinfldp(sys)

end # glinfldp
@testset "gnehari" begin

sys = rss(0,3,2)
@time sysn, s1 = gnehari(sys,eps())
@test glinfnorm(sysn-sys)[1] == 0

sys = rss(0,10,2)
@time sysn, s1 = gnehari(sys)
@test glinfnorm(sysn-sys)[1] == 0

# poles on the imaginary axis
s = rtf('s');
sys = dss(1/s);
@test_throws ErrorException sysn, s1 = gnehari(sys)

# Improper continuous-time system sys
s = rtf('s');
sys = dss(s);
@test_throws ErrorException sysn, s1 = gnehari(sys)

# poles on the unit circle
z = rtf('z');
sys = dss(1/(z-1));
@test_throws ErrorException sysn, s1 = gnehari(sys)

# stable continuous-time system
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);
@time sysn, s1 = gnehari(sys)
@test glinfnorm(sysn-sys)[1] < 1.e-7 && s1 == 0

z = rtf('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 = gir(dss(G),atol=1.e-7);
@time sysn, s1 = gnehari(sys)
@test glinfnorm(sysn-sys)[1] ≈ s1 && s1 ≈ 8.662176191833835

sys = rdss(50,10,14)';
@time sysn, s1 = gnehari(sys); s1
@test glinfnorm(gminreal(sysn-sys,atol=1.e-10))[1] ≈ s1

@time sysn, s1opt = gnehari(sys,s1*1.01); s1
@test glinfnorm(gminreal(sysn-sys,atol=1.e-10))[1] > s1opt

##  Francis (1987) - optimal Nehari, p. 60
s  = rtf('s');
sys = gir(dss([1/(s^2-1) 4; 1/(s^2-s+1) (s+1)/(s-1)]),atol=1.e-7);
@time sysn, s1 = gnehari(sys);
@test glinfnorm(gminreal(sysn-sys,atol=1.e-7))[1] ≈ s1

##  Francis (1987) - sub-optimal Nehari, p. 128

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

@time sysn, s1 = gnehari(sys);
@test glinfnorm(gminreal(sysn-sys,atol=1.e-7))[1] ≈ s1

end # gnehari
end # model_matching
end # module


text:
import Random
import Statistics: mean
import SelectivelyAdaptiveLasso as SAL
using Test

@testset "FeatureVector" begin for i in 1:100
    n = 100
    U = rand(n)
    X = SAL.FeatureVector(U)
    I = Random.shuffle(1:n)[1:10]

    @test all(U[I] == X.sorted[X.raw_to_sorted[I]])
    @test all(U[X.sorted_to_raw[I]] == X.sorted[I])
    @test X.sorted[sort(X.raw_to_sorted[I])] == U[I[SAL.sortperm_subset(X,I)]]
    @test issorted(X.sorted[X.raw_to_sorted[I[SAL.sortperm_subset(X,I)]]], rev=true)
end end

@testset "1-dim basis selection (right)" begin for i in 1:100
    n,c = 100, 0.4
    I = collect(1:2:n) # must be increasing to keep X sorted
    X = sort(rand(n), rev=true)
    Y = Vector{Float64}(X.≥c)

    Ŷ(X,Y) = X*mean(Y[X])
    sses_brute(X,Y) = [sum((Y.-Ŷ(X.≥xi,Y)).^2) for xi in X]
    sses_clever = SAL.sse_indicator_ols(Y)

    @test all(sses_brute(X,Y) .≈ sses_clever)
    @assert issorted(X[I], rev=true)

    Ĩ, Ỹ, _, x̃, l̃ = SAL.best_split_sorted(I,X[I],Y[I])
    @test Set(Ĩ) == Set((1:n)[X.≥c]) ∩ Set(I)
    @test all(Ỹ .== 1)
    @test x̃ == min(X[I][X[I].≥c]...)
    @test l̃ ≈ min(sses_brute(X[I], Y[I])...)
end end

@testset "changepoint inclusion" begin for i in 1:100
    X = [1.0, 0.9, 0.9, 0.5, 0.1]
    I = 1:length(X)
    Y = rand(length(X))
    Ĩ, Ỹ, _, x̃, l̃ = SAL.best_split_sorted(I,X[I],Y[I])
    @test ((2 ∉ Ĩ) & (3 ∉ Ĩ)) | ((2 ∈ Ĩ) & (3 ∈ Ĩ))
end end

@testset "multidim basis selection (right)" begin for i in 1:100
    c = 0.4
    n,p = 100,3
    j° = p
    U = rand(n,p)
    I = collect(1:2:n) # must be increasing to keep X sorted
    X = [SAL.FeatureVector(x) for x in eachcol(U)]
    Y = Vector{Float64}(U[:,j°] .≥ c)

    Ĩ, Ỹ, _, (j̃, x̃, d), l̃ = SAL.interact_basis(I, X, Y[I])
    @test Set(Ĩ) == Set((1:n)[U[:,j°].≥c]) ∩ Set(I)
    @test all(Ỹ .== 1)
    @test x̃ == min(U[I,j°][U[I,j°].≥c]...)
    @test d == :right
    @test j̃ == j°
end end

@testset "multidim basis selection (left)" begin for i in 1:100
    c = 0.4
    n,p = 100,3
    j° = p
    U = rand(n,p)
    I = collect(1:2:n) # must be increasing to keep X sorted
    X = [SAL.FeatureVector(x) for x in eachcol(U)]
    Y = Vector{Float64}(U[:,j°] .≤ c)

    Ĩ, Ỹ, _, (j̃, x̃, d), l̃ = SAL.interact_basis(I, X, Y[I])
    @test Set(Ĩ) == Set((1:n)[U[:,j°].≤c]) ∩ Set(I)
    @test all(Ỹ .== 1)
    @test x̃ == max(U[I,j°][U[I,j°].≤c]...)
    @test d == :left
    @test j̃ == j°
end end

@testset "build basis" begin for i in 1:100
    c = 0.4
    n,p = 100,3
    U = rand(n,p)
    X = [SAL.FeatureVector(x) for x in eachcol(U)]

    index, basis = SAL.basis_search_random(X)
    n_val = 20
    U_val = rand(n_val,p)
    X_val = [SAL.FeatureVector(x) for x in eachcol(U_val)]
    basis_val = SAL.build_basis(X_val, index)

    build_basis_brute(U, index) = reduce(.*, [d==:right ? U[:,j] .≥ c : U[:,j] .≤ c for (j,c,d) in index])
    bool_basis = build_basis_brute(U, index)
    bool_basis_val = build_basis_brute(U_val, index)
    @test Set(Vector{Int}(basis.nzind)) == Set((1:n)[bool_basis])
    @test Set(Vector{Int}(basis_val.nzind)) == Set((1:n_val)[bool_basis_val])
end end

ext: jl


text:
using NetParam
using Test

@time @testset "NetParam Reading" begin
    touchstone_file_names = readdir( "data/touchstone" )

    # Threads.@threads
    for filename in touchstone_file_names
        filepath = string( "data/touchstone/", filename )

        # for excluding mixer files
        if !occursin( "px", lowercase( filename ) )
            # print( filename, '\n' )
            NetParam.read_touchstone( filepath )
        end
    end
end

@time @testset "mixertest" begin
    @test true
end


text:
using Treesitter

parser = Parser()
l = set_language(parser, tree_sitter_c())
tree = Treesitter.parse_filewith(
    parser,
    joinpath(@__DIR__, "corpus/easystruct.h"),
)
rnode = tree_getroot(tree)

nodetype = node_type(rnode)
nodesym = node_symbol(rnode)
nbstart = node_startbyte(rnode)
nbstop = node_endbyte(rnode)

# Length of node text is equal to the number of bytes in the raw sourcecode.
slen = parser.buffer |> length
nlen = nbstop - nbstart |> Int
@assert slen == nlen

root_children = node_children(rnode)
@assert length(root_children) == 4

rchildrentype = node_type.(root_children)
rsymbol = node_symbol.(root_children)
rbstart = node_startbyte.(root_children)
rbstop = node_endbyte.(root_children)
rslen = rbstop .- rbstart

rnodestrings = node_string.(root_children)
rcisnull = node_isnull.(root_children)
rcisnamed = node_isnamed.(root_children)
rcisextra = node_isextra.(root_children)


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

module GeoStatsDevTools

using Reexport

@warn "GeoStatsDevTools.jl is deprecated. Please use GeoStatsBase.jl"

@reexport using GeoStatsBase

end


text:
using SMC, Base.Test

srand(123)
N = 50000
x = randn(N)
w = rand(N); w /= sum(w)

ps1 = Particles(x, w)
ps2 = multinomialresampling(ps1)

mps1 = mean(ps1)
mps2 = mean(ps2)

# mean and variance of resampled set should be similar to initial.
@test isapprox(mean(ps1), mean(ps2), atol=0.1)
@test isapprox(var(ps2.x), sum(ps1.x.^2 .* ps1.w)-mps1^2, atol=0.1)


text:
using Test
using NCDatasets
using DataStructures

function example_file(i,array)
    fname = tempname()
    @debug "fname $fname"

    Dataset(fname,"c") do ds
        # Dimensions

        ds.dim["lon"] = size(array,1)
        ds.dim["lat"] = size(array,2)
        ds.dim["time"] = Inf

        # Declare variables

        ncvar = defVar(ds,"var", Float64, ("lon", "lat", "time"))
        ncvar.attrib["field"] = "u-wind, scalar, series"
        ncvar.attrib["units"] = "meter second-1"
        ncvar.attrib["long_name"] = "surface u-wind component"
        ncvar.attrib["time"] = "time"
        ncvar.attrib["coordinates"] = "lon lat"

        nclat = defVar(ds,"lat", Float64, ("lat",))
        nclat.attrib["units"] = "degrees_north"

        nclon = defVar(ds,"lon", Float64, ("lon",))
        nclon.attrib["units"] = "degrees_east"
        nclon.attrib["modulo"] = 360.0

        nctime = defVar(ds,"time", Float64, ("time",))
        nctime.attrib["long_name"] = "surface wind time"
        nctime.attrib["field"] = "time, scalar, series"
        nctime.attrib["units"] = "days since 2000-01-01 00:00:00 GMT"

        # Global attributes

        ds.attrib["history"] = "foo"

        # Define variables

        g = defGroup(ds,"mygroup")
        ncvarg = defVar(g,"var", Float64, ("lon", "lat", "time"))
        ncvarg.attrib["field"] = "u-wind, scalar, series"
        ncvarg.attrib["units"] = "meter second-1"
        ncvarg.attrib["long_name"] = "surface u-wind component"
        ncvarg.attrib["time"] = "time"
        ncvarg.attrib["coordinates"] = "lon lat"

        ncvar[:,:,1] = array
        ncvarg[:,:,1] = array.+1
        #nclon[:] = 1:size(array,1)
        #nclat[:] = 1:size(array,2)
        nctime[:] = i
    end
    return fname
end

A = randn(2,3,1)

fname = example_file(1,A)

ds = Dataset(fname)

info = NCDatasets.metadata(ds)

#@show info

dds = DeferDataset(fname);
varname = "var"
datavar = variable(dds,varname);
data = datavar[:,:,:]

@test A == datavar[:,:,:]

@test dds.attrib["history"] == "foo"
@test datavar.attrib["units"] == "meter second-1"

@test dimnames(datavar) == ("lon", "lat", "time")
lon = variable(dds,"lon");
@test lon.attrib["units"] == "degrees_east"
@test size(lon) == (size(data,1),)

datavar = dds[varname]
@test A == datavar[:,:,:]

@test dimnames(datavar) == ("lon", "lat", "time")

@test dds.dim["lon"] == size(A,1)
@test dds.dim["lat"] == size(A,2)
@test dds.dim["time"] == size(A,3)

@test dds.dim["time"] == size(A,3)
close(dds)

# show
dds_buf = IOBuffer()
show(dds_buf,dds)

ds_buf = IOBuffer()
show(ds_buf,dds)

@test String(take!(dds_buf)) == String(take!(ds_buf))

#=
# write
dds = Dataset(fnames,"a");
dds["datavar"][2,2,:] = 1:length(fnames)

for n = 1:length(fnames)
    Dataset(fnames[n]) do ds
        @test ds["datavar"][2,2,1] == n
    end
end

dds.attrib["history"] = "foo2"
sync(dds)

Dataset(fnames[1]) do ds
    @test ds.attrib["history"] == "foo2"
end

@test_throws NCDatasets.NetCDFError Dataset(fnames,"not-a-mode")

@test keys(dds) == ["datavar", "lat", "lon", "time"]
@test keys(dds.dim) == ["lon", "lat", "time"]
@test NCDatasets.groupname(dds) == "/"
@test size(dds["var"]) == (2, 3, 3)
@test size(dds["var"].var) == (2, 3, 3)
@test name(dds["var"].var) == "var"
@test NCDatasets.groupname(dds.group["mygroup"]) == "mygroup"

# create new dimension in all files
dds.dim["newdim"] = 123;
sync(dds);
Dataset(fnames[1]) do ds
    @test ds.dim["newdim"] == 123
end
close(dds)
=#
nothing


text:
module Compat

using Base.Meta

if VERSION >= v"0.4.0-dev+3184"
    include("ngenerate.jl")
    using .CompatCartesian
    export @ngenerate, @nsplat
end

if VERSION < v"0.4.0-dev+1624"
    eachindex(A::AbstractArray) = 1:length(A)
    export eachindex
end

if VERSION < v"0.4.0-dev+1419"
    export UInt, UInt8, UInt16, UInt32, UInt64, UInt128
    const UInt = Uint
    const UInt8 = Uint8
    const UInt16 = Uint16
    const UInt32 = Uint32
    const UInt64 = Uint64
    const UInt128 = Uint128
end

if VERSION < v"0.4.0-dev+1387"
    typealias AbstractString String
    export AbstractString
end

if VERSION < v"0.4.0-dev+3324"
    typealias AssertionError ErrorException
    export AssertionError
end

if VERSION < v"0.4.0-dev+412"
    eval(Base, :(const IPAddr = IpAddr))
end

if VERSION < v"0.4.0-dev+2197"
    Base.IPv4(ipstr::AbstractString) = Base.parseipv4(ipstr)
    Base.IPv6(ipstr::AbstractString) = Base.parseipv6(ipstr)
end

if VERSION < v"0.4.0-dev+2200"
    function Base.isless{T<:Base.IPAddr}(a::T, b::T)
        return isless(a.host, b.host)
    end
end

if VERSION < v"0.4.0-dev+980"
    macro Dict(pairs...)
        esc(Expr(:dict, pairs...))
    end
    macro AnyDict(pairs...)
        esc(Expr(:typed_dict, :(Any=>Any), pairs...))
    end

    Base.Dict(kv) = dict_with_eltype(kv, eltype(kv))
    dict_with_eltype{K,V}(kv, ::Type{(K,V)}) = Dict{K,V}(kv)
    dict_with_eltype(kv, t) = Dict{Any,Any}(kv)
else
    macro Dict(pairs...)
        esc(Expr(:call, :Dict, pairs...))
    end
    macro AnyDict(pairs...)
        esc(Expr(:call, :(Base.AnyDict), pairs...))
    end
end

import Base: round, ceil, floor, trunc
if VERSION < v"0.4.0-dev+1827"
    for (fnew,fold) in ((:round, :iround), (:ceil, :iceil), (:floor, :ifloor), (:trunc, :itrunc))
        @eval begin
            ($fnew){T<:Integer}(::Type{T}, x::Integer) = convert(T, x)  # ambiguity resolution with digits/base version, not all old methods defined
            ($fnew){T<:Integer}(::Type{T}, x::Real) = ($fold)(T, x)
            ($fnew){T<:Integer}(::Type{T}, x::Rational) = convert(T, ($fold)(x)) # no e.g. iround(::Type{T}, x::Rational) is defined in 0.3
            function ($fnew){T,R}(::Type{T}, x::AbstractArray{R,1})
                [ ($fnew)(T, x[i])::T for i = 1:length(x) ]
            end
            function ($fnew){T,R}(::Type{T}, x::AbstractArray{R,2})
                [ ($fnew)(T, x[i,j])::T for i = 1:size(x,1), j = 1:size(x,2) ]
            end
            function ($fnew){T}(::Type{T}, x::AbstractArray)
                reshape([ ($fnew)(T, x[i])::T for i = 1:length(x) ], size(x))
            end
        end
    end
end

if VERSION < v"0.4.0-dev+1884"
    randexp(rng::MersenneTwister) = Base.Random.randmtzig_exprnd(rng)
    randexp() = Base.Random.randmtzig_exprnd()
    export randexp
end

if VERSION < v"0.4.0-dev+2014"
    sizehint! = Base.sizehint
    export sizehint!
end

if VERSION < v"0.4.0-dev+3413"
    # based on pipe in base/process.jl:
    import Base: AbstractCmd, Redirectable
    function pipe(cmd::AbstractCmd; stdin=nothing, stdout=nothing, stderr=nothing, append::Bool=false)
        if append && stdout === nothing && stderr === nothing
            error("append set to true, but no output redirections specified")
        end
        if stdin !== nothing
            cmd = stdin |> cmd
        end
        if stdout !== nothing
            cmd = append ? cmd >> stdout : cmd |> stdout
        end
        if stderr !== nothing
            cmd = append ? cmd .>> stderr : cmd .> stderr
        end
        return cmd
    end
    pipe(cmd::AbstractCmd, dest) = pipe(cmd, stdout=dest)
    pipe(src::Union(Redirectable,AbstractString), cmd::AbstractCmd) = pipe(cmd, stdin=src)
    pipe(a, b, c, d...) = pipe(pipe(a,b), c, d...)
    export pipe
end

function rewrite_dict(ex)
    length(ex.args) == 1 && return ex

    f = ex.args[1]
    if isexpr(f, :curly)
        newex = Expr(:typed_dict, :($(f.args[2])=>$(f.args[3])))
    else
        newex = Expr(:dict)
    end

    for i = 2:length(ex.args)
        pair = ex.args[i]
        !isexpr(pair, :(=>)) && return ex
        push!(newex.args, pair)
    end
    newex
end

# rewrite Julia 0.4-style split or rsplit (str, splitter; kws...)
# into 0.2/0.3-style positional arguments
function rewrite_split(ex, f)
    limit = nothing
    keep = nothing
    for i in 4:length(ex.args)
        if isexpr(ex.args[i], :kw)
            kw = ex.args[i].args
            if kw[1] == :limit
                limit = kw[2]
            elseif kw[1] == :keep
                keep = kw[2]
            end
        end
    end
    if limit == nothing
        if keep == nothing
            return Expr(:call, f, ex.args[2], ex.args[3])
        else
            return Expr(:call, f, ex.args[2], ex.args[3], keep)
        end
    else
        if keep == nothing
            return Expr(:call, f, ex.args[2], ex.args[3], limit)
        else
            return Expr(:call, f, ex.args[2], ex.args[3], limit, keep)
        end
    end
end

if VERSION < v"0.4.0-dev+707"
    macro inline(ex)
        esc(ex)
    end
end

if VERSION < v"0.4.0-dev+2056"
    macro noinline(ex)
        esc(ex)
    end
end

if VERSION < v"0.4.0-dev+2440"
    bitrand(r::AbstractRNG, dims::Dims)   = rand!(r, BitArray(dims))
    bitrand(r::AbstractRNG, dims::Int...) = rand!(r, BitArray(dims))
    bitrand(dims::Dims)   = rand!(BitArray(dims))
    bitrand(dims::Int...) = rand!(BitArray(dims))
    export bitrand
    Base.rand(::Type{Bool}) = randbool()
end

if VERSION < v"0.4.0-dev+2485"
    startswith = Base.beginswith
    export startswith
end

if VERSION < v"0.4.0-dev+3609"
    fieldnames(x::DataType) = collect(x.names)
    function fieldnames(v)
        t = typeof(v)
        if !isa(t, DataType)
            throw(ArgumentError("cannot call fieldnames() on a non-composite type"))
        end
        return fieldnames(t)
    end
    export fieldnames
end

if VERSION < v"0.4.0-dev+3874"
    Base.parse{T<:Integer}(::Type{T}, c::Char) = parseint(T, c)
    Base.parse{T<:Integer}(::Type{T}, c::Char, base::Integer) = parseint(T, c, base)
    Base.parse{T<:Integer}(::Type{T}, s::AbstractString) = parseint(T, s)
    Base.parse{T<:Integer}(::Type{T}, s::AbstractString, base::Integer) = parseint(T, s, base)
    Base.parse{T<:Union(Float32,Float64)}(::Type{T}, s::AbstractString) = parsefloat(T, s)
end

if VERSION < v"0.4.0-dev+4539"
    Base.parse(::Type{BigFloat}, s::AbstractString) = BigFloat(s)
end

if VERSION < v"0.4.0-dev+3710"
    const unsafe_convert = Base.convert
else
    import Base.unsafe_convert
end

if VERSION < v"0.4.0-dev+3732"
    const calltypes = Dict{Symbol,Symbol}()
    for (k,v) in [(:Integer,:integer),
                   (:Signed,:signed),
                   (:Unsigned,:unsigned),
                   (:Int,:int),
                   (:Int8,:int8),
                   (:Int16,:int16),
                   (:Int32,:int32),
                   (:Int64,:int64),
                   (:Int128,:int128),
                   (:UInt,:uint),
                   (:UInt8,:uint8),
                   (:UInt16,:uint16),
                   (:UInt32,:uint32),
                   (:UInt64,:uint64),
                   (:UInt128,:uint128),
                   (:Float16,:float16),
                   (:Float32,:float32),
                   (:Float64,:float64),
                   (:Complex32,:complex32),
                   (:Complex64,:complex64),
                   (:Complex128,:complex128),
                   (:Char,:char),
                   (:Bool,:bool)]
        calltypes[k] = v
    end
end

function _compat(ex::Expr)
    if ex.head == :call
        f = ex.args[1]
        if VERSION < v"0.4.0-dev+980" && (f == :Dict || (isexpr(f, :curly) && length(f.args) == 3 && f.args[1] == :Dict))
            ex = rewrite_dict(ex)
        elseif VERSION < v"0.4.0-dev+129" && (f == :split || f == :rsplit) && length(ex.args) >= 4 && isexpr(ex.args[4], :kw)
            ex = rewrite_split(ex, f)
        elseif VERSION < v"0.4.0-dev+3732" && haskey(calltypes, f) && length(ex.args) > 1
            T = ex.args[1]
            if T in (:Complex32, :Complex64, :Complex128)
                ex = Expr(:call, calltypes[T], ex.args[2:end]...)
            else
                ex = Expr(:(::), Expr(:call, :convert, T, ex.args[2:end]...), T)
            end
        elseif VERSION < v"0.4.0-dev+3732" && (f == :map && haskey(calltypes, ex.args[2]) && length(ex.args) > 2)
            ex = Expr(:call, calltypes[ex.args[2]], ex.args[3:end]...)
        elseif VERSION < v"0.4.0-dev+1419" && isexpr(f, :curly) && f.args[1] == :Ptr && length(ex.args) == 2 && ex.args[2] == 0
            ex = Expr(:call, :zero, f)
        elseif VERSION < v"0.4.0-dev+4356" && f == :chol
            s = ex.args[3]
            if isexpr(s, :curly) && s.args[1] == :Val
                ex = Expr(:call, :chol, ex.args[2], s.args[2])
            end
        elseif VERSION < v"0.4.0-dev+4739" && isexpr(f, :curly) && (f.args[1] == :Vector || f.args[1] == :Array)
            if f.args[1] == :Vector && length(ex.args) == 1
                ex = Expr(:call, :Array, f.args[2], 0)
            else
                ex = Expr(:call, :Array, f.args[2], ex.args[2:end]...)
            end
        end
    elseif ex.head == :curly
        f = ex.args[1]
        if VERSION < v"0.4.0-dev+4319" && f == :Tuple
            args = ex.args[2:end]
            has_ellipsis = any(args) do arg
                isa(arg, Expr) && (arg.head == :...)
            end
            ex = if has_ellipsis
                Expr(:call, TopNode(:tuple), args...)
            else
                Expr(:tuple, args...)
            end
        end
    elseif ex.head == :macrocall
        f = ex.args[1]
        if f == symbol("@generated") && VERSION < v"0.4.0-dev+4387"
            f = ex.args[2]
            if isexpr(f, :function)
                ex = Expr(:stagedfunction, f.args...)
            end
        end
    elseif VERSION < v"0.4.0-dev+5322" && ex.head == :(::) && isa(ex.args[end], Symbol)
        # Replace Base.Timer with Compat.Timer2 in type declarations
        if ex.args[end] == :Timer || ex.args[end] == :(Base.Timer)
            ex.args[end] = :(Compat.Timer2)
        end
    end
    return Expr(ex.head, map(_compat, ex.args)...)
end
_compat(ex) = ex

macro compat(ex)
    esc(_compat(ex))
end

export @compat, @inline, @noinline

if VERSION < v"0.4.0-dev+656"
    include("nullable.jl")
end

if VERSION < v"0.4.0-dev+3864"
    function tryparse(T::Type{Float32}, s)
        r = Array(T,1)
        float32_isvalid(s, r) ? Nullable(r[1]) : Nullable{Float32}()
    end

    function tryparse(T::Type{Float64}, s)
        r = Array(T,1)
        float64_isvalid(s, r) ? Nullable(r[1]) : Nullable{Float64}()
    end

    function tryparse{T<:Integer}(::Type{T}, s)
        local ret
        try
            val = parse(T, s)
            ret = Nullable{T}(val)
        catch
            ret = Nullable{T}()
        end
        return ret
    end
    export tryparse
end

if VERSION < v"0.4.0-dev+3844"
    @eval module Libc
        const FILE = Base.CFILE
        const calloc = Base.c_calloc
        const free = Base.c_free
        const malloc = Base.c_malloc
        const realloc = Base.c_realloc
        using Base: munmap, systemsleep
        export FILE, TmStruct, calloc, errno, flush_cstdio, free,
               gethostname, getpid, malloc, msync, realloc, strerror,
               strftime, strptime, systemsleep, time
    end
    @eval module Libdl
        using Base.Sys: dlext, dllist, dlpath
        export DL_LOAD_PATH, RTLD_DEEPBIND, RTLD_FIRST, RTLD_GLOBAL, RTLD_LAZY,
               RTLD_LOCAL, RTLD_NODELETE, RTLD_NOLOAD, RTLD_NOW, dlclose, dlext,
               dllist, dlopen, dlopen_e, dlpath, dlsym, dlsym_e, find_library
    end
    export Libc, Libdl
else
    import Base: Libc, Libdl
end

if VERSION < v"0.4.0-dev+2418"
    function findprev(A, start)
        for i = start:-1:1
            A[i] != 0 && return i
        end
        0
    end
    findlast(A) = findprev(A, length(A))
    function findprev(A, v, start)
        for i = start:-1:1
            A[i] == v && return i
        end
        0
    end
    findlast(A, v) = findprev(A, v, length(A))
    function findprev(testf::Function, A, start)
        for i = start:-1:1
            testf(A[i]) && return i
        end
        0
    end
    findlast(testf::Function, A) = findprev(testf, A, length(A))
    export findprev, findlast
end

if VERSION < v"0.4.0-dev+4524"
    isdiag(A::AbstractMatrix) = istril(A) && istriu(A)
    isdiag(x::Number) = true
    export isdiag
end

if VERSION < v"0.4.0-dev+4502"
    keytype{K,V}(::Associative{K,V}) = K
    valtype{K,V}(::Associative{K,V}) = V
    export keytype, valtype
end

if VERSION < v"0.4.0-dev+2254"
    immutable Val{T} end
    export Val
end

if VERSION < v"0.4.0-dev+4603"
    # used for C string arguments to ccall
    # (in Julia 0.4, these types also check for embedded NUL chars)
    const Cstring = Ptr{Cchar}
    const Cwstring = Ptr{Cwchar_t}
    export Cstring, Cwstring
end

if VERSION < v"0.4.0-dev+2823"
    fma(x::Number,y::Number,z::Number) = x*y + z
    export fma
end

if VERSION < v"0.4.0-dev+2861"
    muladd(x::Number,y::Number,z::Number) = x*y + z
    export muladd
end

if VERSION < v"0.4.0-dev+4734"
    function is_valid_utf32(str::Union(Vector{Char}, Vector{UInt32}))
        for i=1:length(str)
            @inbounds if !is_valid_char(reinterpret(UInt32, str[i])) ; return false ; end
        end
        return true
    end
    is_valid_utf32(str::UTF32String) = is_valid_utf32(str.data)
    export is_valid_utf32
end

if VERSION < v"0.4.0-dev+4939"
    import Base.isvalid
    isvalid(ch::Char)         = is_valid_char(ch)
    isvalid(str::ASCIIString) = is_valid_ascii(str)
    isvalid(str::UTF8String)  = is_valid_utf8(str)
    isvalid(str::UTF16String) = is_valid_utf16(str)
    isvalid(str::UTF32String) = is_valid_utf32(str)
    isvalid(::Type{Char}, ch)         = is_valid_char(ch)
    isvalid(::Type{ASCIIString}, str) = is_valid_ascii(str)
    isvalid(::Type{UTF8String},  str) = is_valid_utf8(str)
    isvalid(::Type{UTF16String}, str) = is_valid_utf16(str)
    isvalid(::Type{UTF32String}, str) = is_valid_utf32(str)
    export isvalid
end

if VERSION < v"0.4.0-dev+5322"
    include("timer.jl")
end

if VERSION < v"0.4.0-dev+5688"
    typealias Irrational MathConst
    @eval const $(symbol("@irrational")) = getfield(Base, symbol("@math_const"))
    export Irrational
else
    import Base.@irrational
end

if VERSION < v"0.4.0-dev+5305"
    function Base.gc_enable(on::Bool)
        state = ccall(:jl_gc_is_enabled, Cint, ()) != 0
        on ? gc_enable() : gc_disable()
        state
    end
end

end # module

max_stars_repo_name: ambujone/DIDebug

max_issues_repo_name: ambujone/DIDebug

max_forks_repo_name: ambujone/DIDebug


text:
# Xpress callbacks

mutable struct CallbackData
    model_root::Model # should not use operations here
    data::Any # data for user
    model::Model# local model # ptr_model::Ptr{Nothing}
end

mutable struct MessageData
    data::Any # data for user
    msg::Cstring
    len::Cint
    msgtype::Cint
end

function setcbintsol_wrapper(ptr_model::Ptr{Nothing}, my_object::Ptr{Nothing})
    callback, model, data = unsafe_pointer_to_objref(my_object)
    callback(CallbackData(model, data, Model(ptr_model;finalize_env = false)))
    return convert(Cint,0)
end

set_callback_intsol!(model::Model, callback::Function) = set_callback_intsol!(model, callback, nothing)
function set_callback_intsol!(model::Model, callback::Function, data::Any)
    xprscallback = cfunction(setcbintsol_wrapper, Cint, (Ptr{Nothing}, Ptr{Nothing}))
    usrdata = (callback, model, data)
    ret = @xprs_ccall(addcbintsol, Cint, (Ptr{Nothing}, Ptr{Nothing}, Any), model.ptr_model, xprscallback, usrdata)
    if ret != 0
        throw(XpressError(model))
    end
    # we need to keep a reference to the callback function
    # so that it isn't garbage collected
    push!(model.callback,usrdata)
    return nothing
end

function addcbpreintsol_wrapper(ptr_model::Ptr{Nothing}, my_object::Ptr{Nothing})
    callback, model, data = unsafe_pointer_to_objref(my_object)#::(Function,Model,Any)
    callback(CallbackData(model, data, Model(ptr_model;finalize_env = false)))
    return convert(Cint,0)
end

set_callback_preintsol!(model::Model, callback::Function) = set_callback_preintsol!(model, callback, nothing)
function set_callback_preintsol!(model::Model, callback::Function, data::Any)

    xprscallback = cfunction(addcbpreintsol_wrapper, Cint, (Ptr{Nothing}, Ptr{Nothing}))
    usrdata = (callback, model, data)
    ret = @xprs_ccall(addcbpreintsol, Cint, (Ptr{Nothing}, Ptr{Nothing}, Any), model.ptr_model, xprscallback, usrdata)
    if ret != 0
        throw(XpressError(model))
    end

    # we need to keep a reference to the callback function
    # so that it isn't garbage collected
    push!(model.callback,usrdata)
    return nothing
end

function addmessage_wrapper(ptr_model::Ptr{Nothing}, my_object::Ptr{Nothing}, msg::Cstring, len::Cint, msgtype::Cint)
    callback, data = unsafe_pointer_to_objref(my_object)#::(Function,Model,Any)
    callback(MessageData(data, msg, len, msgtype))
    return convert(Cint,0)
end

set_callback_message!(model::Model, callback::Function) = set_callback_message!(model, callback, nothing)
function set_callback_message!(model::Model, callback::Function, data::Any)
    xprscallback = cfunction(addmessage_wrapper, Cint, (Ptr{Nothing}, Ptr{Nothing}, Cstring, Cint, Cint))
    usrdata = (callback, data)
    ret = @xprs_ccall(addcbmessage, Cint, (Ptr{Nothing}, Ptr{Nothing}, Any), model.ptr_model, xprscallback, usrdata)
    if ret != 0
        throw(XpressError(model))
    end
    # we need to keep a reference to the callback function
    # so that it isn't garbage collected
    push!(model.callback,usrdata)
    return nothing
end


text:
const drugs = ["Citalopram","Methysergide","Altanserin","Way_100135","SB242084","Optogenetic"]

const Cit_per = Date(2016,03,14):Day(1):Date(2016,03,18)
const Meth_per = Date(2016,03,21):Day(1):Date(2016,03,25)
const Alt_per = Date(2016,03,28):Day(1):Date(2016,04,01)
const WAY_per = Date(2016,04,04):Day(1):Date(2016,04,08)
const SB_per = Date(2016,05,23):Day(1):Date(2016,05,27)
const Opto_per = Date(2016,07,19):Day(1):Date(2016,08,04)

# const Phase_calendar = OrderedDict()
# for (d,p) in zip(drugs,[Cit_per,Meth_per,Alt_per,WAY_per,SB_per, Opto_per])
#     Phase_calendar[d] = p
# end

const Phase_Calendar = OrderedDict()
for (dic,drug) in zip([Cit_per,Meth_per,Alt_per,WAY_per,SB_per, Opto_per],drugs)
    for x in dic
        Phase_Calendar[x] = drug
    end
end
###

const Drug_assignment = OrderedDict()

for d in drugs
    Drug_assignment[d] = OrderedDict()
end
for d in drugs
    for g in ["Group A", "Group B","Group C", "Group D"]
        Drug_assignment[d][g] = OrderedDict()
    end
end

Drug_assignment["Citalopram"]["Group A"] = ["None","PreVehicle","Citalopram","PostVehicle","Saline"]
Drug_assignment["Citalopram"]["Group B"] = ["None","Saline","PreVehicle","Citalopram","PostVehicle"]
Drug_assignment["Citalopram"]["Days"] = Cit_per
Drug_assignment["Methysergide"]["Group A"] = ["None","PreVehicle","Methysergide","PostVehicle","Saline"]
Drug_assignment["Methysergide"]["Group B"] = ["None","Saline","PreVehicle","Methysergide","PostVehicle"]
Drug_assignment["Methysergide"]["Days"] = Meth_per
Drug_assignment["Altanserin"]["Group A"] = ["None","PreVehicle","Altanserin","PostVehicle","Saline"]
Drug_assignment["Altanserin"]["Group B"] = ["None","Saline","PreVehicle","Altanserin","PostVehicle"]
Drug_assignment["Altanserin"]["Days"] = Alt_per
Drug_assignment["Way_100135"]["Group A"] = ["None","PreVehicle","Way_100135","PostVehicle","Saline"]
Drug_assignment["Way_100135"]["Group B"] = ["None","Saline","PreVehicle","Way_100135","PostVehicle"]
Drug_assignment["Way_100135"]["Days"] = WAY_per
Drug_assignment["SB242084"]["Group A"] = ["None","Saline","PreVehicle","SB242084","PostVehicle"]
Drug_assignment["SB242084"]["Group B"] = ["None","PreVehicle","SB242084","PostVehicle","Saline"]
Drug_assignment["SB242084"]["Days"] = SB_per
Drug_assignment["Optogenetic"]["Group C"] = ["Saline", "Saline", "Saline","Saline", "SB242084_opt", "None",
                                          "Saline", "Saline", "Saline", "SB242084_opt", "Saline", "Saline",
                                           "Saline", "SB242084_opt", "Saline", "Saline", "SB242084_opt"]
Drug_assignment["Optogenetic"]["Group D"] = ["Saline", "Saline", "Saline", "SB242084_opt", "Saline", "None",
                                          "Saline", "SB242084_opt","Saline", "Saline", "Saline", "SB242084_opt",
                                          "Saline", "Saline", "Saline","SB242084_opt", "Saline"]
Drug_assignment["Optogenetic"]["Days"] = Opto_per

##
const Procedure =  DataFrame(Phase = String[], Group = String[], Day = Date[],
    Treatment = String[],PhaseDay = Int64[])

for drug in drugs
    for (i,day) in enumerate(Drug_assignment[drug]["Days"])
        for group in ["Group A", "Group B"]
            if isempty(Drug_assignment[drug][group])
                continue
            else
                treatment = Drug_assignment[drug][group][i]
                push!(Procedure,[drug,group,day,treatment,i])
            end
        end
    end
end

##
function get_treatment(drug,group,day)
    phase = get(Drug_assignment,drug,"None")
    if phase == "None"
        return "None"
    else
        v = collect(Drug_assignment[drug]["Days"])
        idx = findall(t-> t == day,v)
        if length(idx)>1
            println("weird amount of dates matching")
        else
            idx = idx[1]
        end
        return collect(Drug_assignment[drug][group])[idx]
    end
end

########################################
# const Sequence = OrderedDict()
# Sequence["Group A"] = ["none","sal","+","o","sal"]
# Sequence["Group B"] = ["none","sal","o","+","sal"]
#
# const Full_Drug_Calendar = OrderedDict()
# for drug in drugs
#     Full_Drug_Calendar[drug] = OrderedDict()
# end
#
# for drug in drugs
#     for g in ["Group A","Group B"]
#         Full_Drug_Calendar[drug][g] = OrderedDict()
#     end
# end
#
# for drug in drugs
#     for g in ["Group A","Group B"]
#         for (day,inj) in zip(Phase_calendar[drug],Sequence[g])
#             Full_Drug_Calendar[drug][g][day] = inj
#         end
#     end
# end

function get_injection(phase,group,day)
    phase = get(Full_Drug_Calendar,phase,"none")
    if phase == "none"
        return "none"
    else
        return get(phase[group],day,"none")
    end
end


text:
mutable struct BufferView
    buffer::Int
    byteOffset::Union{Nothing,Int}
    byteLength::Int
    byteStride::Union{Nothing,Int}
    target::Union{Nothing,Int}
    name::Union{Nothing,String}
    extensions::Union{Nothing,Dict}
    extras::Union{Nothing,Dict}
    function BufferView()
        obj = new()
        obj.byteOffset = nothing
        obj.byteStride = nothing
        obj.target = nothing
        obj.name = nothing
        obj.extensions = nothing
        obj.extras = nothing
        obj
    end
end

function Base.getproperty(obj::BufferView, sym::Symbol)
    x = getfield(obj, sym)
    sym === :byteOffset && x === nothing && return 0
    return x
end

function Base.setproperty!(obj::BufferView, sym::Symbol, x)
    if sym === :buffer
        x ≥ 0 || throw(ArgumentError("the index of the buffer should be ≥ 0"))
    elseif sym === :byteOffset && x !== nothing
        x ≥ 0 || throw(ArgumentError("byteOffset should be ≥ 0"))
    elseif sym === :byteLength
        x ≥ 1 || throw(ArgumentError("byteLength should be ≥ 1"))
    elseif sym === :byteStride && x !== nothing
        4 ≤ x ≤ 252 || throw(ArgumentError("byteStride should be ≥ 4 and ≤ 252"))
    elseif sym === :target && x !== nothing
        x === ARRAY_BUFFER || x === ELEMENT_ARRAY_BUFFER ||
            throw(ArgumentError("target should be ARRAY_BUFFER(34962) or ELEMENT_ARRAY_BUFFER(34963)"))
    end
    setfield!(obj, sym, x)
end

JSON3.StructType(::Type{BufferView}) = JSON3.Mutable()
StructTypes.omitempties(::Type{BufferView}) = (:byteOffset, :byteStride, :target, :name, :extensions, :extras)