text:
module Test_conversions

using DescriptorSystems
using LinearAlgebra
using Polynomials
using Test

println("Test_conversions")

@testset "gbilin" begin

s = rtf('s'); # define the complex variable s
Gc = [s^2 s/(s+1); 0 1/s] # define the 2-by-2 improper Gc(s)
@time sysc = dss(Gc);            # build continuous-time descriptor system realization
g = (s+0.01)/(1+0.01*s);

@time syst, = gbilin(sysc,g,atol = 1.e-7,minimal=true);

@test opnorm(evalfr(sysc-syst,1)) < 1.e-7
# compute the ν-gap
@time nugap = gnugap(sysc,syst,atol=1.e-4)[1];
@test nugap < 0.02

z = rtf('z');
Gd = [z^2 z/(z-2); 0 1/z]     # define the 2-by-2 improper Gd(z)
@time syst, = gbilin(dss(Gd),rtfbilin("cayley")[2],atol=1.e-7,minimal = true);
@time sysi, = gbilin(syst,rtfbilin("c2d")[1],atol = 1.e-7,minimal = true);

@test ghinfnorm(gminreal(dss(Gd)-sysi,atol=1.e-7))[1] < 1.e-7

for Ty in (Float64,Complex{Float64})

type = "c2d";
sys = rss(T = Ty, 3,2,3);
@time g, ginv = rtfbilin(type)
@time syst, gi = gbilin(sys,g,atol=1.e-7);
@time sysi, g1 = gbilin(syst,ginv,atol=1.e-7);
@test iszero(sys-sysi,atol=1.e-7) && gi == ginv && g == g1

type = "c2d";
sys = rdss(T = Ty, 3,2,3);
@time g, ginv = rtfbilin(type)
@time syst, gi = gbilin(sys,g,atol=1.e-7);
@time sysi, g1 = gbilin(syst,ginv,atol=1.e-7);
@test iszero(sys-sysi,atol=1.e-7) && gi == ginv && g == g1

type = "cayley";
sys = rss(T = Ty, 3,2,3,disc = true);
@time g, ginv = rtfbilin(type)
@time syst, g1 = gbilin(sys,ginv,atol=1.e-7);
@time sysi, gi1 = gbilin(syst,g,atol=1.e-7);
@test iszero(sys-sysi,atol=1.e-7) && g1 == g && ginv == gi1

type = "cayley";
sys = rdss(T = Ty, 3,2,3,disc = true);
@time g, ginv = rtfbilin(type)
@time syst, g1 = gbilin(sys,ginv,atol=1.e-7);
@time sysi, gi1 = gbilin(syst,g,atol=1.e-7);
@test iszero(sys-sysi,atol=1.e-7) && g1 == g && ginv == gi1

type = "tustin";
sys = rss(T = Ty, 3,2,3);
@time g, ginv = rtfbilin(type,Ts=0.1)
@time syst, gi = gbilin(sys,g,atol=1.e-7);
@time sysi, g1 = gbilin(syst,ginv,atol=1.e-7);
@test iszero(sys-sysi,atol=1.e-7) && gi == ginv && g == g1

type = "lft";
sys = rss(T = Ty, 3,2,3);
val = rand(Ty,4)
@time g, ginv = rtfbilin(type,a=val[1],b=val[2],c=val[3],d=val[4])
@time syst, gi = gbilin(sys,g,atol=1.e-7);
@time sysi, g1 = gbilin(syst,ginv,atol=1.e-7);
@test iszero(sys-sysi,atol=1.e-7) && gi == ginv && g == g1

type = "lft";
sys = rss(T = Ty, 3,2,3,disc=true);
val = rand(Ty,4)
@time g, ginv = rtfbilin(type,Ts=1, Tsi=1, a=val[1],b=val[2],c=val[3],d=val[4])
@time syst, gi = gbilin(sys,g,atol=1.e-7);
@time sysi, g1 = gbilin(syst,ginv,atol=1.e-7);
@test iszero(sys-sysi,atol=1.e-7) && gi == ginv && g == g1

type = "lft";
sys = rss(T = Ty, 3,2,3,disc=true);
val = rand(Ty,4)
@time g, ginv = rtfbilin(type,Ts=0, Tsi=-1, a=val[1],b=val[2],c=val[3],d=val[4])
@time syst, gi = gbilin(sys,g,atol=1.e-7);
@time sysi, g1 = gbilin(syst,ginv,atol=1.e-7);
@test iszero(sys-sysi,atol=1.e-7) && gi == ginv && g == g1

type = "lft";
sys = rss(T = Ty, 3,2,3,disc=false);
val = rand(Ty,4)
@time g, ginv = rtfbilin(type,Ts=-1, Tsi=0, a=val[1],b=val[2],c=val[3],d=val[4])
@time syst, gi = gbilin(sys,g,atol=1.e-7);
@time sysi, g1 = gbilin(syst,ginv,atol=1.e-7);
@test iszero(sys-sysi,atol=1.e-7) && gi == ginv && g == g1

type = "lft";
sys = rss(T = Ty, 3,2,3,disc=false);
val = rand(Ty,2)
@time g, ginv = rtfbilin(type,Ts=-1, Tsi=0, a=val[1],b=val[2])
@time syst, gi = gbilin(sys,g,atol=1.e-7);
@time sysi, g1 = gbilin(syst,ginv,atol=1.e-7);
@test iszero(sys-sysi,atol=1.e-7) && gi ≈ ginv && g ≈ g1

type = "lft";
sys = rdss(T = Ty, 3,2,3);
val = rand(Ty,4)
g, ginv = rtfbilin(type,a=val[1],b=val[2],c=val[3],d=val[4])
syst, gi = gbilin(sys,g,atol=1.e-7);
sysi, g1 = gbilin(syst,ginv,atol=1.e-7);
@test iszero(sys-sysi,atol=1.e-7) && gi == ginv && g == g1

type = "lft";
sys = rdss(T = Ty, 3,2,3);
val = rand(Ty,2)
@time g, ginv = rtfbilin(type,a=val[1],b=val[2])
@time syst, gi = gbilin(sys,g,atol=1.e-7);
@time sysi, g1 = gbilin(syst,ginv,atol=1.e-7);
@test iszero(sys-sysi,atol=1.e-7) && gi ≈ ginv && g ≈ g1

type = "lft";
sys = rdss(T = Ty, 3,2,3,disc=true);
val = rand(Ty,4)
@time g, ginv = rtfbilin(type,Ts=1, Tsi=1, a=val[1],b=val[2],c=val[3],d=val[4])
@time syst, gi = gbilin(sys,g,atol=1.e-7)
@time sysi, g1 = gbilin(syst,ginv,atol=1.e-7);
@test iszero(sys-sysi,atol=1.e-7) && gi == ginv && g == g1

type = "lft";
sys = rdss(T = Ty, 3,2,3,disc=true);
val = rand(4)
@time g, ginv = rtfbilin(type,Ts=0, Tsi=-1, a=val[1],b=val[2],c=val[3],d=val[4])
@time syst, gi = gbilin(sys,g,atol=1.e-7);
@time sysi, g1 = gbilin(syst,ginv,atol=1.e-7);
@test iszero(sys-sysi,atol=1.e-7) && gi == ginv && g == g1

type = "lft";
sys = rdss(T = Ty, 3,2,3,disc=false);
val = rand(Ty,4)
@time g, ginv = rtfbilin(type,Ts=-1, Tsi=0, a=val[1],b=val[2],c=val[3],d=val[4])
@time syst, gi = gbilin(sys,g,atol=1.e-7);
@time sysi, g1 = gbilin(syst,ginv,atol=1.e-7);
@test iszero(sys-sysi,atol=1.e-7) && gi == ginv && g == g1

end #Ty
end # gbilin

@testset "c2d - descriptor systems" begin

a = [-4 -2;1 0]; b = [2;0]; c = [0.5 1]; d = [0]; x0 = [1,2]; u0 = [1];
sysc = dss(a,b,c,d);
@time sysd, xd0, Mx, Mu = c2d(sysc,1; x0 = [1,2], u0 = [1], state_mapping = true);
EAt = exp(a);
@test dcgain(sysd) ≈ dcgain(sysc) && sort(gpole(sysd)) ≈ sort(exp.(gpole(sysc))) &&
      sysd.A ≈ EAt && sysd.B ≈ a\(EAt-I)*b && x0 ≈ Mx*xd0+Mu*u0
@time sysd_foh, xd0, Mx, Mu = c2d(sysc,1,"foh"; x0, u0, state_mapping = true);
@test dcgain(sysd_foh) ≈ dcgain(sysc) && x0 ≈ Mx*xd0+Mu*u0
@time sysd_imp, xd0, Mx, Mu = c2d(sysc,1,"impulse"; x0, u0, state_mapping = true);
@test sort(gpole(sysd_imp)) ≈ sort(exp.(gpole(sysc))) && x0 ≈ Mx*xd0+Mu*u0
@time sysd1, xd1, Mx1, Mu1 = c2d(sysc,1,"tustin"; x0, u0, standard = true, state_mapping = true)
@time sysd2, xd2, Mx2, Mu2 = c2d(sysc,1,"tustin"; x0, u0, standard = false, state_mapping = true)
syst = gbilin(sysc,rtfbilin("Tustin"; Ts = 1)[1])[1]
@test iszero(sysd1 - sysd2, atol=1.e-7) && iszero(sysd1-syst, atol=1.e-7) &&
      xd1-xd2 ≈ (sysd2.E-I)*x0 && x0 ≈ [Mx1 Mu1]*[xd1;u0] && x0 ≈ [Mx2 Mu2]*[xd2;u0]

@time sysd1, xd1, Mx1, Mu1 = c2d(sysc,1,"tustin"; prewarp_freq=1, x0, u0, standard = true, state_mapping = true)
@time sysd2, xd2, Mx2, Mu2 = c2d(sysc,1,"tustin"; prewarp_freq=1, x0, u0, standard = false, state_mapping = true)
syst = gbilin(sysc,rtfbilin("Tustin"; Ts = 1, prewarp_freq=1)[1])[1]
@test gnrank(sysd1 - sysd2, atol=1.e-7) == 0 && iszero(sysd1-syst, atol=1.e-7) && xd1-xd2 ≈ (sysd2.E-I)*x0 &&
      evalfr(sysc,fval=1) ≈ evalfr(sysd1,fval=1) && evalfr(sysc,fval=1) ≈ evalfr(sysd2,fval=1)  &&
      x0 ≈ [Mx1 Mu1]*[xd1;u0] && x0 ≈ [Mx2 Mu2]*[xd2;u0]

@test_throws ErrorException c2d(sysc,1,"Euler")

a = [-4 -2;1 0]; b = [2;0]; c = [0.5 1]; d = [0]; e = [1 2; 3 0]; x0 = [1,2]; u0 = [1];
sysc = dss(a,e,b,c,d);
@time sysd, xd0, Mx, Mu = c2d(sysc,1; x0, u0, state_mapping = true);
a1, e1, b1, c1, d1 = dssdata(sysc);
EAt = exp(e1\a1);
@test dcgain(sysd) ≈ dcgain(sysc) && sort(gpole(sysd)) ≈ sort(exp.(gpole(sysc))) &&
      sysd.A ≈ EAt && sysd.B ≈ (e1\a1)\(EAt-I)*(e1\b1) && x0 ≈ Mx*xd0+Mu*u0

a = [-4 -2;1 0]; b = [2;0]; c = [0.5 1]; d = [0]; e = [1 2; 0 0]; x0 = [1,2]; u0 = [1];
sysc = dss(a,e,b,c,d);
@time sysd, xd0, Mx, Mu = c2d(sysc,1; x0, u0, state_mapping = true, simple_infeigs = false);
@test norm(Mx*xd0+Mu*u0 - x0) > 0.001

x0 = Mx*xd0+Mu*u0;
@time sysd, xd1, Mx, Mu = c2d(sysc,1; x0, u0, state_mapping = true, simple_infeigs = false);
@test norm(Mx*xd1+Mu*u0 - x0) < 0.001

@time sysd, xd0, Mx, Mu = c2d(sysc,1; x0, u0, state_mapping = true, simple_infeigs = true);
@test norm(Mx*xd0+Mu*u0 - x0) > 0.001

x0 = Mx*xd0+Mu*u0;
@time sysd, xd1, Mx, Mu = c2d(sysc,1; x0, u0, state_mapping = true, simple_infeigs = true);
@test norm(Mx*xd1+Mu*u0 - x0) < 0.001

a = [-4 -2;1 0]; b = [2;1]; c = [0.5 1]; d = [0]; e = [1 2; 0 0]; x0 = [1,2]; u0 = [1];
sysc = dss(a,e,b,c,d);
@time sysd, xd0, Mx, Mu = c2d(sysc,1; x0, u0, state_mapping = true);
a1, e1, b1, c1, d1 = dssdata(gir(sysc,contr=false,obs=false,noseig=true));
EAt = exp(e1\a1);
@test dcgain(sysd) ≈ dcgain(sysc) && sort(gpole(sysd)) ≈ sort(exp.(gpole(sysc))) &&
      sysd.A ≈ EAt && sysd.B ≈ (e1\a1)\(EAt-I)*(e1\b1) && x0 ≈ Mx*xd0+Mu*u0
@time sysd_foh, xd0, Mx, Mu = c2d(sysc,1,"foh"; x0, u0, state_mapping = true);
@test dcgain(sysd_foh) ≈ dcgain(sysc) && x0 ≈ Mx*xd0+Mu*u0
@time sysd_imp, xd0, Mx, Mu = c2d(sysc,1,"impulse"; x0, u0, state_mapping = true);
@test sort(gpole(sysd_imp)) ≈ sort(exp.(gpole(sysc))) && x0 ≈ Mx*xd0+Mu*u0

@time sysd1, xd1, Mx1, Mu1 = c2d(sysc,1,"tustin"; x0, u0, standard = true, state_mapping = true)
@time sysd2, xd2, Mx2, Mu2 = c2d(sysc,1,"tustin"; x0, u0, standard = false, state_mapping = true)
syst = gbilin(sysc,rtfbilin("Tustin"; Ts = 1)[1])[1]
@test iszero(sysd1 - sysd2, atol=1.e-7) && iszero(sysd1-syst, atol=1.e-7) &&
      xd1-xd2 ≈ (sysd2.E-I)*x0 && x0 ≈ [Mx1 Mu1]*[xd1;u0] && x0 ≈ [Mx2 Mu2]*[xd2;u0]

@time sysd1, xd1, Mx1, Mu1 = c2d(sysc,1,"tustin"; prewarp_freq=1, x0, u0, standard = true, state_mapping = true)
@time sysd2, xd2,Mx2, Mu2 = c2d(sysc,1,"tustin"; prewarp_freq=1, x0, u0, standard = false, state_mapping = true)
syst = gbilin(sysc,rtfbilin("Tustin"; Ts = 1, prewarp_freq=1)[1])[1]
@test iszero(sysd1 - sysd2, atol=1.e-7) && iszero(sysd1-syst, atol=1.e-7) &&
      evalfr(sysc,fval=1) ≈ evalfr(sysd1,fval=1) && evalfr(sysc,fval=1) ≈ evalfr(sysd2,fval=1) &&
      x0 ≈ [Mx1 Mu1]*[xd1;u0] && x0 ≈ [Mx2 Mu2]*[xd2;u0]

a = [-4 -2;1 0]; b = [2;0]; c = [0.5 1]; d = [0]; e = [0 0; 0 0]; x0 = [1,2]; u0 = [1];
sysc = dss(a,e,b,c,d);
@time sysd, xd0, Mx, Mu = c2d(sysc,1; x0, u0, state_mapping = true, simple_infeigs = true);
@test norm(Mx*xd0+Mu*u0 - x0) > 0.001
x0 = Mx*xd0+Mu*u0;
@time sysd, xd1, Mx, Mu = c2d(sysc,1; x0, u0, state_mapping = true, simple_infeigs = true);
@test norm(Mx*xd1+Mu*u0 - x0) < 0.001

end # c2d

@testset "c2d - rational transfer functions" begin

s = rtf('s'); z = rtf('z',Ts=1);
rc = 1/s;

rd1 = 1/(z-1);
@time rd = c2d(rc, 1, "zoh");
@test rd1 ≈ rd

rd2 = 0.5*(z+1)/(z-1);
@time rd = c2d(rc, 1, "foh");
@test rd2 ≈ rd

rd3 = z/(z-1);
@time rd = c2d(rc, 1, "impulse");
@test rd3 ≈ rd

rd4 = 0.5*(z+1)/(z-1);
@time rd = c2d(rc, 1, "tustin");
@test rd4 ≈ rd

rd5 = 1.0001/(z-1);
@time rd = c2d(rc, 1, "matched");
@test norm((rd5-rd).num.coeffs,Inf) < 0.001

rc = 1/(s+1);

z1 = exp(-1);
rd1 = (1-z1)/(z-z1);
@time rd = c2d(rc, 1, "zoh");
@test rd1 ≈ rd

z2 = 1-2*z1;
rd2 = (z1*z+z2)/(z-z1);
@time rd = c2d(rc, 1, "foh");
@test rd2 ≈ rd

rd3 = z/(z-z1);
@time rd = c2d(rc, 1, "impulse");
@test rd3 ≈ rd

rd4 = 1/3*(z+1)/(z-1/3);
@time rd = c2d(rc, 1, "tustin");
@test rd4 ≈ rd

rd5 = (1-z1)/(z-z1);
@time rd = c2d(rc, 1, "matched");
@test rd5 ≈ rd

rc = s+1;

rd4 = (3*z-1)/(z+1);
@time rd = c2d(rc, 1, "tustin");
@test rd4 ≈ rd

z1 = exp(-1);
k1 = 1/(1-z1);
rd5 = k1*(z-z1);
@time rd = c2d(rc, 1, "matched");
@test rd5 ≈ rd

s = rtf('s');
Gc = [s^2 s/(s+1); 0 1/s]
@time Gd = c2d(Gc, 1, "matched");

end

@testset "dss2pm & dss2rm" begin

# Example 3: P = [λ^2 λ; λ 1] DeTeran, Dopico, Mackey, ELA 2009
λ = Polynomial([0,1],:s)
P = [λ^2 λ; λ 1]
@test all(P .≈ dss2pm(dss(P)))
@test all(P .≈ dss2rm(dss(P)))

z = rtf('z')
@test_throws ErrorException dss2pm(dss(z+1/z))

s = Polynomial([0, 1],:s);
R = rtf.([s^2+3*s+3 1; -1 2*s^2+7*s+4] ./ [(s+1)^2 s+2; (s+1)^3 (s+1)*(s+2)]);
@test all(R .≈ dss2rm(dss(R)))

end # dss2pm & dss2rm

@testset "timeresp" begin

a = [-4 -2;1 0]; b = [2 1;0 1]; c = [0.5 1]; d = [0 1]; x0 = [1,2]; u0 = [1, 1]; Ts = 1;
sysc = dss(a,b,c,d);
@time sysd, xd0, = c2d(sysc, Ts; x0, u0);
@time y, tout, x = timeresp(sysd, ones(11,2), Int[], xd0; state_history = true);

@time y1, tout1, x1 = timeresp(sysc, ones(11,2), tout, x0; state_history = true, interpolation = "zoh")
@test norm(y-y1,Inf) < 1.e-7 && norm(x-x1,Inf) < 1.e-7

@time y2, tout2, x2 = timeresp(sysc, ones(11,2), tout, x0; state_history = true, interpolation = "foh")
@test norm(y-y2,Inf) < 1.e-7 && norm(x-x2,Inf) < 1.e-7

ed = rand(2,2);
sysdd = dss(ed*sysd.A,ed,ed*sysd.B,sysd.C,sysd.D,Ts=sysd.Ts)
@time y3, tout3, x3 = timeresp(sysdd, ones(11,2), Int[], xd0; state_history = true)
@test norm(y-y3,Inf) < 1.e-7 && norm(x-x3,Inf) < 1.e-7

u = rand(11,2);
@time sysd, xd0, Mx, Mu = c2d(sysc, Ts, "foh"; x0, u0 = u[1,:], state_mapping = true);

@time y, tout, x = timeresp(sysd, u, Int[], xd0; state_history = true)

@time y2, tout2, x2 = timeresp(sysc, u, tout, x0; state_history = true, interpolation = "foh")
@test norm(y-y2,Inf) < 1.e-7 && norm([Mx Mu]*[x';u']-x2',Inf) < 1.e-7

a = [-4 -2;1 0]; b = [2 1;0 1]; c = [0.5 1]; d = [0 1]; x0 = [1,2]; u0 = [1, 1]; Ts = 1;
e = [1 3; 1 1];
u = ones(11,2);
sysc = dss(a,b,c,d); t = 0:Ts:10;
@time y, tout, x = timeresp(sysc, u, t, x0; state_history = true)

sysc1 = dss(e*a,e,e*b,c,d);
@time y1, tout, x1 = timeresp(sysc1, u, t, x0; state_history = true, interpolation = "zoh")
@test norm(y-y1,Inf) < 1.e-7 && norm(x-x1,Inf) < 1.e-7

@time sysd, xd0, Mx, Mu = c2d(sysc, Ts; state_mapping = true, x0, u0);
@time y2, tout, x2 = timeresp(sysd, u, t, xd0; state_history = true) # OK
@test norm(y-y2,Inf) < 1.e-7 && norm(x-x2,Inf) < 1.e-7

sysdd = dss(e*sysd.A,e,e*sysd.B,sysd.C,sysd.D,Ts = 1);
@time yd, tout, xd = timeresp(sysdd, u, t, xd0; state_history = true) # OK
@test norm(y-yd,Inf) < 1.e-7 && norm(x-xd,Inf) < 1.e-7

@time sysd1, xd1, Mx, Mu = c2d(sysc1, Ts; state_mapping = true, x0, u0);
@time y3, tout3, x3 = timeresp(sysd1, u, t, xd1; state_history = true) # OK
@test norm(y-y3,Inf) < 1.e-7 && norm([Mx Mu]*[x3';u']-x',Inf) < 1.e-7

@time y4, tout, x4 = timeresp(sysc, u, t, x0; state_history = true, interpolation = "foh")
@test norm(y-y4,Inf) < 1.e-7 && norm(x-x4,Inf) < 1.e-7

@time y5, tout, x5 = timeresp(sysc1, u, t, x0; state_history = true, interpolation = "foh")
@test norm(y-y5,Inf) < 1.e-7 && norm(x-x5,Inf) < 1.e-7

@time sysd2, xd2, Mx, Mu = c2d(sysc, Ts, "foh"; state_mapping = true, x0, u0);
@time y6, tout, x6 = timeresp(sysd2, u, t, xd2; state_history = true) # OK
@test norm(y-y6,Inf) < 1.e-7 && norm([Mx Mu]*[x6';u']-x',Inf) < 1.e-7

sysdd2 = dss(e*sysd2.A,e,e*sysd2.B,sysd2.C,sysd2.D,Ts = 1);
@time ydd2, tout, xdd2 = timeresp(sysdd2, u, t, xd2; state_history = true) # OK
@test norm(y-ydd2,Inf) < 1.e-7 && norm([Mx Mu]*[xdd2';u']-x',Inf) < 1.e-7

@time sysd3, xd3, Mx, Mu = c2d(sysc1, Ts, "foh"; state_mapping = true, x0, u0);
@time y7, tout, x7 = timeresp(sysd3, u, t, xd3; state_history = true) # OK
@test norm(y-y7,Inf) < 1.e-7 && norm([Mx Mu]*[x7';u']-x',Inf) < 1.e-7

u = rand(11,2);
@time sysd, xd0, Mx, Mu = c2d(sysc, Ts, "foh"; x0, u0 = u[1,:], state_mapping = true);

@time y, tout, x = timeresp(sysd, u, Int[], xd0; state_history = true)

@time y2, tout2, x2 = timeresp(sysc, u, tout, x0; state_history = true, interpolation = "foh")
@test norm(y-y2,Inf) < 1.e-7 && norm([Mx Mu]*[x';u']-x2',Inf) < 1.e-7

a = [-4 -2;1 0]; b = [2 1;0 1]; c = [0.5 1]; d = [0 1]; x0 = [1,2]; u0 = [1, 1]; Ts = 1;
sysc = dss(a,b,c,d); t = 0:Ts:10;
u = ones(11,2);
@time y, tout, x = timeresp(sysc, u, t, x0; state_history = true);
e = rand(Complex{Float64},2,2);
syscc = dss(e*a,e,e*b,c,d); t = 0:Ts:10;
@time sysdc, xd0, Mx, Mu = c2d(syscc, Ts; x0, u0, state_mapping = true);
@time yc, tout, xc = timeresp(sysdc, u, t, xd0; state_history = true)
@test norm(y-yc,Inf) < 1.e-7 && norm([Mx Mu]*[xc';u']-x',Inf) < 1.e-7

@time yc1, tout, xc1 = timeresp(syscc, u, t, x0; state_history = true)
@test norm(y-yc1,Inf) < 1.e-7 && norm(x-xc1,Inf) < 1.e-7

@time yc2, tout1, xc2 = timeresp(syscc, u, t, x0; state_history = true, interpolation = "foh")
@test norm(y-yc2,Inf) < 1.e-7 && norm(x-xc2,Inf) < 1.e-7

u = exp.(im*rand(11,2));
@time y, tout, x = timeresp(sysc, u, t, x0; state_history = true, interpolation = "foh");

@time sysdc, xd0, Mx, Mu = c2d(syscc, Ts, "foh"; x0, u0 = u[1,:], state_mapping = true);

@time yc1, tout, xc1 = timeresp(sysdc, u, t, xd0; state_history = true)
@test norm(y-yc1,Inf) < 1.e-7 && norm([Mx Mu]*[xc1';u']-x',Inf) < 1.e-7

@time yc2, tout2, xc2 = timeresp(syscc, u, tout, x0; state_history = true, interpolation = "foh")
@test norm(y-yc2,Inf) < 1.e-7 && norm(x-xc2,Inf) < 1.e-7

a = [-4 -2;1 0]; b = [2 1;0 1]; c = [0.5 1]; d = [0 1]; x0 = [1,2]; u0 = [1, 1]; Ts = 1;
e = [1 2; 0 0];
sysc = dss(a,e,b,c,d);
@time sysd, xd0, Mx, Mu = c2d(sysc, Ts; x0, u0, state_mapping = true, atol=1.e-7);

@time y, tout, x = timeresp(sysd, ones(11,2), Int[],xd0; state_history = true, atol=1.e-7)

@time y1, tout1, x1 = timeresp(sysc, ones(11,2), tout, x0; state_history = true, interpolation = "zoh", atol=1.e-7)
@test norm(y-y1,Inf) < 1.e-7

@time y2, tout2, x2 = timeresp(sysc, ones(11,2), tout, x0; state_history = true, interpolation = "foh", atol=1.e-7)
@test norm(y-y2,Inf) < 1.e-7

sysd = rdss(3,2,2,disc=true,stable=true,id=ones(Int,2))
sysdr = gss2ss(sysd)[1]
u = rand(11,2);
x0 = rand(5);

@time y, tout, x = timeresp(sysd, u, Int[], x0; state_history = true, atol=1.e-7)
sys, xt, Mx, Mu = dss2ss(sysd, x0; state_mapping = true)
@time y1, tout1, x1 = timeresp(sys, u, Int[], xt; state_history = true, atol=1.e-7)
@test norm(y-y1,Inf) < 1.e-7 && [Mx Mu]*[x1';u'] ≈ x'

@time y2, tout2, x2 = timeresp(sys, u, Int[]; state_history = false, atol=1.e-7)
@time y3, tout3, x3 = timeresp(sysdr, u, Int[]; state_history = false, atol=1.e-7)
@test norm(y2-y3,Inf) < 1.e-7

end # timeresp

end #module


text:
# First-order irreversible series reaction
id = :Sahlodin2011a

pL = [-1.0]
pU = [1.0]
Sahlodin2011a_x0(p::Vector{T}) where T = T[9.0*one(T)]
function Sahlodin2011a_f!(du::Vector{T}, u::Vector{T}, p::Vector{T}, t) where T
    du[1] = -u[1]^2 + p[1]
    return
end
tspan = (0.0, 1.0)

prob = ODERelaxProb(Sahlodin2011a_f!, tspan, Sahlodin2011a_x0, pL, pU)
url = "https://www.sciencedirect.com/science/article/abs/pii/S0168927411000316"
source = "Sahlodin2011"
desc = "First-order irreversible series reaction"

ext: jl


text:
"""
Acquisition geometry for each supersource.
"""
mutable struct AGeomss
	sx::Vector{Float64}
	sz::Vector{Float64}
	rx::Vector{Float64}
	rz::Vector{Float64}
	ns::Int64 # change to +ve only integer later
	nr::Int64 # change to +ve only integer later
	"adding conditions that are to be false while construction"
	AGeomss(sx, sz, rx, rz, ns, nr) =
		any([
                  length(sx) != ns,
                  length(sz) != ns,
                  length(rx) != nr,
                  length(rz) != nr,
		  ]) ?
		error("AGeomss construct") : new(sx, sz, rx, rz, ns, nr)
end # type

include("gallery.jl")

"""
Randomly place a given number of source and receivers in a `mgrid`.
"""
function AGeomss(mgrid::AbstractArray{T}, s::Srcs, r::Recs) where {T<:StepRangeLen}
	xmin=mgrid[2][1]
	xmax=mgrid[2][end]
	zmin=mgrid[1][1]
	zmax=mgrid[1][end]
	return AGeomss(Random.rand(Uniform(xmin,xmax),s.n), Random.rand(Uniform(zmin,zmax),s.n),
		    Random.rand(Uniform(xmin,xmax),r.n), Random.rand(Uniform(zmin,zmax),r.n), s.n, r.n)
end

"""
```julia
ageom=AGeom(mgrid, SSrcs(2), [Srcs(3), Srcs(2)], [Recs(10), Recs(20)])
```

Initialize an acquisition, with 2 supersources, by randomly placing a
sources and receivers in `mgrid`.
"""
function Array{AGeomss,1}(mgrid::AbstractVector{T}, ss::SSrcs, s::Vector{Srcs}, r::Vector{Recs}) where {T<:StepRangeLen}
	return [AGeomss(mgrid,s[i],r[i]) for i in 1:ss.n]
end

"""
```julia
ageom=AGeom(mgrid, SSrcs(10), Srcs(10), Recs(20))
```
Initialize an acquisition, with 10 supersources, by randomly placing a
sources and receivers in `mgrid`. All the supersources have same source and receiver positions.
"""
function Array{AGeomss,1}(mgrid::AbstractVector{T}, ss::SSrcs, s::Srcs, r::Recs)  where {T<:StepRangeLen}
	return Array{AGeomss,1}(mgrid, ss, fill(s,ss.n), fill(r,ss.n))
end

function Base.isequal(ageom1::AGeomss, ageom2::AGeomss, attrib=:all)
	srcfnames = [:sx, :sz, :ns]
	recfnames = [:rx, :rz, :nr]
	srcvec=[(isequal(getfield(ageom1, name),getfield(ageom2, name))) for name in srcfnames]
	recvec=[(isequal(getfield(ageom1, name),getfield(ageom2, name))) for name in recfnames]
	if(attrib == :sources)
		return all(srcvec)
	elseif(attrib == :receivers)
		return all(recvec)
	else
		return all(vcat(srcvec, recvec))
	end
end

"""
```julia
isequal(ageom1, ageom2)
```
Assert if `ageom1` equals `ageom2`.
"""
function Base.isequal(ageom1::AGeom, ageom2::AGeom, attrib=:all)
	return (length(ageom1)==length(ageom2)) ? all([isequal(ageom1[iss],ageom2[iss]) for iss in 1:length(ageom1)]) : false
end

function Base.show(io::Base.IO, ageomss::AGeomss)
	print(io, "\t> acquisition with ",ageomss.ns," source(s) and ",ageomss.nr, " receiver(s)")
end

"""
Output `n` interpolated 2-D points between `p1` and `p2`.
"""
function spread(n, p1::Vector{T1}, p2::Vector{T2}, rand_flag::Bool=false) where {T1<:Real,T2<:Real}
	@assert (length(p1)==2) && (length(p2)==2)

	if(n==1)
		return [p1[1]], [p1[2]]
	else
		knots = ([1,n],)
		itpz = Interpolations.interpolate(knots, [p1[1],p2[1]],  Gridded(Linear()));
		itpx = Interpolations.interpolate(knots, [p1[2],p2[2]],  Gridded(Linear()));
		return [itpz(i) for i in 1:n], [itpx(i) for i in 1:n]
	end
end

function spread(n, p0::Vector{T1}, rad::T2, angles::Vector{T3}, rand_flag::Bool=false)  where {T1<:Real,T2<:Real,T3<:Real}
	th,ni=spread(n, [angles[1], 1],[angles[2], n], rand_flag)
	x=zeros(n); z=zeros(n)
	for i in 1:n
		x[i]=p0[2]+rad*cos(th[i])
		z[i]=p0[1]+rad*sin(th[i])
	end
	return z,x
end

"""
```julia
update!(ageom[1], Srcs(10), p1, p2)
```
Update source geometry for the first supersource, with 10 sources placed regularly between `p1=[z1,x1]` and `p2=[z2,x2]`

```julia
update!(ageom[2], Srcs(10), p0, rad, angles)
```
Update source geometry for the second supersource, with 10 sources placed equidistant from `p0=[z0,x0]` given radius `rad`
and angles `angles=[0,2pi]`.
"""
function update!(ageomss::AGeomss, ::Srcs, args...) where {T<:Real}
	z,x=spread(ageomss.ns,args...)
	copyto!(ageomss.sx,x)
	copyto!(ageomss.sz,z)
	return ageomss
end

"""
```julia
update!(ageom[1], Recs(10), args...)
```
Similar to updating source positions, but for receivers.
"""
function update!(ageomss::AGeomss, ::Recs, args...)
	z,x=spread(ageomss.nr,args...)
	copyto!(ageomss.rx,x)
	copyto!(ageomss.rz,z)
	return ageomss
end

"""
```julia
update!(ageom, SSrcs(), p1, p2)
update!(ageom, SSrcs(), p0, rad, angles)
```
Update source positions for all supersources.
"""
function update!(ageom::AGeom, ::SSrcs, args...)
	z,x=spread(length(ageom), args...)
	for iss in 1:length(ageom)
		update!(ageom[iss], Srcs(), [z[iss],x[iss]], [z[iss],x[iss]])
	end
	return ageom
end

"""
```julia
update!(ageom, Recs(), p1, p2)
update!(ageom, Recs(), p0, rad, angles)
```
Update receiver positions for all supersources.
"""
function update!(ageom::AGeom, ::Recs, args...)
	for iss in 1:length(ageom)
		update!(ageom[iss], Recs(), args...)
	end
	return ageom
end

"""
```julia
ageom ∈ mgrid
in(ageom, mgrid)
```
Assert if all the sources and receivers in `AGeom` are within bounds of `mgrid`.
"""
function Base.in(ageom::AGeom, mgrid::AbstractVector{T}) where {T<:StepRangeLen}
	xmin, zmin, xmax, zmax = mgrid[2][1], mgrid[1][1], mgrid[2][end], mgrid[1][end]
	checkvec=fill(false, length(ageom))
	for iss=1:length(ageom)
		checkvec[iss] = any(
		      vcat(
			   (((ageom[iss].sx .- xmin).*(xmax .- ageom[iss].sx)) .< 0.0),
			   (((ageom[iss].sz .- zmin).*(zmax .- ageom[iss].sz)) .< 0.0),
			   (((ageom[iss].rx .- xmin).*(xmax .- ageom[iss].rx)) .< 0.0),
			   (((ageom[iss].rz .- zmin).*(zmax .- ageom[iss].rz)) .< 0.0),
			   isnan.(ageom[iss].sx),
			   isnan.(ageom[iss].sz),
			   isnan.(ageom[iss].rx),
			   isnan.(ageom[iss].rz),
		     ) )
	end
	return !(any(checkvec))
end

"""
```julia
A=SparseArrays.SparseMatrixCSC(ageom[1],mgrid)
D=A*P
```
Return a sparse matrix that restricts the field `P` on `mgrid` to receiver positions, to give `D`.
"""
function SparseArrays.SparseMatrixCSC(ageomss::AGeomss,mgrid::AbstractVector{T}) where {T<:StepRangeLen}
	nz,nx=length.(mgrid)
	ACQ=spzeros(ageomss.nr,prod(length.(mgrid)))
	for ir = 1:ageomss.nr
		irx=argmin(abs.(mgrid[2].-ageomss.rx[ir]))
		irz=argmin(abs.(mgrid[1].-ageomss.rz[ir]))
		ACQ[ir,irz+(irx-1)*nz]=1.0
	end
	return ACQ
end

#=

"""
Advance either source or receiver array in an acquisition ageometry in horizontal or vertical directions.

# Arguments
* `ageom::AGeom` : acquisition ageometry that is updated
* `advances::Vector{Float64}=[[0.,0.], [0.,0.,]]` : source and receiver advancements
"""
function AGeom_advance(ageom::AGeom, advances::Vector{Vector{Float64}}=[[0.,0.], [0.,0.]])
	ageom_out=deepcopy(ageom)
	for iss=1:ageom.nss
		ageom_out.sx[iss][:] .+= advances[1][2]
		ageom_out.sz[iss][:] .+= advances[1][1]
		ageom_out.rx[iss][:] .+= advances[2][2]
		ageom_out.rz[iss][:] .+= advances[2][1]
	end
	return ageom_out
end

"""
return a vector of the order
"""
function AGeom_getvec(ageom::Vector{AGeom}, field::Symbol)
	np = length(ageom);
	vect = [getfield(ageom[ip],field) for ip=1:np]
	return vec(vcat(vcat(vect...)...))
end

function save(ageom, folder)
	!(isdir(folder)) && error("invalid directory")
	for iss in 1:ageom.nss
		file=joinpath(folder, string("s", iss, ".csv"))
		CSV.write(file,DataFrame(hcat(ageom.sx[iss], ageom.sz[iss])))
		file=joinpath(folder, string("r", iss, ".csv"))
		CSV.write(file,DataFrame(hcat(ageom.rx[iss], ageom.rz[iss])))
	end
end

"""
Appends the input a vector of acquisition ageometries.
"Merge" `AGeom` objects in `ageomv` to `ageom1`.
Typically you have one acquisition,
you change it and want to create the macro acquisition
where the supersources of both acquisitions are appended.
"""
function AGeom_add!(ageom1::AGeom, ageomv::Vector{AGeom})
	for iageom=1:length(ageomv)
		append!(ageom1.sx, ageomv[iageom].sx)
		append!(ageom1.sz, ageomv[iageom].sz)
		append!(ageom1.rx, ageomv[iageom].rx)
		append!(ageom1.rz, ageomv[iageom].rz)
		append!(ageom1.ns, ageomv[iageom].ns)
		append!(ageom1.nr, ageomv[iageom].nr)
		ageom1.nss += ageomv[iageom].nss
	end
	return ageom1
end

"""
Adds input positions as either sources or receivers of every supershot.
"""
function AGeom_add!(ageom::AGeom, zpos::Vector{Float64}, xpos::Vector{Float64}, attrib::Symbol)
	length(xpos) == length(zpos) ? nothing : error("same length vectors needed")
	for iss=1:ageom.nss
		if(attrib == :receivers)
			ageom.rx[iss] = vcat(ageom.rx[iss][:], xpos)
			ageom.rz[iss] = vcat(ageom.rz[iss][:], zpos)
			ageom.nr[iss] += length(xpos)
		elseif(attrib == :sources)
			ageom.sx[iss] = vcat(ageom.sx[iss][:], xpos)
			ageom.sz[iss] = vcat(ageom.sz[iss][:], zpos)
			ageom.ns[iss] += length(xpos)
		end
	end
end

"""
A fixed spread acquisition has same set of sources and
receivers for each supersource.
This method constructs a
fixed spread acquisition ageometry using either a
horizontal or vertical array of supersources/ receivers.
Current implementation has only one source for every supersource.

# Arguments

* `ssmin::Float64` : minimum coordinate for sources
* `ssmax::Float64` : maximum coordinate for sources
* `ss0::Float64` : consant coordinate for sources
* `rmin::Float64` : minimum coordinate for receivers
* `rmax::Float64` : maximum coordinate for receivers
* `r0::Float64` : consant coordinate for receivers
* `nss::Int64` : number of supersources
* `nr::Int64` : number of receivers
* `ssattrib::Symbol=:horizontal` : supersource array kind
  `=:vertical` : vertical array of supersources
  `=:horizontal` horizontal array of supersources
* `rattrib::Symbol=:horizontal` : receiver array kind
  `=:vertical` : vertical array of receivers
  `=:horizontal` horizontal array of receivers
* `rand_flags::Vector{Bool}=[false, false]` : decide placement of supersources and receivers
  `=[true, false]` : randomly place supersources for regularly spaced receivers
  `=[true, true]` : randomly place supersources and receivers
  `=[false, false]` : regularly spaced supersources and receivers
  `=[false, true]` : randomly place receivers for regularly spaced supersources

# Return
* a fixed spread acquisition ageometry `AGeom`
"""
function AGeom_fixed(
	      ssmin::Real,
	      ssmax::Real,
	      ss0::Real,
	      rmin::Real,
	      rmax::Real,
	      r0::Real,
	      nss::Int64,
	      nr::Int64,
	      ssattrib::Symbol=:horizontal,
	      rattrib::Symbol=:horizontal,
 	      rand_flags::Vector{Bool}=[false, false];
	      ssα::Real=0.0,
	      rα::Real=0.0,
	      ns::Vector{Int64}=ones(Int,nss),
	      srad::Real=0.0
	     )

	ssmin=Float64(ssmin); rmin=Float64(rmin)
	ssmax=Float64(ssmax); rmax=Float64(rmax)
	ss0=Float64(ss0); r0=Float64(r0)
	ssα=Float64(ssα)*pi/180.
	rα=Float64(rα)*pi/180.

	ssarray = isequal(nss,1) ? fill(ssmin,1) : (rand_flags[1] ?
					     Random.rand(Uniform(minimum([ssmin,ssmax]),maximum([ssmin,ssmax])),nss) : range(ssmin,stop=ssmax,length=nss))
	if(ssattrib==:horizontal)
		ssz = ss0.+(ssarray.-minimum(ssarray)).*sin(ssα)/cos(ssα); ssx=ssarray
	elseif(ssattrib==:vertical)
		ssx = ss0.+(ssarray.-minimum(ssarray)).*sin(ssα)/cos(ssα); ssz=ssarray
	else
		error("invalid ssattrib")
	end

	rarray = isequal(nr,1) ? fill(rmin,1) : (rand_flags[2] ?
					  Random.rand(Uniform(minimum([rmin,rmax]),maximum([rmin,rmax])),nr) : range(rmin,stop=rmax,length=nr))
	if(rattrib==:horizontal)
		rz = r0.+(rarray.-minimum(rarray)).*sin(rα)/cos(rα); rx = rarray
	elseif(rattrib==:vertical)
		rx = r0.+(rarray.-minimum(rarray)).*sin(rα)/cos(rα); rz = rarray
	else
		error("invalid rattrib")
	end
	rxall = [rx for iss=1:nss];
	rzall = [rz for iss=1:nss];
	ssxall=[zeros(ns[iss]) for iss=1:nss];
	sszall=[zeros(ns[iss]) for iss=1:nss];
	for iss in 1:nss
		for is in 1:ns[iss]
			θ=Random.rand(Uniform(-Float64(pi),Float64(pi)))
			r = iszero(srad) ? 0.0 : Random.rand(Uniform(0,srad))
			x=r*cos(θ)
			z=r*sin(θ)
			ssxall[iss][is]=ssx[iss]+x
			sszall[iss][is]=ssz[iss]+z
		end
	end
	return AGeom(ssxall, sszall, rxall, rzall, nss, ns, fill(nr,nss))
end

=#


text:
module CentroidalTrajOpt

export
    CentroidalTrajectoryProblem,
    ContactRegion,
    solve!,
    normals,
    disallow_jumping!,
    ObjectiveTypes

# Re-exports from other packages
export
    SVector,
    RotMatrix,
    AffineMap,
    with_optimizer,
    derivative,
    breaks

# Should move
export
    map_subfunctions

# Visualization
export
    CentroidalTrajectoryVisualizer,
    set_com_trajectory!,
    set_state!,
    save_result,
    load_result

using JuMP
using LinearAlgebra
using StaticArrays
using StaticUnivariatePolynomials
using CoordinateTransformations
using Rotations
using AxisArrays

import QPControl

using QPControl.Trajectories: BezierCurve, Piecewise, Constant
using QPControl.Trajectories: breaks
using StaticUnivariatePolynomials: derivative

const SUP = StaticUnivariatePolynomials

import Polyhedra
using Polyhedra: hrep, vrep, polyhedron, Mesh

import SCIP

include("util.jl")
include("region.jl")
include("problem.jl")
include("result.jl")

using GeometryTypes: Point, Vec, HyperSphere, HyperRectangle, GLNormalMesh

import MeshCat
using MeshCat: Visualizer
using MeshCat: RGB, RGBA
using MeshCat: ArrowVisualizer, HyperSphere, PointCloud, LineSegments, Cone
using MeshCat: MeshLambertMaterial, MeshBasicMaterial, Texture, PngImage, LineBasicMaterial
using MeshCat: setobject!, settransform!, setprop!, setvisible!
using MeshCat: Animation, atframe, setanimation!

include("visualization.jl")
include("serialization.jl")

using .Serialization

end # module


text:
module PageDown

export Page, next!, prev!, go!, first!, last!
export clear, proper, imgcat
export @page, @current, @next, @prev, @go, @first, @last, @step, @open

__precompile__(true)

type Page
  pages::AbstractArray
  start::Int
  step::Int
  option::Dict
  count::Int
  Page(pages::AbstractArray) = new(pages, 1, 1, Dict(), length(pages))
  Page(pages::AbstractArray, start::Int, step::Int) = new(pages, start, step, Dict(), length(pages))
  Page(pages::AbstractArray, start::Int, step::Int, option::Dict) = new(pages, start, step, option, length(pages))
  Page(pages::AbstractArray, option::Pair...) = new(pages, 1, 1, Dict(option), length(pages))
  Page(pages::AbstractArray, option::Dict) = new(pages, 1, 1, option, length(pages))
end

function next!(page::Page)
  if (page.start + page.step) <= page.count
    page.start += page.step
  else
    #println("last page")
  end
  page
end

function prev!(page::Page)
  if page.start > page.step
    page.start -= page.step
  else
    #println("first page")
  end
  page
end

function go!(page::Page, number::Int)
  if number >= 1 && page.count >= number
    page.start = number
  else
    throw(ArgumentError("not found"))
  end
  page
end

function first!(page::Page)
  go!(page, 1)
end

function last!(page::Page)
  if page.count > page.step
    go!(page, page.count - (page.step - page.count % page.step - 1))
  else
    go!(page, page.count)
  end
end

function clear()
  @windows ? nothing : run(`clear`)
end

function proper(page::Page)
  if page.start + page.step <= page.count
    slides = page.pages[page.start:page.start+page.step-1]
  else
    slides = page.pages[page.start:end]
  end
  page,slides
end

function slides(tup::Tuple{Page,Array})
  page,slides = tup
  slides
end

function zzal(tup::Tuple{Page,Array})
  page,slides = tup
  img = haskey(page.option, :IMG) ? page.option[:IMG] : "IMG"
  imgcmd = haskey(page.option, :imgcmd) ? page.option[:imgcmd] : "go run ../images/imgcat.go "
  imgfolder = haskey(page.option, :imgfolder) ? page.option[:imgfolder] : "../images/"
  for slide in slides
    isa(slide, Base.Markdown.MD) && for content in slide.content
      isa(content, Base.Markdown.List) && for item in content.items
        for el in item
          if startswith(el, img)
            zz,image = split(el, "$img ")
            run(`$imgcmd $imgfolder$image`)
          end
        end
      end
    end
  end
  page,slides
end

function open(page::Page, site::Int)
  counter = 0
  opencmd = haskey(page.option, :opencmd) ? page.option[:opencmd] : "open"
  for slide in page |> proper |> slides
    isa(slide, Base.Markdown.MD) && for content in slide.content
      isa(content, Base.Markdown.List) && for item in content.items
        for el in item
          m = match(r"http.*", el)
          if isa(m, RegexMatch)
            counter += 1
            if counter == site
              @osx_only run(`$opencmd $(m.match)`)
              return
            end
          end
        end
      end
    end
  end
end

### macros

macro page()
  p = esc(:page)
  :( StepRange($p.start:$p.step:$p.count) )
end

macro current()
  :( $(esc(:page)) |> proper |> zzal |> slides )
end

macro next()
  clear()
  println()
  :( next!($(esc(:page))) |> proper |> zzal |> slides )
end

macro prev()
  :( prev!($(esc(:page))) |> proper |> slides )
end

macro go(number::Int)
  :( go!($(esc(:page)),$number) |> proper |> zzal |> slides )
end

macro first()
  :( first!($(esc(:page))) |> proper |> zzal |> slides )
end

macro last()
  :( last!($(esc(:page))) |> proper |> zzal |> slides )
end

macro step(pace::Int)
  :( $(esc(:page)).step = $pace )
end

macro open(site::Int)
  :( open($(esc(:page)),$site) )
end

# term
Base.Markdown.term(t::Union{Base.PipeEndpoint,Base.TTY}, a::Array{Base.Markdown.MD,1}) =
  [Base.Markdown.term(t, e) for e in a]
Base.Markdown.term(t::Union{Base.PipeEndpoint,Base.TTY}, a::Array{ASCIIString,1}) = a

end # module


text:
# Algorithm 2.2
# Fixed point iteration, solves g(p) = p

function bisect_approx(f, p0, tolerance, maxIter)
	for i ∈ 1:maxIter
		p = f(p0)
		if abs(p - p0) < tolerance
			return p
		end
		p0 = p
	end
	throw(ArgumentError("Method failed"))
end

println(bisect_approx(x -> ((10/(4 + x))^(1/2)), 1.5, 0.1, 4))


text:
function vikor(decisionMat::DataFrame, weights::Array{Float64,1}; v::Float64=0.5)::VikorResult
    w = unitize(weights)

    nalternatives, ncriteria = size(decisionMat)

    col_max = colmaxs(decisionMat)
    col_min = colmins(decisionMat)

    A = similar(decisionMat)

    for i in 1:nalternatives
        for j in 1:ncriteria
            @inbounds A[i, j] = abs((col_max[j] - decisionMat[i, j]) / (col_max[j] - col_min[j]))
        end
    end

    weightedA = w * A

    s = Array{Float64,1}(undef, nalternatives)
    r = similar(s)
    q = similar(s)

    for i in 1:nalternatives
        s[i] = sum(weightedA[i,:])
        r[i] = maximum(weightedA[i,:])
    end

    smin = minimum(s)
    smax = maximum(s)
    rmin = minimum(r)
    rmax = maximum(r)
    q = v .* ((s .- smin ./ (smax .- smin))) + (1 - v) .* ((r .- rmin ./ (rmax .- rmin)))

    best_index = sortperm(q) |> first

    result = VikorResult(
        decisionMat,
        w,
        weightedA,
        best_index,
        q
    )
    return result
end

max_stars_repo_name: ersagunkuruca/JMcDM

max_issues_repo_name: ersagunkuruca/JMcDM

max_forks_repo_name: ersagunkuruca/JMcDM


text:
# full test of fast and naive validation
using Dates
using ArgParse
using Random

include("../btg.jl")
include("../datasets/load_abalone.jl")
# shuffle data
ind_shuffle = randperm(MersenneTwister(1234), size(data, 1))
data = data[ind_shuffle, :]
target = target[ind_shuffle]

s = ArgParseSettings()
@add_arg_table! s begin
    "--fast"
        help = "use fast or not"
        action = :store_true
end
parsed_args = parse_args(ARGS, s)

if filesize("test_validation7_partial.txt") == 0
    # write the setting headers
    io = open("test_validation7_partial.txt", "w")
    write(io, "Time         ;    ind      ;       θ    ;    λ   ;    x   ;   cov  ;  fast  ;   elapsedmin \n")
    close(io)
end

ind = 1:200
posx = 1:7
posc = 1:3
x = data[ind, posx]
Fx = data[ind, posc]
y = float(target[ind])
trainingData1 = trainingData(x, Fx, y) #training data used for testing various functions
d = getDimension(trainingData1); n = getNumPts(trainingData1); p = getCovDimension(trainingData1)

#prediction data
pind = 1:1
Fx0 = reshape(data[pind, posc], 1, length(posc))
x0 = reshape(data[pind, posx], 1, length(posx))

#parameter setting
rangeθ = [10.0 1000]
rangeθm = repeat(rangeθ, d, 1)
rangeλ = [0. 3.] #we will always used 1 range scale for lambda

btg1 = btg(trainingData1, rangeθm, rangeλ; quadtype = ["SparseGrid", "Gaussian"])
(pdf, cdf, dpdf) = solve(btg1); #initialize training_buffer_dicts, solve once so can use fast techiques to extrapolate submatrix determinants, etc.

"""
leave-one-out-training-data
"""
function lootd(td::AbstractTrainingData, i::Int64)
    x = getPosition(td)
    Fx = getCovariates(td)
    z = getLabel(td)
    x_minus_i = x[[1:i-1;i+1:end], :]
    Fx_minus_i = Fx[[1:i-1;i+1:end], :]
    z_minus_i = z[[1:i-1;i+1:end]]
    x_i = x[i:i, :]
    Fx_i = Fx[i:i, :]
    z_i = z[i:i, :]
    return trainingData(x_minus_i, Fx_minus_i, z_minus_i), x_i, Fx_i, z_i
end

nrow = 6; ncol = 8
y1_temp = Array{Float64, 1}(undef, 100)
y2_temp = Array{Float64, 1}(undef, 100)
y1_set = Array{Float64, 2}(undef, 100, ncol*nrow)
y2_set = Array{Float64, 2}(undef, 100, ncol*nrow)
xgrid = range(.01, stop=1.2, length=100)
before = Dates.now()
for j = 1:nrow*ncol
    mod(j, 10) == 0 ? (@info j) : nothing
    if parsed_args["fast"]
        (pdf, cdf, dpdf) = solve(btg1, validate = j);
        b1 = y0 -> pdf(x0, Fx0, y0);
        c1 = y0 -> cdf(x0, Fx0, y0);
    else
        (trainingdata_minus_i, x_i, Fx_i, z_i) = lootd(trainingData1, j)
        btg2 = btg(trainingdata_minus_i, rangeθm, rangeλ; quadtype = ["SparseGrid", "Gaussian"])
        (pdf_minus_i, cdf_minus_i, dpdf_minus_i) = solve(btg2)
        b1 = y -> pdf_minus_i(x_i, Fx_i, y)
        c1 = y -> cdf_minus_i(x_i, Fx_i, y)
    end
    y1_temp = b1.(xgrid)
    y2_temp = c1.(xgrid)
    if j <= ncol*nrow
        # save for plot
        y1_set[:, j] = y1_temp
        y2_set[:, j] = y2_temp
    end
end
after = Dates.now()
elapsedmin = round(((after - before) / Millisecond(1000))/60, digits=5)

io = open("test_validation7_partial.txt", "a")
write(io, "$(Dates.now())  ;    $ind      ;       $rangeθ    ;    $rangeλ   ;    $posx   ;   $posc  ;  $(parsed_args["fast"])  ;   $elapsedmin \n")
close(io)

## Plot
PyPlot.close("all") #close existing windows
plt, axs = PyPlot.subplots(nrow, ncol)
PyPlot.suptitle("Cross Validation $(parsed_args["fast"])", fontsize=10)
for j = 1:nrow*ncol
    ind1 = Int64(ceil(j/ncol))
    ind2 = Int64(j - ncol*(floor((j-.1)/ncol)))
    axs[ind1, ind2].plot(xgrid, y1_set[:, j], color = "red", linewidth = 1.0, linestyle = "-")
    axs[ind1, ind2].plot(xgrid, y2_set[:, j], color = "orange", linewidth = 1.0, linestyle = "-")
    axs[ind1, ind2].axvline(x = getLabel(btg1.trainingData)[j])
end
# Hide x labels and tick labels for top plots and y ticks for right plots.
for ax in axs
    ax.label_outer()
end
PyPlot.savefig("figure/test_v7_ind_$(ind[1])_$(ind[end])_rθ_$(Int(rangeθ[1]))_$(Int(rangeθ[2]))_rλ_$(Int(rangeλ[1]))_$(Int(rangeλ[2]))_$(parsed_args["fast"]).pdf")

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
ENV["LINES"] = 1000

# %%
n = 100
[(a, b, isqrt(a^2+b^2)) for a in 3:n for b in a+1:a^2÷2
    if gcd(a, b) == 1 && a^2 + b^2 == isqrt(a^2 + b^2)^2]

# %%
n = 100
[(a, b, isqrt(a^2+b^2)) for a in 3:n for b in a+1:max(n, a^2÷2)
    if gcd(a, b) == 1 && a^2 + b^2 == isqrt(a^2 + b^2)^2 && b != a^2÷2]

# %%

max_stars_repo_path: 0020/isqrt example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0020/isqrt example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0020/isqrt example.jl

max_forks_repo_name: genkuroki/public


text:
"""

# read_variational

Read variational sample output files created by cmdstan.

### Method
```julia
read_variational(m::VariationalModelodel)
```

### Required arguments
```julia
* `m::VariationalModel`    : VariationalModel object
```

"""
function read_variational(m::VariationalModel)

  local a3d, index, idx, indvec

  ftype = "chain"

  for i in 1:StanBase.get_n_chains(m)
    if isfile("$(m.output_base)_$(ftype)_$(i).csv")
      instream = open("$(m.output_base)_$(ftype)_$(i).csv")
      skipchars(isspace, instream, linecomment='#')
      line = Unicode.normalize(readline(instream), newline2lf=true)
      idx = split(strip(line), ",")
      index = [idx[k] for k in 1:length(idx)]
      indvec = 1:length(index)
      if i == 1
        a3d = fill(0.0, m.method.output_samples, length(indvec),
          StanBase.get_n_chains(m))
      end
      skipchars(isspace, instream, linecomment='#')
      for j in 1:m.method.output_samples
        skipchars(isspace, instream, linecomment='#')
        line = Unicode.normalize(readline(instream), newline2lf=true)
        if eof(instream) && length(line) < 2
          close(instream)
          break
        else
          flds = parse.(Float64, (split(strip(line), ",")))
          flds = reshape(flds[indvec], 1, length(indvec))
          a3d[j,:,i] = flds
        end
      end
    end
  end

  cnames = convert.(String, idx[indvec])

  (a3d, cnames)

end


text:
using Gadfly, DataArrays, RDatasets

plot(dataset("ggplot2", "diamonds"), x="Price", Geom.histogram(bincount=30))


text:
__precompile__()

module ShowItLikeYouBuildIt

export showarg, summary_build, type_complexity

# This is just a stub package on Julia after the merger of #23389
if VERSION < v"0.7.0-DEV.1790"
    include("code.jl")
else
    include("stub.jl")
end

dimstring(inds) = Base.inds2string(inds)
dimstring(inds::Tuple{Vararg{Base.OneTo}}) = Base.dims2string(map(length, inds))

end # module


text:
using Compat

mutable struct MonomialBasis{T,Degree,NF} <: RefSpace{T,NF} end

valuetype(::MonomialBasis{T}) where {T} = T
degree(::MonomialBasis{T,D}) where {T,D} = D
numfunctions(x::MonomialBasis) = degree(x)+1
function (x::MonomialBasis)(p)
    t = cartesian(p)[1]
    valuetype(x)[t^d for d in 0 : degree(x)]
end

abstract type AbstractTimeBasisFunction end

"""
    TimeBasisFunction{N,D}

T: the value type of the time basis function
N: the number of intervals in the support (this included the semi infinite interval
    stretching to +∞)
D1: the degree of the TBF restricted to each of the intervals **plus one**
"""
mutable struct TimeBasisFunction{T,N,D1,D} <: AbstractTimeBasisFunction
    timestep::T
    numfunctions::Int
    polys::SVector{N,Polynomial{D1,T}}
end

function derive(tbf::TimeBasisFunction{T,N,D1,D}) where {T,N,D1,D}
    dpolys = Polynomial{D1-1,T}[ derive(p) for p in tbf.polys ]
    TimeBasisFunction{T,N,D1-1,D-1}(
        tbf.timestep,
        tbf.numfunctions,
        SVector{N,Polynomial{D1-1,T}}(dpolys...)
    )
end

function integrate(tbf::TimeBasisFunction{T,NI,D1,D}) where {T,NI,D1,D}

    timestep = tbf.timestep
    numfunctions = tbf.numfunctions

    polys = Vector{Polynomial{D1+1,T}}(undef, length(tbf.polys))
    c = zero(timestep)
    for i in 1:length(polys)
        t0 = (i-2) * timestep
        t1 = (i-1) * timestep
        polys[i] = integrate(tbf.polys[i],t0,c)
        c = polys[i](t1)
    end

    TimeBasisFunction{T,NI,D1+1,D+1}(
        tbf.timestep,
        tbf.numfunctions,
        SVector(polys...)
    )
end

mutable struct TimeBasisDelta{T} <: AbstractTimeBasisFunction
    timestep::T
    numfunctions::Int
end

mutable struct DiracBoundary{T} <: RefSpace{T,1} end
numfunctions(x::DiracBoundary) = 1

scalartype(x::TimeBasisDelta{T}) where {T} = T
numfunctions(x::TimeBasisDelta) = x.numfunctions
refspace(x::TimeBasisDelta{T}) where {T} = DiracBoundary{T}()
timestep(x::TimeBasisDelta) = x.timestep
numintervals(x::TimeBasisDelta) = 1

timebasisdelta(dt, Nt)= TimeBasisDelta(dt, Nt)

numintervals(tbf::TimeBasisFunction{T,N,D1,D}) where {T,N,D1,D} = N
degree(tbf::TimeBasisFunction{T,N,D1,D}) where {T,N,D1,D} = D
timestep(tbf::TimeBasisFunction) = tbf.timestep
scalartype(tbf::TimeBasisFunction{T,N,D1,D}) where {T,N,D1,D} = T

numfunctions(t::TimeBasisFunction) = t.numfunctions
refspace(t::TimeBasisFunction{T,N,D1,D}) where {T,N,D1,D} = MonomialBasis{T,D,D1}()

function truncatetail(t::TimeBasisFunction{T,N,D1,D})  where {T,N,D1,D}
    P = Polynomial{D1,T}
    S = SVector{N,P}
    R = SVector{D1,T}

    out = deepcopy(t)
    polys = [p for p in t.polys]
    polys[end] = P(zero(R))
    out.polys = S(polys)
    return out
end

geometry(t::AbstractTimeBasisFunction) = [SegmentedAxis(timestep(t), numfunctions(t)+numintervals(t)-3)]

"""
    timebasisc0d1(type, timestep, numfunctions)

Build the space of continuous, piecewise linear time basis functions. The DoFs
are the time steps. `numfunctions` basis functions will be built in total.
"""
function timebasisc0d1(timestep, numfunctions, T::Type=Float64)
    i, z = one(T), zero(T)
    polys = SVector(
        Polynomial(SVector(i,  i/timestep)), # 1 + t
        Polynomial(SVector(i, -i/timestep)), # 1 - t
        Polynomial(SVector(z,  z))
    )
    TimeBasisFunction{T,3,2,1}(timestep, numfunctions, polys)
end

"""
    timebasiscxd0(timestep, numfunctions, T::Type=Float64)

Create a temporal basis based on shifted copies of the nodal continuous, piecewise
linear interpolant.
"""
function timebasiscxd0(timestep, numfunctions, T::Type=Float64)
    i, z = one(T), zero(T)
    polys = SVector(
        Polynomial(i), # 1
        Polynomial(z), # 0
    )
    TimeBasisFunction{T,2,1,0}(timestep, numfunctions, polys)
end

"""
    timebasisspline2(timestep, numfunctions, T::Type=Float64)

Create a temporal basis based on shifted copies of the quadratic spline. The
spline is the convolution of a cxd0 and a c0d1 basis function.
"""
function timebasisspline2(dt, numfunctions, T::Type=Float64)
    i, z = one(T), zero(T)
    polys = SVector{4,Polynomial{3,T}}(
        Polynomial(SVector(i/2, i/dt, i/2/dt/dt)),
        Polynomial(SVector(i/2, i/dt, -i/dt/dt)),
        Polynomial(SVector(2*i, -2*i/dt, i/2/dt/dt)),
        Polynomial(SVector(z, z, z))
    )
    TimeBasisFunction{T,4,3,2}(dt, numfunctions, polys)
end

function timebasisshiftedlagrange(dt, numfunctions, degree, T::Type=typeof(dt))
    z, i = zero(T), one(T)
    c = Polynomial(SVector(i))
    t = Polynomial(SVector(z,i/dt))
    polys = Polynomial{degree+1,T}[]
    for k = 0:degree
        f = c
        for i in 1:k
            f = (1/i) * f * (i - t)
        end
        g = c
        for i in 1:(degree-k)
            g = (1/i) * g * (i + t)
        end
        push!(polys, f*g)
    end
    push!(polys, 0*Base.power_by_squaring(t, degree))
    @assert length(polys) == degree+2
    polys = SVector{degree+2,Polynomial{degree+1,T}}(polys...)
    TimeBasisFunction{T,degree+2,degree+1,degree}(dt, numfunctions, polys)
end

function (f::TimeBasisFunction)(t::Real)
    dt = timestep(f)
    ni = numintervals(f)
    t < -dt && return zero(t)
    i = floor(Int, (t/dt)) + 2
    p = i <= ni ? f.polys[i] : f.polys[end]
    return p(t)
end

function assemblydata(tbf::TimeBasisFunction)

    T = scalartype(tbf)
    Δt = timestep(tbf)
    z = (0, zero(T))

    t = Polynomial(zero(T), one(T))

    num_cells = numfunctions(tbf)
    num_refs  = degree(tbf)+1

    max_num_funcs = numintervals(tbf)
    numfuncs = zeros(Int, num_cells, num_refs)
    data = fill(z, max_num_funcs, num_refs, num_cells)

    els = [ simplex(point((i-1)*Δt),point(i*Δt)) for i in 1:num_cells ]

    for k in 1 : numfunctions(tbf)
        tk = (k-1) * Δt
        for i in 1 : numintervals(tbf)
            # Focus on interval [(i-2)Δt,(i-1)Δt]
            p = tbf.polys[i]
            q = substitute(p,t-tk)

            c = k + i - 2
            1 <= c <= num_cells || continue
            for d = 0 : degree(q)
                r = d + 1
                w = q[d]

                j = (numfuncs[c,r] += 1)
                data[j,r,c] = (k,w)
            end
        end
    end

    return els, AssemblyData(data)
end

function temporalassemblydata(tbf; kmax=typemax(Int))

    T = scalartype(tbf)
    Δt = timestep(tbf)

    t = Polynomial(zero(T), one(T))

    num_cells = numfunctions(tbf)
    num_refs  = degree(tbf)+1

    has_zero_tail = all(tbf.polys[end].data .== 0)
    # @show has_zero_tail

    if has_zero_tail
        max_num_funcs = numintervals(tbf) # 3
    else
        max_num_funcs = min(kmax, numfunctions(tbf))
    end

    numfuncs = zeros(Int, num_cells, num_refs)
    data = fill((0,zero(T)), max_num_funcs, num_refs, num_cells)
    for k in 1 : min(kmax, numfunctions(tbf))
        tk = (k-1) * Δt
        for i in 1 : numintervals(tbf)-1
            p = tbf.polys[i]
            q = substitute(p,t+tk)

            c = k - i + 1
            c < 1 && continue
            for d = 0 : degree(q)
                r = d + 1
                w = q[d]

                j = (numfuncs[c,r] += 1)
                data[j,r,c] = (k,w)
            end
        end

        # Treat the half-open interval separately
        if !has_zero_tail
            p = tbf.polys[numintervals(tbf)]
            q = substitute(p, t+tk)
            c = k - numintervals(tbf) + 1
            while c > 0
                for d = 0 : degree(q)
                    r = d+1
                    w = q[d]
                    j = (numfuncs[c,r] += 1)
                    data[j,r,c] = (k,w)
                end
                c = c-1
            end
        end
    end

    return AssemblyData(data)
end

function assemblydata(tbf::TimeBasisDelta)

    T = scalartype(tbf)
    Δt = timestep(tbf)

    z = zero(scalartype(tbf))
    w = one(scalartype(tbf))

    num_cells = numfunctions(tbf)
    num_refs  = 1

    max_num_funcs = 1
    num_funcs = zeros(Int, num_cells, num_refs)
    data = fill((0,z), max_num_funcs, num_refs, num_cells)

    els = [ simplex(point((i-0)*Δt),point((i+1)*Δt)) for i in 1:num_cells ]

    for k in 1 : numfunctions(tbf)-1
        data[1,1,k] = (k+1,w)
    end

    return els, AssemblyData(data)
end

function convolve(f::TimeBasisFunction, g::TimeBasisFunction)
    # @info "BEAST.convolve: computing quadratic spline out of a pulse and a hat"
    dt = timestep(f)
    fg = timebasisspline2(dt, numfunctions(f), scalartype(f))
    fg.polys = dt * fg.polys
    return fg
end

# function convolve(f::TimeBasisFunction{T,2,1,0}, g::TimeBasisFunction{T,2,1,0}) where {T}
#     @info "Convolving two pulses into a hat"
#     dt = timestep(f)
#     fg = timebasisshiftedlagrange(dt, numfunctions(f), 1, scalartype(f))
#     fg.polys = dt * fg.polys
#     return fg
# end

function convolve(δ::TimeBasisDelta, g::TimeBasisFunction)
    # @info "Computing the trivial convolution between a delta and something else"
    return g
end

function convolve(f::TimeBasisFunction{T,2,1,0}, g::TimeBasisFunction{T,2,1,0}) where {T}
    # @info "Convolving two pulses into a hat"
    dt = BEAST.timestep(f)
    fg = timebasisshiftedlagrange(dt, numfunctions(f), 1, scalartype(f))
    fg.polys = dt * fg.polys
    return fg
end

function BEAST.convolve(
    f::TimeBasisFunction{T,2,1,0},
    g::TimeBasisFunction{T,N,D1,D}) where {T,N,D1,D}

    # @info "Convolving a left-of-origin pulse with a piecewise poly"

    @assert degree(f.polys[1]) == 0
    @assert timestep(f) ≈ timestep(g)

    # The result of convolving with a LoO pulse extends the support
    # by Δt to the right.
    Δt = timestep(f)
    P = Polynomial{D1+1,T}
    polys = Vector{P}(undef, N+1)
    t = Polynomial{2,T}(SVector{2,T}(0, 1))

    # Treat the leftmost interval separately
    polys[1] = integrate(g.polys[1], -Δt, 0)

    for i in 2:N
        q1 = -1 * substitute(integrate(g.polys[i-1], (i-2)*Δt, 0), t-Δt)
        q2 = integrate(g.polys[i], (i-2)*Δt, 0)
        polys[i] = q1 + q2
    end

    q1 = -1 * substitute(integrate(g.polys[N], (N-1)*Δt, 0), t-Δt)
    q2 = integrate(g.polys[N], (N-1)*Δt, 0)
    polys[N+1] = q1 + q2

    TimeBasisFunction{T,N+1,D1+1,D+1}(
        g.timestep,
        g.numfunctions,
        polys)
end

"""
    TimeBasisDeltaShifted{T}

Represents a TimeBasisDelta{T} retarded by a fraction of the time step.
"""
struct TimeBasisDeltaShifted{T} <: AbstractTimeBasisFunction
	tbf   :: TimeBasisDelta{T}
	shift :: T
end
scalartype(x::TimeBasisDeltaShifted) = scalartype(x.tbf)
numfunctions(x::TimeBasisDeltaShifted) = numfunctions(x.tbf)
refspace(x::TimeBasisDeltaShifted) = refspace(x.tbf)
timestep(x::TimeBasisDeltaShifted) = timestep(x.tbf)
numintervals(x::TimeBasisDeltaShifted) = numintervals(x.tbf)

function assemblydata(tbds::TimeBasisDeltaShifted)
	tbf = tbds.tbf

    T = scalartype(tbf)
    Δt = timestep(tbf)

    z = zero(scalartype(tbf))
    w = one(scalartype(tbf))

    num_cells = numfunctions(tbf)
    num_refs  = 1

    max_num_funcs = 1
    num_funcs = zeros(Int, num_cells, num_refs)
    data = fill((0,z), max_num_funcs, num_refs, num_cells)

    els = [ simplex(point((i-0+tbds.shift)*Δt),point((i+1+tbds.shift)*Δt)) for i in 1:num_cells ]

    for k in 1 : numfunctions(tbf)
        data[1,1,k] = (k,w)
    end

    return els, AssemblyData(data)
end


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the linear advection equation

advection_velocity = (0.2, -0.7)
equations = LinearScalarAdvectionEquation2D(advection_velocity)

# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)

# First mesh is the left half of a [-1,1]^2 square
coordinates_min1 = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
coordinates_max1 = ( 0.0,  1.0) # maximum coordinates (max(x), max(y))

cells_per_dimension = (4, 8)

mesh1 = StructuredMesh(cells_per_dimension, coordinates_min1, coordinates_max1)

# A semidiscretization collects data structures and functions for the spatial discretization
semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition_convergence_test, solver,
                                     boundary_conditions=(
                                       # Connect left boundary with right boundary of right mesh
                                       x_neg=BoundaryConditionCoupled(2, (:end, :i_forward), Float64),
                                       # Connect right boundary with left boundary of right mesh
                                       x_pos=BoundaryConditionCoupled(2, (:begin, :i_forward),  Float64),
                                       y_neg=boundary_condition_periodic,
                                       y_pos=boundary_condition_periodic))

# Second mesh is the right half of a [-1,1]^2 square
coordinates_min2 = (0.0, -1.0) # minimum coordinates (min(x), min(y))
coordinates_max2 = (1.0,  1.0) # maximum coordinates (max(x), max(y))

mesh2 = StructuredMesh(cells_per_dimension, coordinates_min2, coordinates_max2)

semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition_convergence_test, solver,
                                     boundary_conditions=(
                                       # Connect left boundary with right boundary of left mesh
                                       x_neg=BoundaryConditionCoupled(1, (:end, :i_forward), Float64),
                                       # Connect right boundary with left boundary of left mesh
                                       x_pos=BoundaryConditionCoupled(1, (:begin, :i_forward),  Float64),
                                       y_neg=boundary_condition_periodic,
                                       y_pos=boundary_condition_periodic))

# Create a semidiscretization that bundles semi1 and semi2
semi = SemidiscretizationCoupled((semi1, semi2))

###############################################################################
# ODE solvers, callbacks etc.

# Create ODE problem with time span from 0.0 to 1.0
ode = semidiscretize(semi, (0.0, 1.0));

# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
# and resets the timers
summary_callback = SummaryCallback()

# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
analysis_callback = AnalysisCallback(semi, interval=100)

# The SaveSolutionCallback allows to save the solution to a file in regular intervals
save_solution = SaveSolutionCallback(interval=100,
                                     solution_variables=cons2prim)

# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step
stepsize_callback = StepsizeCallback(cfl=1.6)

# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
callbacks = CallbackSet(summary_callback, stepsize_callback, analysis_callback, save_solution)

###############################################################################
# run the simulation

# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);

# Print the timer summary
summary_callback()


text:
using FileIO
using Images
using ImageSegmentation
using Statistics

"""
Use Laplace of Gaussian(LoG) filter extract cell from raw 3d image
Version Comment
0.1		initial
0.2		only LoG
"""

"""
Generate border form watershed result
"""
function watershedborder(watershed_segments)
    marker_border = BitArray(undef, size(watershed_segments.image_indexmap));
    marker_border .= false
    for label in watershed_segments.segment_labels
        marker_border .|= ((watershed_segments.image_indexmap.==label) .⊻ erode(watershed_segments.image_indexmap .==label));
    end
    marker_border;
end

"""
Use LoG fiter raw image to extract cell
"""
function split_cell_LoG(stack::Array{Gray{Normed{UInt16,16}},4};
                        LoG::Integer=40, thres = -1e-7)
	println("Applying LoG(40) at Maximum Z Projection")
    #img_edge = zeros(N0f16, 1900, 1300, time);
    #mask_edge = zeros(Int16, 1900, 1300, time);
    h, w, d, time = size(stack)
    #local mask_markers = zeros(Bool, h, w, time);
    local mask_markers = Array{Bool}(undef, h, w, time);
	GC.gc() # garbage clean imediately to avoid double free insize threads.@threads
    #@inbounds Threads.@threads for t in 1:time  #use 40 threads slow down speed. may due to gc time
    @inbounds for t in 1:time  #use 40 threads slow down speed. may due to gc time
		# remove possion noise with median filter on maximum z-project image
		#local imgx = mapwindow(median!,
        # maximum(view(stack, :, :, :, t), dims=3)[:,:,1], (5,5));
        # LoG will blur image to low image struture under \sigma,
        # So I remove median filter
        local imgx = maximum(view(stack, :, :, :, t), dims=3)[:,:,1];
		# extract intensity info with LoG
        mask_markers[:,:,t] = imfilter(imgx, Kernel.LoG(LoG)) .< thres ;
        #imgx_dist = distance_transform(feature_transform(imgx_log));
		# filter markers for watershed
        #imgx_markers = label_components( imgx_dist .> 50);
		#mask_markers[:,:,i] = imgx_markers
        #imgx_segments = watershed( imfilter(1 .- imgx, Kernel.gaussian(9)), imgx_markers);
        #img_edge[:,:,i] = .~watershedborder(imgx_segments).*imgx;
        #mask_clear[:,:,i] = extract_nucleus( imgx, imgx_segments) .> 0;
		#mask_edge[:,:,i] = imgx_segments.image_indexmap;
        #print("*");
    end
	println("Done")
    mask_markers;
end

#data_dir = "/datahub/rawdata/tandeng/mRNA_imaging/mRNA_confocal_hamamatsu-60X-TIRF";
#img_16_2 = load(File(format"TIFF", "$data_dir/20200316/HE7-11-1-80uw-PWM_1_s2.ome.tiff"));

#@time edge, clear = split_cell_LoG(img_16_2, 137);
#res_dir = "/datahub/rawdata/tandeng/mRNA_imaging/CoutingmRNA.jl"
#save("output/img_16_2_edge_all.tiff", edge);
#save("output/img_16_2_clear_all.tiff", clear);
#h5write("output/img_16_2_clear_all.h5", "img", rawview(clear));


text:
"""
Perform a union between two HyperRectangles.
"""
function union(h1::HyperRectangle{N, T}, h2::HyperRectangle{N, T}) where {T,N}
    m = min.(minimum(h1), minimum(h2))
    mm = max.(maximum(h1), maximum(h2))
    HyperRectangle{N, T}(m, mm-m)
end

"""
Perform a difference between two HyperRectangles.
"""
diff(h1::HyperRectangle, h2::HyperRectangle) = h1

"""
Perform a intersection between two HyperRectangles.
"""
function intersect(h1::HyperRectangle{N, T}, h2::HyperRectangle{N, T}) where {T,N}
    m = max.(minimum(h1), minimum(h2))
    mm = min.(maximum(h1), maximum(h2))
    HyperRectangle{N, T}(m, mm-m)
end

function intersect(a::SimpleRectangle, b::SimpleRectangle)
    min_n = max.(minimum(a), minimum(b))
    max_n = min.(maximum(a), maximum(b))
    w = max_n - min_n
    SimpleRectangle(min_n[1], min_n[2], w[1], w[2])
end


text:
using JuMP
using Gurobi
include("parser.jl")

function FastBendersGurobi(MyFileName::String, nn = 0::Int64, bendersbool = false::Bool)
  n, f, c, w, d = read_instance(MyFileName)
  if nn != 0
    n = nn
    d = n÷2
    f = Array{Int64,1}(zeros(n))
    c = Array{Int64,1}(zeros(n))
    f[1], c[1] = 7, 8
    for i in 2:n
      f[i] = (f[i-1]*f[1])%159
      c[i] = (c[i-1]*c[1])%61
    end
  end
  if bendersbool
    w = 1
  end
  b = 1
  v = Array{Int64,1}(zeros(n-1))

  #### probleme maitre ####
  # Create the model
  m = Model(CPLEX.Optimizer)

  ## Variables
  @variable(m, y[1:n-1], Bin)
  @variable(m, w >= 0)

  ## Constraints
  @constraint(m, w >= d*b - sum(y[i]*v[i] for i in 1:n-1))
  @constraint(m, y[1] >= y[2])
  @constraint(m, y[1] >= y[3])

  ## Objective
  @objective(m, Min, w + sum(f[i]*y[i] for i in 1:n-1))

  #resolution
  optimize!(m)

  y_star = JuMP.getvalue.( m[:y] )
  y_star = Array{Int64,1}(ones(n-1))
  w_star = JuMP.getvalue.( m[:w] )
  b_star = 1000
  v_star = v

  while w_star < d*b_star - sum(y_star[i]*v_star[i] for i in 1:n-1)
    if sum(y_star[i] for i in 1:n-1) > d
      #### sous-probleme ####
      # Create the model
      m1 = Model(CPLEX.Optimizer)

      ## Variables
      @variable(m1, v[1:n-1] >= 0)
      @variable(m1, b)

      ## Constraints
      @constraint(m1, [i in 1:n-1], b - v[i] <= c[i])

      ## Objective
      @objective(m1, Max, d*b - sum(y_star[i]*v[i] for i in 1:n-1))

      #resolution
      optimize!(m1)

      w_star  = JuMP.objective_value.(m1)
      v_star = JuMP.getvalue.( m1[:v] )
      b_star = JuMP.getvalue.( m1[:b] )

      #### adding Constraints
      @constraint(m, w >= d*b_star - sum(y[i]*v_star[i] for i in 1:n-1))

      optimize!(m)
      y_star = JuMP.getvalue.( m[:y] )
      w_star = JuMP.getvalue.( m[:w] )
    end
    if sum(y_star[i] for i in 1:n-1) == d
      # sans PL
      b_star = max(c[i], i in 1:n-1)
      for i in 1:n-1
        v_star[i] = b_star - c[i]
      end
    end
  end
end

max_stars_repo_name: boumpteryx/MLA

max_issues_repo_name: boumpteryx/MLA

max_forks_repo_name: boumpteryx/MLA


text:
# TODO: make a nice package
# some docs

using Nullables
import Statistics: mean, std
import DataStructures: heappop!, heappush!
import Base: isless, push!

const HEAP_GC_FAC = 0.2 # fraction of used cells before gc'ing the heap

include("data_structures.jl")
include("heap.jl")
include("tree.jl")


text:
const ERROR_MSG = "Pretend Error"

struct AsyncUser end

hear(s::Scene{AsyncUser}, ::Val{:die!}) = async(s) do s
    error(ERROR_MSG)
end

struct AsyncPlay end

function hear(s::Scene{AsyncPlay}, ::Genesis!)
    user = invite!(s, AsyncUser())

    async(s) do s
        @test my(s) isa AsyncPlay

        say(s, user, Val(:die!))
    end

    msg = expect(s, Died!)
    @test msg.corpse.task.exception.task.exception.msg == ERROR_MSG
    leave!(s)
end

@testset LuvvyTestSet expect=2 "Async" begin
    play!(AsyncPlay())
end


text:
abstract type Conductance end
abstract type IonChannel <: Conductance end
abstract type Synapse <: Conductance end

"""
return unparameterised type as symbol
"""
function get_name(ch::Conductance)
    Base.typename(ch |> typeof).name |> Symbol
end

"""
fallback option for channels without a calcium current
"""
calcium_current(ch::IonChannel, sys::ODESystem) = Num(0)
voltage_hook(V, ch::IonChannel, sys::ODESystem) = V ~ sys.V
calcium_hook(Ca, ch::IonChannel, sys::ODESystem) = Ca ~ sys.Ca

"""
does the channel have a reversal
"""
get_reversal(c::Conductance) = nothing
#################### Channels ###############################


text:
# ==============================================================================
# The gravity machine (Man of Steel) -> to terraform the world

println("""\nAlgorithme "Gravity machine" --------------------------------\n""")

const verbose = true
const graphic = true

println("-) Active les packages requis\n")
using JuMP, GLPK, PyPlot, Printf, Random
verbose ? println("  Fait \n") : nothing

generateurVisualise = -1

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

include("GMdatastructures.jl") # types, datastructures and global variables specially defined for GM
include("GMparsers.jl")        # parsers of instances and non-dominated points
include("GMgenerators.jl")     # compute the generators giving the L bound set
include("GMjumpModels.jl")     # JuMP models for computing relaxed optima of the SPA
include("GMrounding.jl")       # Startegies for rounding a LP-solution to a 01-solution
include("GMprojection.jl")     # JuMP models for computing the projection on the polytope of the SPA
include("GMmopPrimitives.jl")  # usuals algorithms in multiobjective optimization
include("GMperturbation.jl")   # routines dealing with the perturbation of a solution when a cycle is detected
include("GMquality.jl")        # quality indicator of the bound set U generated

# ==============================================================================
# Ajout d'une solution relachee initiale a un generateur

function ajouterX0!(vg::Vector{tGenerateur}, k::Int64, s::tSolution{Float64})

    vg[k].sRel = deepcopy(s) # met en place le generateur \bar{x}^k
    vg[k].sPrj = deepcopy(s) # le generateur est la premiere projection \bar{x}^{k,0}
    return nothing
end

# ==============================================================================
# Ajout d'une solution entiere (arrondie ou perturbee) a un generateur

function ajouterXtilde!(vg::Vector{tGenerateur}, k::Int64, x::Vector{Int64}, y::Vector{Int64})

    vg[k].sInt.x = copy(x)
    vg[k].sInt.y = copy(y)
    return nothing
end

# ==============================================================================
# Ajout d'une solution fractionnaire (projetee) a un generateur

function ajouterXbar!(vg::Vector{tGenerateur}, k::Int64, x::Vector{Float64}, y::Vector{Float64})

    vg[k].sPrj.x = copy(x)
    vg[k].sPrj.y = copy(y)
    return nothing
end

# ==============================================================================
# Elabore 2 ensembles d'indices selon que xTilde[i] vaut 0 ou 1

function split01(xTilde::Array{Int,1})

   indices0 = (Int64)[]
   indices1 = (Int64)[]

   for i=1:length(xTilde)
       if xTilde[i] == 0
           push!(indices0,i)
       else
           push!(indices1,i)
       end
    end

   return indices0, indices1
end

# ==============================================================================
# test si une solution est admissible en verifiant si sa relaxation lineaire
# conduit a une solution entiere

function estAdmissible(x::Vector{Float64})

    admissible = true
    i=1
    while admissible && i<=length(x)
        if round(x[i], digits=3)!=0.0 && round(x[i], digits=3)!=1.0
            admissible = false
        end
        i+=1
    end
    return admissible
end

# ==============================================================================
# calcule la performance z d'une solution x sur les 2 objectifs

function evaluerSolution(x::Vector{Float64}, c1::Array{Int,1}, c2::Array{Int,1})

    z1 = 0.0; z2 = 0.0
    for i in 1:length(x)
        z1 += x[i] * c1[i]
        z2 += x[i] * c2[i]
    end
    return round(z1, digits=2), round(z2, digits=2)
end

# ==============================================================================
# nettoyage des valeurs des variables d'une solution x relachee sur [0,1]

function nettoyageSolution!(x::Vector{Float64})

    nbvar=length(x)
    for i in 1:nbvar
        if     round(x[i], digits=3) == 0.0
                   x[i] = 0.0
        elseif round(x[i], digits=3) == 1.0
                   x[i] = 1.0
        else
                   x[i] = round(x[i], digits=3)
        end
    end
end

# ==============================================================================
# predicat : verifie si une solution entiere est realisable
function isFeasible(vg::Vector{tGenerateur}, k::Int64)
    #verbose && vg[k].sFea == true ? println("   feasible") : nothing
    return (vg[k].sFea == true)
end

# ==============================================================================
# predicat : verifie si le nombre d'essai maximum a ete tente
function isFinished(trial::Int64, maxTrial::Int64)
#    verbose && trial > maxTrial ? println("   maxTrial") : nothing
    return (trial > maxTrial)
end

# ==============================================================================
# predicat : verifie si le budget de calcul maximum a ete consomme
function isTimeout(temps, maxTime)
#    verbose && time()- temps > maxTime ? println("   maxTime") : nothing
    return (time()- temps > maxTime)
end

# ==============================================================================
# elabore pC le pointeur du cone ouvert vers L

function elaborePointConeOuvertversL(vg::Vector{tGenerateur}, k::Int64, pB::tPoint, pA::tPoint)

    # recupere les coordonnees du point projete
    pC=tPoint(vg[k].sPrj.y[1], vg[k].sPrj.y[2])

    # etablit le point nadir pN au depart des points pA et pB adjacents au generateur k
    pN = tPoint( pA.x , pB.y )

#    print("Coordonnees du cone 2 : ")
#    @show pC, pN

    # retient pN si pC domine pN (afin d'ouvrir le cone)
    if pC.x < pN.x  &&  pC.y < pN.y
        # remplace pC par pN
        pC=tPoint( pA.x , pB.y )
    end

    return pC
end

# ==============================================================================
#= Retourne un booléen indiquant si un point se trouve dans un secteur défini dans
  le sens de rotation trigonométrique (repère X de gauche à droite, Y du haut vers
  le bas).
  https://www.stashofcode.fr/presence-dun-point-dans-un-secteur-angulaire/#more-328
  M    Point dont la position est à tester (point resultant a tester)
  O    Point sommet du secteur (point generateur)
  A    Point de départ du secteur (point adjacent inferieur)
  B    Point d'arrivée du secteur (point adjacent superieur)
  sortie : Booléen indiquant si le point est dans le secteur ou non.

  Exemple :

  B=point(2.0,1.0)
  O=point(2.5,2.5)
  A=point(5.0,5.0)

  M=point(5.0,4.0)
  inSector(M, O, A, B)
=#

function inSector(M, O, A, B)

    cpAB = (A.y - O.y) * (B.x - O.x) - (A.x - O.x) * (B.y - O.y)
    cpAM = (A.y - O.y) * (M.x - O.x) - (A.x - O.x) * (M.y - O.y)
    cpBM = (B.y - O.y) * (M.x - O.x) - (B.x - O.x) * (M.y - O.y)

    if (cpAB > 0)
        if ((cpAM > 0) && (cpBM < 0))
            return true
        else
            return false
        end
    else
        if (!((cpAM < 0) && (cpBM > 0)))
            return true
        else
            return false
        end
    end
end

function inCone(pOrg, pDeb, pFin, pCur)
    # pOrg : point origine du cone (la ou il est pointe)
    # pDeb : point depart du cone (point du rayon [pOrg,pDeb])
    # pFin : point final du cone (point du rayon [pOrg,pFin])
    # pCur : point courant a tester
    # retourne VRAI si pCur est dans le cone pDeb-pFin-pOrg, FAUX sinon

    cp_pDeb_pFin = (pDeb.x - pOrg.x) * (pFin.y - pOrg.y) - (pDeb.y - pOrg.y) * (pFin.x - pOrg.x)
    cp_pDeb_pCur = (pDeb.x - pOrg.x) * (pCur.y - pOrg.y) - (pDeb.y - pOrg.y) * (pCur.x - pOrg.x)
    cp_pFin_pCur = (pFin.x - pOrg.x) * (pCur.y - pOrg.y) - (pFin.y - pOrg.y) * (pCur.x - pOrg.x)

    if (cp_pDeb_pFin > 0)
        if ((cp_pDeb_pCur >= 0) && (cp_pFin_pCur <= 0))
            return true
        else
            return false
        end
    else
        if (!((cp_pDeb_pCur < 0) && (cp_pFin_pCur > 0)))
            return true
        else
            return false
        end
    end
end

function inCone1VersZ(pOrg, pDeb, pFin, pCur)
    return inCone(pOrg, pDeb, pFin, pCur)
end

function inCone2Vers0(pOrg, pDeb, pFin, pCur)
    return !inCone(pOrg, pDeb, pFin, pCur)
end

# ==============================================================================
# Selectionne les points pour le cone pointe sur le generateur k (pCour) et ouvert vers Y
function selectionPoints(vg::Vector{tGenerateur}, k::Int64)
    nbgen = size(vg,1)
    if k==1
        # premier generateur (point predecesseur fictif)
        pPrec = tPoint(vg[k].sRel.y[1], vg[k].sRel.y[2]+1.0)
        pCour = tPoint(vg[k].sRel.y[1], vg[k].sRel.y[2])
        pSuiv = tPoint(vg[k+1].sRel.y[1], vg[k+1].sRel.y[2])
    elseif k==nbgen
        # dernier generateur (point suivant fictif)
        pPrec = tPoint(vg[k-1].sRel.y[1], vg[k-1].sRel.y[2])
        pCour = tPoint(vg[k].sRel.y[1], vg[k].sRel.y[2])
        pSuiv = tPoint(vg[k].sRel.y[1]+1.0, vg[k].sRel.y[2])
    else
        # generateur non extreme
        pPrec = tPoint(vg[k-1].sRel.y[1], vg[k-1].sRel.y[2])
        pCour = tPoint(vg[k].sRel.y[1], vg[k].sRel.y[2])
        pSuiv = tPoint(vg[k+1].sRel.y[1], vg[k+1].sRel.y[2])
    end
#    print("Coordonnees du cone 1 : ")
#    @show pPrec, pCour, pSuiv
    return pPrec, pCour, pSuiv
end

# ==============================================================================
# Calcule la direction d'interet du nadir vers le milieu de segment reliant deux points generateurs
function calculerDirections(L::Vector{tSolution{Float64}}, vg::Vector{tGenerateur})
    # function calculerDirections(L, vg::Vector{tGenerateur})

     nbgen = size(vg,1)
     λ1=Vector{Float64}(undef, nbgen - 1)
     λ2=Vector{Float64}(undef, nbgen - 1)
     for k in 2:nbgen

         n1 = L[end].y[1]
         n2 = L[1].y[2]

         x1,y1 = vg[k-1].sRel.y[1], vg[k-1].sRel.y[2]
         x2,y2 = vg[k].sRel.y[1], vg[k].sRel.y[2]
         xm=(x1+x2)/2.0
         ym=(y1+y2)/2.0
         Δx = abs(n1-xm)
         Δy = abs(n2-ym)
         λ1[k-1] =  1 - Δx / (Δx+Δy)
         λ2[k-1] =  1 - Δy / (Δx+Δy)
         @printf("  k= %3d   \n",k)
         @printf("  x1= %7.2f   y1= %7.2f \n",x1,y1)
         @printf("  x2= %7.2f   y2= %7.2f \n",x2,y2)
         @printf("  Δx= %7.2f    Δy= %7.2f \n",Δx,Δy)
         @printf("  λ1= %6.5f    λ2= %6.5f \n",λ1[k-1],λ2[k-1])
         plot(n1, n2, xm, ym, linestyle="-", color="blue", marker="+")
         annotate("",
                  xy=[xm;ym],# Arrow tip
                  xytext=[n1;n2], # Text offset from tip
                  arrowprops=Dict("arrowstyle"=>"->"))
         println("")
     end
     return λ1, λ2
 end

# ==============================================================================
# Calcule la direction d'interet du nadir vers un point generateur
function calculerDirections2(L::Vector{tSolution{Float64}}, vg::Vector{tGenerateur})
    #function calculerDirections2(L, vg::Vector{tGenerateur})

    nbgen = size(vg,1)
    λ1=Vector{Float64}(undef, nbgen)
    λ2=Vector{Float64}(undef, nbgen)
    for k in 1:nbgen

        n1 = L[end].y[1]
        n2 = L[1].y[2]

        xm=vg[k].sRel.y[1]
        ym=vg[k].sRel.y[2]
        Δx = abs(n1-xm)
        Δy = abs(n2-ym)
        λ1[k] =  1 - Δx / (Δx+Δy)
        λ2[k] =  1 - Δy / (Δx+Δy)
        @printf("  k= %3d   ",k)
        @printf("  xm= %7.2f   ym= %7.2f ",xm,ym)
        @printf("  Δx= %8.2f    Δy= %8.2f ",Δx,Δy)
        @printf("  λ1= %6.5f    λ2= %6.5f \n",λ1[k],λ2[k])
        if generateurVisualise == -1
            # affichage pour tous les generateurs
            plot(n1, n2, xm, ym, linestyle="-", color="blue", marker="+")
            annotate("",
                     xy=[xm;ym],# Arrow tip
                     xytext=[n1;n2], # Text offset from tip
                     arrowprops=Dict("arrowstyle"=>"->"))
        elseif generateurVisualise == k
            # affichage seulement pour le generateur k
            plot(n1, n2, xm, ym, linestyle="-", color="blue", marker="+")
            annotate("",
                     xy=[xm;ym],# Arrow tip
                     xytext=[n1;n2], # Text offset from tip
                     arrowprops=Dict("arrowstyle"=>"->"))
        end
        #println("")
    end
    return λ1, λ2
end

# ==============================================================================
# point d'entree principal

function GM( fname::String,
             tailleSampling::Int64,
             maxTrial::Int64,
             maxTime::Int64
           )

    @assert tailleSampling>=3 "Erreur : Au moins 3 sont requis"

    @printf("0) instance et parametres \n\n")
    verbose ? println("  instance = $fname | tailleSampling = $tailleSampling | maxTrial = $maxTrial | maxTime = $maxTime\n\n") : nothing

    # chargement de l'instance numerique ---------------------------------------
    c1, c2, A = loadInstance2SPA(fname) # instance numerique de SPA
    nbctr = size(A,1)
    nbvar = size(A,2)
    nbobj = 2

    # structure pour les points qui apparaitront dans l'affichage graphique
    d = tListDisplay([],[], [],[], [],[], [],[], [],[], [],[], [],[])

    # --------------------------------------------------------------------------
    # --------------------------------------------------------------------------
    @printf("1) calcule les etendues de valeurs sur les 2 objectifs\n\n")

    # calcule la valeur optimale relachee de f1 seule et le point (z1,z2) correspondant
    f1RL, xf1RL = computeLinearRelax2SPA(nbvar, nbctr, A, c1, c2, typemax(Int), 1) # opt fct 1
    minf1RL, maxf2RL = evaluerSolution(xf1RL, c1, c2)

    # calcule la valeur optimale relachee de f2 seule et le point (z1,z2) correspondant
    f2RL, xf2RL = computeLinearRelax2SPA(nbvar, nbctr, A, c1, c2, typemax(Int), 2) # opt fct 2
    maxf1RL, minf2RL = evaluerSolution(xf2RL, c1, c2)

    verbose ? @printf("  f1_min=%8.2f ↔ f1_max=%8.2f (Δ=%.2f) \n",minf1RL, maxf1RL, maxf1RL-minf1RL) : nothing
    verbose ? @printf("  f2_min=%8.2f ↔ f2_max=%8.2f (Δ=%.2f) \n\n",minf2RL, maxf2RL, maxf2RL-minf2RL) : nothing

    # --------------------------------------------------------------------------
    # --------------------------------------------------------------------------
    @printf("2) calcule les generateurs par e-contrainte alternant minimiser z1 et z2\n\n")

    nbgen, L = calculGenerateurs(A, c1, c2, tailleSampling, minf1RL, maxf2RL, maxf1RL, minf2RL, d)

    # --------------------------------------------------------------------------
    # --------------------------------------------------------------------------
    # allocation de memoire pour la structure de donnees -----------------------

    vg = allocateDatastructure(nbgen, nbvar, nbobj)

    # --------------------------------------------------------------------------
    # --------------------------------------------------------------------------
    @printf("3) place L dans structure et verifie l'admissibilite de chaque generateur\n\n")

    for k=1:nbgen

        verbose ? @printf("  %2d  : [ %8.2f , %8.2f ] ", k, L[k].y[1], L[k].y[2]) : nothing

        # copie de l'ensemble bornant inferieur dans la stru de donnees iterative ---
        ajouterX0!(vg, k, L[k])

        # test d'admissibilite et marquage de la solution le cas echeant -------
        if estAdmissible(vg[k].sRel.x)
            ajouterXtilde!(vg, k, convert.(Int, vg[k].sRel.x), convert.(Int, L[k].y))
            vg[k].sFea   = true
            verbose ? @printf("→ Admissible \n") : nothing
            # archive le point obtenu pour les besoins d'affichage
            if generateurVisualise == -1
                # archivage pour tous les generateurs
                push!(d.XFeas,vg[k].sInt.y[1])
                push!(d.YFeas,vg[k].sInt.y[2])
            elseif generateurVisualise == k
                # archivage seulement pour le generateur k
                push!(d.XFeas,vg[k].sInt.y[1])
                push!(d.YFeas,vg[k].sInt.y[2])
            end
        else
            vg[k].sFea   = false
            verbose ? @printf("→ x          \n") : nothing
        end

    end
    verbose ? println("") : nothing

    # --------------------------------------------------------------------------
    # --------------------------------------------------------------------------
    # Sortie graphique

    figure("Gravity Machine",figsize=(6.5,5))
    #xlim(25000,45000)
    #ylim(20000,40000)
    xlabel(L"z^1(x)")
    ylabel(L"z^2(x)")
    PyPlot.title("Cone | 1 rounding | 2-$fname")

    # --------------------------------------------------------------------------
    # --------------------------------------------------------------------------
    # calcule les directions (λ1,λ2) pour chaque generateur a utiliser lors des projections
    λ1,λ2 = calculerDirections2(L,vg)

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

    @printf("4) terraformation generateur par generateur \n\n")

    nbIterTotal = 0
    nbgenNotFeasible = 0
    for k in [i for i in 1:nbgen if !isFeasible(vg,i)]
        temps = time()
        trial = 0
        H =(Vector{Int64})[]

#perturbSolution30!(vg,k,c1,c2,d)

        # rounding solution : met a jour sInt dans vg --------------------------
        #roundingSolution!(vg,k,c1,c2,d)  # un cone
        #roundingSolutionnew24!(vg,k,c1,c2,d) # deux cones
        roundingSolutionNew23!(vg,k,c1,c2,d) # un cone et LS sur generateur

        push!(H,[vg[k].sInt.y[1],vg[k].sInt.y[2]])
        println("   t=",trial,"  |  Tps=", round(time()- temps, digits=4))

        while !(t1=isFeasible(vg,k)) && !(t2=isFinished(trial, maxTrial)) && !(t3=isTimeout(temps, maxTime))

            trial+=1

            # projecting solution : met a jour sPrj, sInt, sFea dans vg --------
            projectingSolution!(vg,k,A,c1,c2,λ1,λ2,d)
            println("   t=",trial,"  |  Tps=", round(time()- temps, digits=4))

            if !isFeasible(vg,k)

                # rounding solution : met a jour sInt dans vg --------------------------
                #roundingSolution!(vg,k,c1,c2,d)
                #roundingSolutionnew24!(vg,k,c1,c2,d)
                roundingSolutionNew23!(vg,k,c1,c2,d)
                println("   t=",trial,"  |  Tps=", round(time()- temps, digits=4))

                # test detection cycle sur solutions entieres ------------------
                cycle = [vg[k].sInt.y[1],vg[k].sInt.y[2]] in H
                if (cycle == true)
                    println("CYCLE!!!!!!!!!!!!!!!")
                    # perturb solution
                    perturbSolution30!(vg,k,c1,c2,d)
                end
                push!(H,[vg[k].sInt.y[1],vg[k].sInt.y[2]])

            end
        end
        if t1
            println("   feasible \n")
        elseif t2
            println("   maxTrial \n")
        elseif t3
            println("   maxTime \n")
        end
        nbIterTotal += trial
        nbgenNotFeasible += 1
    end
    verbose ? @printf("   Nombre d'itération moyenne afin de trouver une solution admissible : %5.3f", nbIterTotal/nbgenNotFeasible) : nothing
    println("");

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

    @printf("5) Extraction des resultats\n\n")

    for k=1:nbgen
        verbose ? @printf("  %2d  : [ %8.2f , %8.2f ] ", k, vg[k].sInt.y[1],vg[k].sInt.y[2]) : nothing
        # test d'admissibilite et marquage de la solution le cas echeant -------
        if vg[k].sFea
            verbose ? @printf("→ Admissible \n") : nothing
        else
            verbose ? @printf("→ x          \n") : nothing
        end
    end

    # allocation de memoire pour les ensembles bornants ------------------------
    U = Vector{tSolution{Int64}}(undef,nbgen)
    for j = 1:nbgen
        U[j] = tSolution{Int64}(zeros(Int64,nbvar),zeros(Int64,nbobj))
    end
    #--> TODO : stocker l'EBP dans U proprement

    # ==========================================================================
    @printf("6) Edition des resultats \n\n")

#    figure("Gravity Machine",figsize=(6.5,5))
    #xlim(25000,45000)
    #ylim(20000,40000)
#    xlabel(L"z^1(x)")
#    ylabel(L"z^2(x)")
    # Donne les points relaches initiaux ---------------------------------------
#    scatter(d.xLf1,d.yLf1,color="blue", marker="x")
#    scatter(d.xLf2,d.yLf2,color="red", marker="+")
    graphic ? scatter(d.xL,d.yL,color="blue", marker="x", label = L"y \in L") : nothing

    # Donne les points entiers -------------------------------------------------
    graphic ? scatter(d.XInt,d.YInt,color="orange", marker="s", label = L"y"*" rounded") : nothing
#    @show d.XInt
#    @show d.YInt

    # Donne les points apres projection Δ(x,x̃) ---------------------------------
    graphic ? scatter(d.XProj,d.YProj, color="red", marker="x", label = L"y"*" projected") : nothing
#    @show d.XProj
#    @show d.YProj

    # Donne les points admissibles ---------------------------------------------
    graphic ? scatter(d.XFeas,d.YFeas, color="green", marker="o", label = L"y \in F") : nothing
#    @show d.XFeas
#    @show d.YFeas

    # Donne l'ensemble bornant primal obtenu + la frontiere correspondante -----
    #--> TODO : stocker l'EBP dans U proprement
    X_EBP_frontiere, Y_EBP_frontiere, X_EBP, Y_EBP = ExtractEBP(d.XFeas, d.YFeas)
    plot(X_EBP_frontiere, Y_EBP_frontiere, color="green", markersize=3.0, marker="x")
    scatter(X_EBP, Y_EBP, color="green", s = 150, alpha = 0.3, label = L"y \in U")
    @show X_EBP
    @show Y_EBP

    # Donne les points qui ont fait l'objet d'une perturbation -----------------
     scatter(d.XPert,d.YPert, color="magenta", marker="s", label ="pertub")

    # Donne les points non-domines exacts de cette instance --------------------
     XN,YN = loadNDPoints2SPA(fname)
     plot(XN, YN, color="black", linewidth=0.75, marker="+", markersize=1.0, linestyle=":", label = L"y \in Y_N")
     scatter(XN, YN, color="black", marker="+")
    @show XN
    @show YN

    # Affiche le cadre avec les legendes des differents traces -----------------
    legend(bbox_to_anchor=[1,1], loc=0, borderaxespad=0, fontsize = "x-small")
    #PyPlot.title("Cone | 1 rounding | 2-$fname")

    # Compute the quality indicator of the bound set U generated ---------------
    # Need at least 2 points in EBP to compute the quality indicator
    if length(X_EBP) > 1
        quality = qualityMeasure(XN,YN, X_EBP,Y_EBP)
        @printf("Quality measure: %5.2f %%\n", quality*100)
    end

end

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

#@time GM("sppaa02.txt", 6, 20, 20)
#@time GM("sppnw03.txt", 6, 20, 20) #pb glpk
#@time GM("sppnw10.txt", 6, 20, 20)
@time GM("didactic5.txt", 5, 5, 10)
#@time GM("sppnw29.txt", 6, 30, 20)
nothing

max_stars_repo_name: ducanh-le/GravityMachine

max_issues_repo_name: ducanh-le/GravityMachine

max_forks_repo_name: ducanh-le/GravityMachine


text:
#=
main:
- Julia version: 1.3.1
- Author: abd_alazez_ahmed
- Date: 2020-04-22
=#
using StatsBase
#--2-a--#
println("\n\nQusetion 2 a\n");
n=sample(1:100, 5, replace = false);
# n= [64 93 70 33 55];
t_max=81;
t_min=31;
println(string(length(n)," random numbers between 0 and 100 are : ",n));
t=round.(Int,t_min.+n.*(t_max-t_min)./100);
println(string("The corresponding time is : ",t));
#--2-a--#
println("\n\nQusetion 2 b\n");
n=sample(1:100, 5, replace = false);
# n=[45 93 96 1 73];
t_mean=12;
println(string(length(n)," random numbers between 0 and 100 are : ",n));
t=round.(Int,t_mean.*log.(100 ./(100 .-n)));
println(string("The corresponding time is : ",t));


text:
@testset "$(rpad("Splitting Tableaus",80))" begin

    # instatiate all splitting tableaus
    @test typeof(getTableauLieA()) <: TableauSplittingNS{Float64}
    @test typeof(getTableauLieB()) <: TableauSplittingNS{Float64}
    @test typeof(getTableauStrang()) <: TableauSplittingNS{Float64}
    @test typeof(getTableauMcLachlan2()) <: TableauSplittingNS{Float64}
    @test typeof(getTableauMcLachlan4()) <: TableauSplittingNS{Float64}

    @test typeof(getTableauTripleJump()) <: TableauSplittingSS{Float64}
    @test typeof(getTableauSuzukiFractal()) <: TableauSplittingSS{Float64}

end


text:
include("../utils/format/Heading.jl")

Heading.main("Default type")
Heading.sub("type for an integer literal depends on whether the target system is 32/64 bit")

println(Sys.WORD_SIZE)
#: 64

println(typeof(1))
#: Int64

# explicitly give make 32
println(typeof(Int32(1)))
#: Int32

# unsigned UInt is given by hex.
# size of UInt increases with num of hex
println(typeof(UInt(1)))
println(typeof(0x1)) # small => UInt8
println(typeof(0x0123456789abcdef)) # big => UInt64
#: UInt64
#: UInt8
#: UInt64

Heading.main("All types with ranges")
rpadval = 50
for T in [Int8, Int16, Int32, Int64, Int128, UInt8, UInt16, UInt32, UInt64, UInt128]
    println("$(T)\t:\t\t $(rpad(typemin(T), rpadval)) $(typemax(T))")
end
#: Int8    :                -128                                               127
#: Int16   :                -32768                                             32767
#: Int32   :                -2147483648                                        2147483647
#: Int64   :                -9223372036854775808                               9223372036854775807
#: Int128  :                -170141183460469231731687303715884105728           170141183460469231731687303715884105727
#: UInt8   :                0                                                  255
#: UInt16  :                0                                                  65535
#: UInt32  :                0                                                  4294967295
#: UInt64  :                0                                                  18446744073709551615
#: UInt128 :                0                                                  340282366920938463463374607431768211455

Heading.main("Modular arithematic overflow")
max_int8 = typemax(Int8)
a = Int8(max_int8) + Int8(1)
print(a == typemin(Int8))
#: true

# Note: Int8 + Int64 => Int64
a = Int8(max_int8) + 1 # type of `1` is `Sys.WORD_SIZE` i.e 64
print(typeof(a))
#: Int64

Heading.main("Solution to Overflow")
Heading.sub("In logics where overflow is possible, do not use Int!!!! Use BigInt")
# BigInt is arbitary precision used as `big`

println(typemax(Int64))
#: 9223372036854775807

println(10^19)
#! -8446744073709551616 (wrong)

println(big(10)^20)
#: 100000000000000000000


text:
# ============================================================
# Geometrical Methods
# ============================================================

export global_frame,
       local_frame,
       unit_normal,
       point_distance
export PSpace1D,
       PSpace2D,
       CSpace2D
export uniform_mesh,
       ndgrid,
       meshgrid,
       find_idx
export UnstructPSpace

include("geo_general.jl")
include("geo_struct.jl")
include("geo_unstruct.jl")


text:
import TPT

using Base.Test

# using Plots; pyplot()

println("--- Unit NFE ---")

#
# Single component: Liquid Na
# Ref: I. H. Umar et al.: J. Phys. F: Metal Phys., Vol. 4 (1974), 1691-1706.
#
ahs1 = TPT.AHS(σ = 6.087, ρ = 3.564e-3)
pp1 = TPT.Ashcroft(1.0, 1.67)
nfe1 = TPT.NFE(ahs1, pp1)
sys1 = TPT.TPTSystem(ahs1, nfe1, T = 373, m = 4.191e+4)

kF1 = TPT.fermiwavenumber(nfe1)
rs = TPT.electrondistance(nfe1)

G1_IU = TPT.localfield_IU(nfe1)
G1_VS = TPT.localfield_VS(nfe1)
# plot([q -> G1_IU(q*kF1), q -> G1_VS(q*kF1)], 0, 10, ylims=(0, 1.5))

ϵ1 = TPT.dielectric(nfe1)
# plot(q -> ϵ(q*kF), 0, 6, ylims=(-1, 10))

ω1₀ = TPT.formfactor(nfe1)[1]
# plot(q -> ω1₀(q*kF), 0, 6, ylims=(-1,1))

F1 = TPT.wnechar(nfe1)[1,1]
# plot(q -> 27.21*F(q*kF), 0, 6, ylims = (-0.1, 0))

u1 = TPT.pairpotential(nfe1)[1,1]
u1′ = TPT.pairpotential_derivative(nfe1)[1,1]
# plot([u1, u1′], 2, 20, ylims = (-2.1e-3, 2.1e-3))

g = TPT.paircorrelation(sys1)[1,1]
# plot(g, 1, 20)

T = 373.
kB = TPT.kB
U1_eg = TPT.internal_eg(nfe1)
S1_eg = TPT.entropy(nfe1, ahs1, T)
U1_es = TPT.internal_es(nfe1, ahs1)
U1_bs = TPT.internal(ahs1, nfe1)
K1 = TPT.kinetic(sys1)
S1_hs = TPT.entropy(sys1) - S1_eg
F1_hs = -T*kB*S1_hs

F1_total = TPT.helmholtz(sys1)

#
# Binary
#
σ_Na = 6.235
σ_K = 7.354

ahs2 = TPT.AHS(σ = [σ_Na, σ_K], ρ = 2.461e-3, c = [0.5, 0.5])
pp2 = TPT.Ashcroft([1.0, 1.0], [1.67, 2.12])
nfe2 = TPT.NFE(ahs2, pp2)
sys2 = TPT.TPTSystem(ahs2, nfe2, T = 373, m = [4.191e+4, 4.191e+4])

ω2₀ = TPT.formfactor(nfe2)
# plot(q -> ω2₀(q*kF), 0, 6, ylims=(-1,1))

F2 = TPT.wnechar(nfe2)
# plot(q -> 27.21*F(q*kF), 0, 6, ylims = (-0.1, 0))

u2 = TPT.pairpotential(nfe2)
# plot(u2, 2, 20, ylims = (-2.1e-3, 2.1e-3))

F2_total = TPT.helmholtz(sys2)

@testset "Unit NFE" begin
  @test isapprox(kF1, 0.473, atol=1e-3)
  @test isapprox(G1_VS(kF1), 0.267, atol=1e-3)
  @test isapprox(ϵ1(kF1), 3.46, atol=1e-2)
  @test isapprox(ω1₀(kF1), -0.141, atol=1e-3)
  @test isapprox(F1(kF1), -0.0437, atol=1e-4)
  @test isapprox(u1(6.087), -0.00078, atol=1e-5)

  @test isapprox(U1_eg, -0.0816, atol=1e-4)
  @test isapprox(S1_eg, 0.0522, atol=1e-4)
  @test isapprox(U1_es, -0.147, atol=1e-3)
  @test isapprox(U1_bs, -0.00833, atol=1e-5)
  @test isapprox(K1, 0.00177, atol=1e-5)
  @test isapprox(F1_hs, -0.00921, atol=1e-5)
end


text:
using OrdinaryDiffEq, ModelingToolkit, Test
@parameters t
@variables y[1:3](t)
@parameters k[1:3]
@derivatives D'~t

eqs = [D(y[1]) ~ -k[1]*y[1] + k[3]*y[2]*y[3],
       D(y[2]) ~  k[1]*y[1] - k[3]*y[2]*y[3] - k[2]*y[2]^2,
       0 ~  y[1] + y[2] + y[3] - 1]

sys = ODESystem(eqs,t,y,k)
M = calculate_massmatrix(sys)
@test M == [1 0 0
            0 1 0
            0 0 0]

f = ODEFunction(sys)
prob_mm = ODEProblem(f,[1.0,0.0,0.0],(0.0,1e5),(0.04,3e7,1e4))
sol = solve(prob_mm,Rodas5(),reltol=1e-8,abstol=1e-8)

function rober(du,u,p,t)
  y₁,y₂,y₃ = u
  k₁,k₂,k₃ = p
  du[1] = -k₁*y₁ + k₃*y₂*y₃
  du[2] =  k₁*y₁ - k₃*y₂*y₃ - k₂*y₂^2
  du[3] =  y₁ + y₂ + y₃ - 1
  nothing
end
f = ODEFunction(rober,mass_matrix=M)
prob_mm2 = ODEProblem(f,[1.0,0.0,0.0],(0.0,1e5),(0.04,3e7,1e4))
sol2 = solve(prob_mm2,Rodas5(),reltol=1e-8,abstol=1e-8)

@test Array(sol) == Array(sol2)


text:
using Conda
using PyCall
import Pkg

env = "py38"

# Change python
if false
    ENV["PYTHON"] = joinpath(Conda.ROOTENV, "envs", env, "python.exe")
    Pkg.build("PyCall")
end

# Default python
if false
    ENV["PYTHON"] = ""
    Pkg.build("PyCall")
end

begin
    # Keyword: VTK, numpy, Julia, import _vtkmodules_static, vtkmodules.all
    # Keyword: Conda, non-root, env
    #
    # From Conda document,
    # NOTE: If you are installing Python packages for use with PyCall,
    # you must use the root environment.
    #
    # The following code allows using env other than root.
    # Add "...\envs\py38\Library\bin" to path
    path = Conda.bin_dir(dirname(PyCall.pyprogramname))
    if isnothing(findfirst(path, ENV["PATH"])) && Sys.iswindows()
        ENV["PATH"] = path * ";" * ENV["PATH"]
    end
    # VTK, numpy etc now will work
end

function learn_basic()
    @show PyCall.pyversion
    @show PyCall.libpython
    @show PyCall.pyprogramname
    @show PyCall.conda

    @show joinpath(dirname(pathof(PyCall)), "deps", "PYTHON")

    sys = pyimport("sys")
    @show sys.version

    math = pyimport("math")
    @show math.sin(math.pi / 4)   # 1/sqrt(2) = 0.70710678

    # use PyCall
    # problem with Qt
    # plt = pyimport("matplotlib.pyplot")
    # x = range(0;stop=2*pi,length=1000); y = sin.(3*x + 4*cos.(2*x));
    # plt.plot(x, y, color="red", linewidth=2.0, linestyle="--") # not working
end

function learn_qt()
    sys = pyimport("sys")
    qws = pyimport("PySide2.QtWidgets")

    app = qws.QApplication([])
    label = qws.QLabel("Qt for Python!")
    label.show()
    #sys.exit(app.exec_())
end

function learn_gtk()
    # conda install pygobject gtk3
    gi = pyimport("gi")
    gi.require_version("Gtk", "3.0")

    grepo = pyimport("gi.repository")
    gtk = pyimport("gi.repository.Gtk")

    win = gtk.Window(title="Hello World")
    win.connect("destroy", gtk.main_quit)
    win.show_all()
    gtk.main()
end

function learn_tkinter()
    tkinter = pyimport(:tkinter)
end

function learn_wx()
    # For graphic heavy library such as wxpython,
    # use py""" """ to import before pyimport(:wx).

    # Must load in Python before loading in Julia
    py"""
    import wx
    """
    wx = pyimport(:wx)
end

#learn_basic()
#learn_qt()
#learn_gtk()

#learn_tkinter()
#learn_wx()

nothing


text:
module Benitor

function bubbleSort!(a)
    n = length(a)

    for  i = 1:n
        for j = 1:n-1
            if a[j] === missing
                a[1] = a[j]
                continue
            end
            if a[j] > a[j+1]
                a[j],a[j+1] = a[j+1],a[j]
            end
        end
    end
return a
end

end # module