text:
module Test_gsfstab

using Random
using LinearAlgebra
using MatrixPencils
using Test

Random.seed!(2351);
function evsym!(ev)
   ev[imag.(ev) .> 0] = conj(ev[imag.(ev) .< 0])
   return ev
end

@testset "Spectrum allocation functions" begin

println("saloc for standard pair (A,B)")
@testset "saloc for standard pair (A,B)" begin

## simple cases: B = [];
a = ones(1,1);  b = zeros(1,0);  evals = [-1];
@time f, SF, blkdims = saloc(a,b, evals =  evals)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f && blkdims == [0, 0, 1]
 evals,eigvals(a+b*f)

## simple cases
a = ones(1,1);  b = zeros(1,5);  evals = [-1];
@time f, SF, blkdims = saloc(a,b, evals =  evals)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f && blkdims == [0, 0, 1]
evals,eigvals(a+b*f)

## simple cases
a = ones(1,1);  b = ones(1,1);  evals = [-1];
@time f, SF, blkdims = saloc(a,b, evals =  evals)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f && blkdims == [0, 1, 0] &&
sort(real( evals)) ≈ sort(real(SF.values)) && sort(imag( evals)) ≈ sort(imag(SF.values))
evals,eigvals(a+b*f)

## simple cases
a = complex(ones(1,1));  b = ones(1,1);  evals = [-im];
@time f, SF, blkdims = saloc(a,b, evals =  evals)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f && blkdims == [0, 1, 0] &&
      sort(real( evals)) ≈ sort(real(SF.values)) && sort(imag( evals)) ≈ sort(imag(SF.values))
evals,eigvals(a+b*f)

## simple cases (no real eigenvalue)
a = ones(1,1);  b = ones(1,1);  evals = [-1+im*0.5;-1-im*0.5];
@time f, SF, blkdims = saloc(a,b, evals =  evals)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f && blkdims == [0, 1, 0]
evals,eigvals(a+b*f)

## simple cases (no real eigenvalue, use sdeg)
a = ones(1,1);  b = ones(1,1);  evals = [-1+im*0.5;-1-im*0.5]; sdeg = -0.5;
@time f, SF, blkdims = saloc(a,b, evals =  evals, sdeg = sdeg)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f && blkdims == [0, 1, 0]
evals,eigvals(a+b*f)

## simple cases
a = ones(1,1);  b = ones(1,1);  evals = [-1;-0.5];
@time f, SF, blkdims = saloc(a,b, evals =  evals)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f && blkdims == [0, 1, 0]
evals,eigvals(a+b*f)

## simple cases
a = [1 0;0 1];  b = rand(2,5);  evals = [-1;-0.5];
@time f, SF, blkdims = saloc(a,b, evals =  evals)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f  && blkdims == [0, 2, 0]
@test sort(real( evals)) ≈ sort(real(SF.values)) && sort(imag( evals)) ≈ sort(imag(SF.values))
evals,eigvals(a+b*f)

 ## simple cases  uncontrollable pair
a = [1 0;0 1];  b = rand(2,1);  evals = [-1;-0.5];
@time f, SF, blkdims = saloc(a,b, evals =  evals)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f && blkdims == [0, 1, 1]
evals,eigvals(a+b*f)

 ## simple cases  controllable pair
 a = [1 0;0 2];  b = rand(2,1);  evals = [-1;-0.5];
 @time f, SF, blkdims = saloc(a,b, evals =  evals)
 @test SF.Z*SF.T*SF.Z' ≈ a+b*f && blkdims == [0, 2, 0] &&
       sort(real( evals)) ≈ sort(real(SF.values)) && sort(imag( evals)) ≈ sort(imag(SF.values))
evals,eigvals(a+b*f)

 ## simple cases  uncontrollable pair
 a = [1 0;0 1];  b = rand(2,1);  evals = [-1+im*0.5;-1-im*0.5];
 @time f, SF, blkdims = saloc(a,b, evals =  evals)
 @test SF.Z*SF.T*SF.Z' ≈ a+b*f  && blkdims == [0, 1, 1]
 evals,eigvals(a+b*f)

## simple cases  controllable pair
a = [1 0;0 2];  b = rand(2,1);  evals = [-1+im*0.5;-1-im*0.5];
@time f, SF, blkdims = saloc(a,b, evals =  evals)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f && blkdims == [0, 2, 0] &&
      sort(real( evals)) ≈ sort(real(SF.values)) && sort(imag( evals)) ≈ sort(imag(SF.values))
evals,eigvals(a+b*f)

## simple cases
a = [1 1;0 2];  b = rand(2,5);  evals = [-1+im*0.5;-1-im*0.5];
@time f, SF, blkdims = saloc(a,b, evals =  evals)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f && blkdims == [0, 2, 0] &&
      sort(real( evals)) ≈ sort(real(SF.values)) && sort(imag( evals)) ≈ sort(imag(SF.values))
evals,eigvals(a+b*f)

## simple cases
a = [1 1 1;-1 1 0;0 0 2];  b = rand(3,5);  evals = [-1+im*0.5;-1-im*0.5;-2];
@time f, SF, blkdims = saloc(a,b, evals =  evals)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f && blkdims == [0, 3, 0] &&
      sort(real( evals)) ≈ sort(real(SF.values)) && sort(imag( evals)) ≈ sort(imag(SF.values))
evals,eigvals(a+b*f)

## simple cases
a = [1 1 0;-1 1 0;0 0 2];  b = [0 0; 0 0; 1 1];  evals = [-1+im*0.5;-1-im*0.5];
@time f, SF, blkdims = saloc(a,b, evals =  evals)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f && blkdims == [0, 1, 2]
evals,eigvals(a+b*f)

## simple cases
a = [1 1; 0 1 ];  b = [0 0; 1 1];  evals = [-1+im*0.5;-1-im*0.5];
@time f, SF, blkdims = saloc(a,b, evals =  evals)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f && blkdims == [0, 2, 0]
evals,eigvals(a+b*f)

## simple cases
a = [1 1 1;-1 1 0;0 0 2];  b = rand(3,5);  evals = [-1;-0.5;-2];
@time f, SF, blkdims = saloc(a,b, evals =  evals)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f && blkdims == [0, 3, 0] &&
      sort(real( evals)) ≈ sort(real(SF.values)) && sort(imag( evals)) ≈ sort(imag(SF.values))
evals,eigvals(a+b*f)

## simple cases
a = [1 1 1;-1 1 0;0 0 2];  b = rand(3,5); sdeg = -0.5;
@time f, SF, blkdims = saloc(a,b,sdeg = sdeg)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f && blkdims == [0, 3, 0]  &&
      all(sdeg .≈ (real(SF.values)))
sdeg,eigvals(a+b*f)

## simple cases
a = [1 1 1;-1 1 0;0 0 2];  b = rand(3,5);  evals = [-1,-0.5]; sdeg = -0.5;
@time f, SF, blkdims = saloc(a,b, evals =  evals, sdeg = sdeg)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f && blkdims == [0, 3, 0]
sdeg, evals, eigvals(a+b*f)

## simple cases
a = [1 1 1;-1 1 0;0 0 2];  b = rand(3,5);  evals = [-3]; sdeg = -2;
@time f, SF, blkdims = saloc(a,b, evals =  evals, sdeg = sdeg)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f && blkdims == [0, 3, 0]
sdeg, evals, eigvals(a+b*f)

## SB01BD EXAMPLE PROGRAM DATA
a = [
  -6.8000   0.0000  -207.0000   0.0000
   1.0000   0.0000     0.0000   0.0000
  43.2000   0.0000     0.0000  -4.2000
   0.0000   0.0000     1.0000   0.0000];
b = [
   5.6400   0.0000
   0.0000   0.0000
   0.0000   1.1800
   0.0000   0.0000];
 evals = [
  -0.5000+im*0.1500
  -0.5000-im*0.1500
  -2.0000
  -0.4000];
sdeg = -0.4;
e = I;
@time f, SF, blkdims = saloc(a,b, evals =  evals, sdeg = sdeg)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f && blkdims == [2, 2, 0]
sdeg, evals, eigvals(a+b*f)

sdeg = missing;
@time f, SF, blkdims = saloc(a,b, evals =  evals, sdeg = sdeg)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f && blkdims == [0, 4, 0] &&
      sort(real( evals)) ≈ sort(real(SF.values)) && sort(imag( evals)) ≈ sort(imag(SF.values))
sdeg, evals, eigvals(a+b*f)

##  SB01DD EXAMPLE PROGRAM DATA
#   for this example PLACE fails
a = [
  -1.0  0.0  2.0 -3.0
   1.0 -4.0  3.0 -1.0
   0.0  2.0  4.0 -5.0
   0.0  0.0 -1.0 -2.0];
b = [
   1.0  0.0
   0.0  0.0
   0.0  0.0
   0.0  1.0];
 evals = [
  -1.0, -1.0, -1.0, -1.0];
sdeg = [];
@time f, SF, blkdims = saloc(a,b, evals =  evals)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f && blkdims == [0, 4, 0] &&
      sort(real( evals)) ≈ sort(real(SF.values)) && sort(imag( evals)) ≈ sort(imag(SF.values))
 evals, eigvals(a+b*f)

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

a = rand(Ty,6,6); b = rand(Ty,6,3);   evals = eigvals(a);
@time f, SF, blkdims = saloc(a,b, evals =  evals)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f  && norm(f) < 1.e-10 && blkdims == [0, 6, 0] &&
      sort(real( evals)) ≈ sort(real(SF.values)) && sort(imag( evals)) ≈ sort(imag(SF.values))

# random example
a = rand(Ty,6,6); b = [zeros(Ty,2,3);rand(Ty,2,3);zeros(Ty,2,3)];   evals = eigvals(a); evals = evals[sortperm(rand(6))];
@time f, SF, blkdims = saloc(a,b, evals =  evals)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f  && norm(f) < 1.e-10 && blkdims == [0, 6, 0] &&
      sort(real( evals)) ≈ sort(real(SF.values)) && sort(imag( evals)) ≈ sort(imag(SF.values))

# random example
a = rand(Ty,6,6); c = [zeros(Ty,3,2) rand(Ty,3,2) zeros(Ty,3,2)];   evals = eigvals(rand(Ty,6,6));
@time k, SF, blkdims = salocd(a, c, evals =  evals)
@test SF.Z*SF.T*SF.Z' ≈ a+k*c  && blkdims == [0, 6, 0] &&
      sort(real( evals)) ≈ sort(real(SF.values)) && sort(imag( evals)) ≈ sort(imag(SF.values))

# random example  - stabilization
a = randn(Ty,6,6); b = rand(Ty,6,3);   sdeg = -0.2;
@time f, SF, blkdims = saloc(a,b, sdeg = sdeg)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f  && sum(blkdims[1:2]) == 6 &&
      all(real(SF.values) .<= sdeg+sqrt(eps()))

# random example  - stabilization
a = randn(Ty,6,6); b = rand(Ty,6,3);
@time f, SF, blkdims = saloc(a,b)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f  && sum(blkdims[1:2]) == 6 &&
      all(real(SF.values) .<= 0)

# random example  - stabilization
a = rand(Ty,6,6); b = rand(Ty,6,3);   sdeg = 0.2;
@time f, SF, blkdims = saloc(a,b, sdeg = sdeg, disc = true)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f  && sum(blkdims[1:2]) == 6 &&
      all(abs.(SF.values) .<= sdeg+sqrt(eps()))

# random example  - stabilization
a = randn(Ty,6,6); b = rand(Ty,6,3);
@time f, SF, blkdims = saloc(a,b, disc = true)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f  && sum(blkdims[1:2]) == 6 &&
      all(abs.(SF.values) .<= 1)

# random example - stabilization
a = randn(Ty,6,6); b = rand(Ty,6,3);   sdeg = -0.2;
@time f, SF, blkdims = saloc(a,b, sdeg = sdeg)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f  && sum(blkdims[1:2]) == 6 &&
      all(real(SF.values) .<= sdeg+sqrt(eps()))

# random example - stabilization
a = randn(Ty,6,6); b = rand(Ty,6,3);
@time f, SF, blkdims = saloc(a,b)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f  && sum(blkdims[1:2]) == 6 &&
      all(real(SF.values) .<= 0)

# random example - stabilization
a = randn(Ty,6,6); b = rand(Ty,6,3);
@time f, SF, blkdims = saloc(a,b,evals = [-1, -2])
@test SF.Z*SF.T*SF.Z' ≈ a+b*f  && sum(blkdims[1:2]) == 6 &&
      all(real(SF.values) .<= 0)

# random example - stabilization
a = randn(Ty,6,6); b = rand(Ty,6,3);
@time f, SF, blkdims = saloc(a,b,evals = [-1])
@test SF.Z*SF.T*SF.Z' ≈ a+b*f  && sum(blkdims[1:2]) == 6 &&
      all(real(SF.values) .<= 0)

# random example - stabilization
a = rand(Ty,6,6); b = rand(Ty,6,3);  sdeg = 0.2
@time f, SF, blkdims = saloc(a,b, sdeg = sdeg, disc = true)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f  && sum(blkdims[1:2]) == 6 &&
      all(abs.(SF.values) .<= sdeg+sqrt(eps()))

a = rand(Ty,6,6); b = rand(Ty,6,3);
@time f, SF, blkdims = saloc(a,b, disc = true)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f  && sum(blkdims[1:2]) == 6 &&
      all(abs.(SF.values) .<= 1)

##  uncontrollable system
n = 10; nc = 6; nu = n-nc; m = 4;
au = rand(Ty,nu,nu); ac = rand(Ty,nc,nc);
a = [ac rand(Ty,nc,nu); zeros(Ty,nu,nc) au]; b = [rand(Ty,nc,m);zeros(Ty,nu,m)];
q = qr(rand(Ty,n,n)).Q;
poluc = eigvals(au);
a = q'*a*q; b = q'*b;
#evals = complex(rand(nc))
evals = eigvals(rand(Ty,nc,nc))
@time f, SF, blkdims = saloc(a,b, evals =  evals, atol1 = 1.e-7, atol2 = 1.e-7)
@test SF.Z*SF.T*SF.Z' ≈ a+b*f  && blkdims == [0, nc, nu] &&
      sort(real(poluc)) ≈ sort(real(SF.values[nc+1:n])) && sort(imag(poluc)) ≈ sort(imag(SF.values[nc+1:n])) &&
      sort(real(evals)) ≈ sort(real(SF.values[1:nc])) && sort(imag( evals)) ≈ sort(imag(SF.values[1:nc]))
end # Ty loop

end #begin

println("saloc for generalized pair (A-λE,B)")
@testset "saloc for generalized pair (A-λE,B)" begin

## simple cases: B = [];
a = ones(1,1); e = I; b = zeros(1,0); evals = [-1];
@time f, SF, blkdims = saloc(a,e,b,evals = evals)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e  && blkdims == [0, 0, 0, 1]
evals,eigvals(a+b*f)

## simple cases: B = [];
a = ones(1,1); e = 2*ones(1,1); b = zeros(1,0); evals = [-1];
@time f, SF, blkdims = saloc(a,e,b,evals = evals)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e  && blkdims == [0, 0, 0, 1]
evals,eigvals(a+b*f,e)

## simple cases
a = ones(1,1); e = 2*ones(1,1); b = zeros(1,5); evals = [-1];
@time f, SF, blkdims = saloc(a,e,b,evals = evals)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e  && blkdims == [0, 0, 0, 1]
evals,eigvals(a+b*f,e)

## simple cases
a = ones(1,1); e = 2*ones(1,1); b = ones(1,1); evals = [-1];
@time f, SF, blkdims = saloc(a,e,b,evals = evals)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e  && blkdims == [0, 0, 1, 0]
@test sort(real(evals)) ≈ sort(real(SF.α ./SF.β)) && sort(imag(evals)) ≈ sort(imag(SF.α ./SF.β))
evals,eigvals(a+b*f,e)

## simple cases
a = ones(1,1); e = 2*ones(1,1); b = ones(1,1); evals = [-1];
@time f, SF, blkdims = saloc(a,e,b,evals = evals,sepinf = false)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e  && blkdims == [0, 0, 1, 0] &&
    sort(real(evals)) ≈ sort(real(SF.α ./SF.β)) && sort(imag(evals)) ≈ sort(imag(SF.α ./SF.β))
evals,eigvals(a+b*f,e)

## simple cases
a = complex(ones(1,1)); e = 2*ones(1,1); b = ones(1,1); evals = [-im];
@time f, SF, blkdims = saloc(a,e,b,evals = evals)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e  && blkdims == [0, 0, 1, 0] &&
      sort(real(evals)) ≈ sort(real(SF.α ./SF.β)) && sort(imag(evals)) ≈ sort(imag(SF.α ./SF.β))
evals,eigvals(a+b*f,e)

## simple cases
a = complex(ones(1,1)); e = 2*ones(1,1); b = ones(1,1); evals = [-im];
@time f, SF, blkdims = saloc(a,e,b,evals = evals,sepinf = false)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e  && blkdims == [0, 0, 1, 0] &&
      sort(real(evals)) ≈ sort(real(SF.α ./SF.β)) && sort(imag(evals)) ≈ sort(imag(SF.α ./SF.β))
evals,eigvals(a+b*f,e)

## simple cases (no real eigenvalue)
a = ones(1,1); e = 2*ones(1,1); b = ones(1,1); evals = [-1+im*0.5;-1-im*0.5];
@time f, SF, blkdims = saloc(a,e,b,evals = evals)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e  && blkdims == [0, 0, 1, 0]
evals,eigvals(a+b*f,e)

## simple cases (no real eigenvalue, use sdeg)
a = ones(1,1); e = 2*ones(1,1);  b = ones(1,1); evals = [-1+im*0.5;-1-im*0.5]; sdeg = -0.5;
@time f, SF, blkdims = saloc(a,e,b,evals = evals, sdeg = sdeg)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e  && blkdims == [0, 0, 1, 0]
evals,eigvals(a+b*f,e)

## simple cases
a = ones(1,1); e = 2*ones(1,1);  b = ones(1,1); evals = [-1;-0.5];
@time f, SF, blkdims = saloc(a,e,b,evals = evals)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e  && blkdims == [0, 0, 1, 0]
evals,eigvals(a+b*f,e)

## simple cases
a = [1 0;0 1]; e = rand(2,2); b = rand(2,5); evals = [-1;-0.5];
e =  [0.0568026  0.988833
0.0651635  0.862077];
b = [0.338582  0.74765  0.314715  0.265301  0.749282
0.221894  0.44414  0.910163  0.994525  0.0763453]
@time f, SF, blkdims = saloc(a,e,b,evals = evals)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e && blkdims == [0, 0, 2, 0]
@test sort(real(evals)) ≈ sort(real(SF.α ./SF.β)) && sort(imag(evals)) ≈ sort(imag(SF.α ./SF.β))
evals,eigvals(a+b*f,e)

## simple cases
a = [1 0;0 1]; e = rand(2,2); b = rand(2,5); evals = [-1;-0.5];
e =  [0.0568026  0.988833
0.0651635  0.862077];
b = [0.338582  0.74765  0.314715  0.265301  0.749282
0.221894  0.44414  0.910163  0.994525  0.0763453]
@time f, SF, blkdims = saloc(a,e,b,evals = evals, sepinf = false)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e && blkdims == [0, 0, 2, 0] &&
      sort(real(evals)) ≈ sort(real(SF.α ./SF.β)) && sort(imag(evals)) ≈ sort(imag(SF.α ./SF.β))

## simple cases
a = [1 0;0 1]; e = rand(2,2); b = rand(2,5); evals = [-1;-0.5];
e =  [0.0568026  0.988833
0  0.862077];
b = [0.338582  0.74765  0.314715  0.265301  0.749282
0.221894  0.44414  0.910163  0.994525  0.0763453]
@time f, SF, blkdims = saloc(a,e,b,evals = evals, sepinf = false)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e && blkdims == [0, 0, 2, 0] &&
      sort(real(evals)) ≈ sort(real(SF.α ./SF.β)) && sort(imag(evals)) ≈ sort(imag(SF.α ./SF.β))

## simple cases  uncontrollable pair
a = rand(2,2); e = a;  b = rand(2,1);  evals = [-1+im*0.5;-1-im*0.5];
@time f, SF, blkdims = saloc(a, e, b, evals =  evals)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e && blkdims == [0, 0, 1, 1]
evals,eigvals(a+b*f,e)

## simple cases
a = [1 1;0 2]; e = rand(2,2); b = rand(2,5); evals = [-1+im*0.5;-1-im*0.5];
@time f, SF, blkdims = saloc(a,e,b,evals = evals)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e && blkdims == [0, 0, 2, 0] &&
      sort(real(evals)) ≈ sort(real(SF.α ./SF.β)) && sort(imag(evals)) ≈ sort(imag(SF.α ./SF.β))
evals,eigvals(a+b*f,e)

## simple cases
e = rand(3,3); a = e*[1 1 0;-1 1 0;0 0 2];  b = e*[0 0; 0 0; 1 1];  evals = [-1+im*0.5;-1-im*0.5];
@time f, SF, blkdims = saloc(a, e, b, evals =  evals)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e && blkdims == [0, 0, 1, 2]
evals,eigvals(a+b*f,e)

## simple cases
e = rand(3,3); a = e*[1 1 1;-1 1 0;0 0 2]; b = rand(3,5); evals = [-1+im*0.5;-1-im*0.5;-2];
@time f, SF, blkdims = saloc(a,e,b,evals = evals)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e && blkdims == [0, 0, 3, 0] &&
      sort(real(evals)) ≈ sort(real(SF.α ./SF.β)) && sort(imag(evals)) ≈ sort(imag(SF.α ./SF.β))
evals,eigvals(a+b*f,e)

## simple cases
e = [1 1 1;-1 1 0;0 0 0]; a = [1 1 1;-1 1 0;0 0 2]; b = rand(3,5); evals = [-1+im*0.5;-1-im*0.5;-2];
@time f, SF, blkdims = saloc(a,e,b,evals = evals)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e && blkdims == [1, 0, 2, 0]
evals,eigvals(a+b*f,e)

## simple cases
e = [1 1 1;-1 1 0;0 0 0]; a = [1 1 1;-1 1 0;0 0 2]; b = rand(3,5); evals = [-1+im*0.5;-1-im*0.5;-2];
@time f, SF, blkdims = saloc(a,e,b,evals = evals, sepinf = false)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e && blkdims == [1, 0, 2, 0]
evals,eigvals(a+b*f,e)

## simple cases
#e = rand(3,3);
e = [
 0.521866  0.369606  0.456723
 0.740103  0.357063  0.280376
 0.397995  0.362346  0.885047];
a = e*[1 1 1;-1 1 0;0 0 2]; b = rand(3,5); evals = [-1;-0.5;-2];
@time f, SF, blkdims = saloc(a,e,b,evals = evals)
#@time f, SF, blkdims = saloc(a,e,b,evals = evals, sepinf = false)

@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e && blkdims == [0, 0, 3, 0] &&
      sort(real(evals)) ≈ sort(real(SF.α ./SF.β)) && sort(imag(evals)) ≈ sort(imag(SF.α ./SF.β))

## simple cases
e = [1 1 1;-1 1 0;0 0 0]; a = [1 1 1;-1 1 0;0 0 2]; b = rand(3,5); evals = [-1;-0.5;-2];
@time f, SF, blkdims = saloc(a,e,b,evals = evals, sepinf = false)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e && blkdims == [1, 0, 2, 0]
evals,eigvals(a+b*f,e)

## simple cases
e = rand(3,3); a = e*[1 1 1;-1 1 0;0 0 2]; b = rand(3,5); sdeg = -0.5;
@time f, SF, blkdims = saloc(a,e,b,sdeg = sdeg)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e  && blkdims == [0, 0, 3, 0] &&
      all(sdeg .≈ (real(SF.values)))
sdeg, eigvals(a+b*f,e)

## simple cases
e = rand(3,3); a = e*[1 1 1;-1 1 0;0 0 2];  b = rand(3,5); evals = [-1,-0.5]; sdeg = -0.5;
@time f, SF, blkdims = saloc(a,e,b,evals = evals, sdeg = sdeg)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e  && blkdims == [0, 0, 3, 0]
sdeg, eigvals(a+b*f,e)

## simple cases
e = rand(3,3); a = e*[1 1 1;-1 1 0;0 0 2]; b = rand(3,5); evals = [-3]; sdeg = -2;
@time f, SF, blkdims = saloc(a,e,b,evals = evals, sdeg = sdeg)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e && blkdims == [0, 0, 3, 0]
sdeg, eigvals(a+b*f,e)

## simple cases
e = rand(3,3); a = e*[1 1 0;-1 1 0;0 0 2];  b = e*[0 0; 0 0; 1 1];  evals = [-1+im*0.5;-1-im*0.5];
@time f, SF, blkdims = saloc(a, e, b, atol1 = 1.e-7, atol2 = 1.e-7, atol3 = 1.e-7, evals =  evals)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e  && blkdims == [0, 0, 1, 2]
evals,eigvals(a+b*f,e)

## SB01BD EXAMPLE PROGRAM DATA
a = [
  -6.8000   0.0000  -207.0000   0.0000
   1.0000   0.0000     0.0000   0.0000
  43.2000   0.0000     0.0000  -4.2000
   0.0000   0.0000     1.0000   0.0000];
b = [
   5.6400   0.0000
   0.0000   0.0000
   0.0000   1.1800
   0.0000   0.0000];
evals = [
  -0.5000+im*0.1500
  -0.5000-im*0.1500
  -2.0000
  -0.4000];
sdeg = -0.4;
e = rand(4,4); a = e*a; b = e*b;
@time f, SF, blkdims = saloc(a,e,b,evals = evals, sdeg = sdeg)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e  && blkdims == [0, 2, 2, 0]
sdeg,evals, eigvals(a+b*f,e)

sdeg = missing;
@time f, SF, blkdims = saloc(a,e,b,evals = evals, sdeg = sdeg)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e && blkdims == [0, 0, 4, 0] &&
      sort(real(evals)) ≈ sort(real(SF.values)) && sort(imag(evals)) ≈ sort(imag(SF.values))
sdeg, evals, eigvals(a+b*f,e)

##  SB01DD EXAMPLE PROGRAM DATA
#   for this example PLACE fails
a = [
  -1.0  0.0  2.0 -3.0
   1.0 -4.0  3.0 -1.0
   0.0  2.0  4.0 -5.0
   0.0  0.0 -1.0 -2.0];
b = [
   1.0  0.0
   0.0  0.0
   0.0  0.0
   0.0  1.0];
evals = [ -1.0, -1.0, -1.0, -1.0];
sdeg = missing;
e = rand(4,4); a = e*a; b = e*b;
@time f, SF, blkdims = saloc(a,e,b,evals = evals)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e  && blkdims == [0, 0, 4, 0]
@test sort(real(evals)) ≈ sort(real(SF.values)) && sort(imag(evals)) ≈ sort(imag(SF.values))
sdeg, evals, eigvals(a+b*f,e)

# random example
fast = true
for fast in (true,false)

Ty = Float64
for Ty in (Float64, Complex{Float64})
complx = Ty <: Complex

a = rand(Ty,6,6); b = rand(Ty,6,3); e = rand(Ty,6,6);
evals = complx ? eigvals(a,e) : evsym!(eigvals(a,e));
@time f, SF, blkdims = saloc(a,e,b,evals = evals, fast = fast)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e  && norm(f) < 1.e-10  && blkdims == [0, 0, 6, 0] &&
      sort(real(evals)) ≈ sort(real(SF.values)) && sort(imag(evals)) ≈ sort(imag(SF.values))

# random example
a = rand(Ty,6,6); b = [zeros(Ty,2,3);rand(Ty,2,3);zeros(Ty,2,3)]; e = rand(Ty,6,6);
evals = complx ? eigvals(a,e) : evsym!(eigvals(a,e));
@time f, SF, blkdims = saloc(a,e,b,evals = evals, fast = fast)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e  && norm(f) < 1.e-10  && blkdims == [0, 0, 6, 0] &&
      sort(real(evals)) ≈ sort(real(SF.values)) && sort(imag(evals)) ≈ sort(imag(SF.values))

# random example  - stabilization
a = randn(Ty,6,6); e = randn(Ty,6,6); b = rand(Ty,6,3);
@time f, SF, blkdims = saloc(a,e,b,evals = [-1, -2])
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e  && sum(blkdims[2:3]) == 6  &&
      all(real(SF.values) .<= 0)

@time f, SF, blkdims = saloc(a,b,evals = [-1, -2])

# random example  - stabilization
a = randn(Ty,6,6); e = randn(Ty,6,6); b = rand(Ty,6,3);
@time f, SF, blkdims = saloc(a,e,b,disc = true, sepinf = false)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e  && sum(blkdims[2:3]) == 6  &&
      all(abs.(SF.values) .<= 1)

# random example  - stabilization
a = randn(Ty,6,6); e = randn(Ty,6,6); b = rand(Ty,6,3);
@time f, SF, blkdims = saloc(a,e,b,disc = true, sepinf = true)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e  && sum(blkdims[2:3]) == 6  &&
      all(abs.(SF.values) .<= 1)

# random example  - infinite eigenvalue assignment
a = randn(Ty,6,6); e = randn(Ty,6,6); b = rand(Ty,6,3);
@time f, g, SF, blkdims = salocinf(a,e,b,atol1=1.e-7,atol2=1.e-7)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e+b*g  && sum(blkdims[2:3]) == 6  &&
      any(isfinite.(SF.values)) == false

# random example  - infinite eigenvalue assignment
a = triu(randn(Ty,6,6),1); e = randn(Ty,6,6); b = rand(Ty,6,3);
@time f, g, SF, blkdims = salocinf(a,e,b,atol1=1.e-7,atol2=1.e-7)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e+b*g  && blkdims == [0, 6, 0]  &&
      any(isfinite.(SF.values)) == false

# random example
a = rand(Ty,6,6); c = [zeros(Ty,3,2) rand(Ty,3,2) zeros(Ty,3,2)]; e = rand(Ty,6,6); evals = eigvals(rand(Ty,6,6),rand(Ty,6,6));
complx || evsym!(evals)
@time k, SF, blkdims = salocd(a,e,c,evals = evals, fast = fast)
@test SF.Q*SF.S*SF.Z' ≈ a+k*c && SF.Q*SF.T*SF.Z' ≈ e  && #norm(k) < 1.e-10  &&
      blkdims == [0, 6, 0, 0] &&
      sort(real(evals)) ≈ sort(real(SF.values)) && sort(imag(evals)) ≈ sort(imag(SF.values))

# random example
a = triu(rand(Ty,6,6),1); c = [zeros(Ty,3,2) rand(Ty,3,2) zeros(Ty,3,2)]; e = rand(Ty,6,6);
@time k, l, SF, blkdims = salocinfd(a, e, c, fast = fast, atol1=1.e-7, atol2 = 1.e-7, atol3=1.e-7)
@test SF.Q*SF.S*SF.Z' ≈ a+k*c && SF.Q*SF.T*SF.Z' ≈ e+l*c  &&  blkdims == [1, 5, 0] &&
      count(isfinite.(SF.values)) == blkdims[1]

##  uncontrollable system
n = 10; nc = 6; nu = n-nc; m = 4;
au = rand(Ty,nu,nu); ac = rand(Ty,nc,nc); eu = rand(Ty,nu,nu); ec = rand(Ty,nc,nc);
a = [ac rand(Ty,nc,nu); zeros(Ty,nu,nc) au]; b = [rand(Ty,nc,m);zeros(Ty,nu,m)];
e = [ec rand(Ty,nc,nu); zeros(Ty,nu,nc) eu];
q = qr(rand(Ty,n,n)).Q; z = qr(rand(Ty,n,n)).Q;
poluc = eigvals(au,eu);
a = q'*a*z; e = q'*e*z; b = q'*b;
evals = eigvals(rand(Ty,nc,nc));
@time f, SF, blkdims = saloc(a,e,b,evals = evals, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7, atol3 = 1.e-7)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e  && blkdims == [0, 0, nc, nu] &&
      sort(real(poluc)) ≈ sort(real(SF.values[nc+1:n])) && sort(imag(poluc)) ≈ sort(imag(SF.values[nc+1:n])) &&
      sort(real(evals)) ≈ sort(real(SF.values[1:nc])) && sort(imag(evals)) ≈ sort(imag(SF.values[1:nc]))

# n = 10; nc = 6; nu = n-nc; m = 4;
# au = rand(Ty,nu,nu); ac = rand(Ty,nc,nc); eu = rand(Ty,nu,nu); ec = rand(Ty,nc,nc);
# a = [ac rand(Ty,nc,nu); zeros(Ty,nu,nc) au]; b = [rand(Ty,nc,m);zeros(Ty,nu,m)];
# e = [ec rand(Ty,nc,nu); zeros(Ty,nu,nc) eu];
# q = qr(rand(Ty,n,n)).Q; z = qr(rand(Ty,n,n)).Q;
# poluc = eigvals(au,eu);
# a = q'*a*z; e = q'*e*z; b = q'*b;
# @time f, g, SF, blkdims = salocinf(a,e,b,fast = fast, atol1 = 1.e-7, atol2 = 1.e-7, atol3 = 1.e-7)
# @test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e+b*g  && blkdims == [0, nc, nu] &&
#       sort(real(poluc)) ≈ sort(real(SF.values[nc+1:n])) && sort(imag(poluc)) ≈ sort(imag(SF.values[nc+1:n])) &&
#       any(isfinite.(SF.values[1:nc])) == false

end # Ty loop

for sepinf in (true,false)

## proper system
n = 10;  ns = 3; m = 4;
e = [zeros(n,ns) randn(n,n-ns)]'; a = randn(n,n); b = randn(n,m);
evals = complex(randn(n-ns))
nc = 2*floor(Int,floor((n-ns+2)/2)*rand())
for ii = 1 : 2 : nc,
   evals[ii]   = evals[ii] + im*evals[ii+1];
   evals[ii+1] = conj( evals[ii] );
end
@time f, SF, blkdims = saloc(a,e,b,evals = evals, sepinf = sepinf, fast = fast)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e  && blkdims == [ns, 0, n-ns, 0] &&
      sort(real(evals)) ≈ sort(real(SF.values[ns+1:n])) && sort(imag(evals)) ≈ sort(imag(SF.values[ns+1:n]))

## Simas' paper (standard system)
n = 10; m = 4;
a = randn(n,n); b = randn(n,m);
evals = complex(randn(n));
nc = 2*floor(Int,floor((n+2)/2)*rand())
for ii = 1 : 2 : nc,
   evals[ii]   = evals[ii] + im*evals[ii+1];
   evals[ii+1] = conj( evals[ii] );
end
@time f, SF, blkdims = saloc(a,I,b,evals = evals, sepinf = sepinf, fast = fast);
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ I  && blkdims == [0, 0, n, 0] &&
      sort(real(evals)) ≈ sort(real(SF.values[1:n])) && sort(imag(evals)) ≈ sort(imag(SF.values[1:n]))

## Sima's paper (descriptor system)
n = 10;  m = 4;
a = randn(n,n);  e = randn(n,n);  b = randn(n,m);
evals = complex(randn(n));
nc = 2*floor(Int,floor((n+2)/2)*rand())
for ii = 1 : 2 : nc,
   evals[ii]   = evals[ii] + im*evals[ii+1];
   evals[ii+1] = conj( evals[ii] );
end
@time f, SF, blkdims = saloc(a,e,b,evals = evals, sepinf = sepinf, fast = fast)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e  && blkdims == [0, 0, n, 0] &&
      sort(real(evals)) ≈ sort(real(SF.values[1:n])) && sort(imag(evals)) ≈ sort(imag(SF.values[1:n]))

## improper system R-stabilization
ni = [3, 2, 1, 1]; nf = 5; ninf = sum(ni); n = ninf+nf; m = 4;
a = [triu(rand(ninf,n)); zeros(nf,ninf) rand(nf,nf)];
ei = rand(ninf,n); b = randn(n,m);
k1 = 0;
k2 = k1 + ni[1]
ei[k1+1:ninf,1:k2] .= 0
k1 = k2
k2 = k1 + ni[2]
ei[k1+1:ninf,1:k2] .= 0
k1 = k2
k2 = k1 + ni[3]
ei[k1+1:ninf,1:k2] .= 0
k1 = k2
k2 = k1 + ni[4]
ei[k1+1:ninf,1:k2] .= 0

e = [ei; zeros(nf,ninf) rand(nf,nf)];
q = qr(rand(n,n)).Q; z = qr(rand(n,n)).Q;
a = q*a*z; e = q*e*z;

evals = complex(-rand(nf) .- 1);
nc = 2*floor(Int,floor((n-ninf+2)/2)*rand())
for ii = 1 : 2 : nc,
   evals[ii]   = evals[ii] + im*evals[ii+1];
   evals[ii+1] = conj( evals[ii] );
end
@time f, SF, blkdims = saloc(a,e,b,evals = evals, fast = fast, sepinf = sepinf, atol1 = 1.e-7, atol2 = 1.e-7, atol3 = 1.e-7)
@test SF.Q*SF.S*SF.Z' ≈ a+b*f && SF.Q*SF.T*SF.Z' ≈ e  && blkdims == [ninf, 0, nf, 0] &&
      sort(real(evals)) ≈ sort(real(SF.values[ninf+1:n])) && sort(imag(evals)) ≈ sort(imag(SF.values[ninf+1:n]))

end # sepinf loop
end # fast loop

end #begin

end #test

end #module


text:
module NamedDims
using Base: @propagate_inbounds
using Base.Broadcast:
    Broadcasted, BroadcastStyle, DefaultArrayStyle, AbstractArrayStyle, Unknown
using LinearAlgebra
using Statistics

export NamedDimsArray, dim, rename

include("name_core.jl")
include("wrapper_array.jl")
include("broadcasting.jl")
include("functions.jl")
include("functions_dims.jl")
include("functions_math.jl")

include("compat.jl")
end # module

ext: jl


text:
using LocalRegistry, Pkg

const __topdir__ = dirname(@__DIR__)

"""
    bump_version(version::VersionNumber)

Bump the version number of all packages in the repository to `version`.
"""
function bump_version(version::VersionNumber)
    for dir in (__topdir__,
                joinpath(__topdir__, "SnoopCompileCore"))
        projfile = joinpath(dir, "Project.toml")
        lines = readlines(projfile)
        idxs = findall(str->startswith(str, "version"), lines)
        idx = only(idxs)
        lines[idx] = "version = \"$version\""
        # Find any dependencies on versions within the repo
        idxs = findall(str->str == "[compat]", lines)
        if !isempty(idxs)
            idxcompat = only(idxs)
            idxend = findnext(str->!isempty(str) && str[1] == '[', lines, idxcompat+1)  # first line of next TOML section
            if idxend === nothing
                idxend = length(lines) + 1
            end
            for i = idxcompat:idxend-1
                if startswith(lines[i], "SnoopCompile")
                    strs = split(lines[i], '=')
                    lines[i] = strs[1] * "= \"~" * string(version) * '"'
                end
            end
        end
        open(projfile, "w") do io
            for line in lines
                println(io, line)
            end
        end
    end
    return nothing
end

function register_all()
    for pkg in ("SnoopCompileCore", "SnoopCompile")
        pkgsym = Symbol(pkg)
        @eval Main using $pkgsym
        register(getfield(Main, pkgsym)::Module)
    end
end


text:
# Toy model from paper
include("../../src/SFFM.jl")
using LinearAlgebra
# import SFFM

## Define the model
T = [-1.0 1.0; 1.0 -1.0]
C = [1.0; -2.0]
r = (
    r = function (x)
        [(x .<= 1)-(x .> 1) ((x .> 1))-2*(x .<= 1).*(x .> 0)-(x.==0)]
    end,
    R = function (x)
        [x.*((x .<= 1)-(x .> 1)) x.*(((x .> 1))-2*(x .<= 1).*(x .> 0)-(x.==0))]
    end,
)

Bounds = [0 1.8; -Inf Inf]
model = SFFM.Model(T, C, r, Bounds = Bounds)

## Define mesh
Nodes = [0;1;1.8]
nBases = 2
Basis = "lagrange"
mesh = SFFM.DGMesh(model, Nodes, nBases, Basis=Basis)

## Make matrices
# All = SFFM.MakeAll( model, mesh)
# M = SFFM.MakeMatrices( model, mesh)
B = SFFM.MakeFullGenerator(model, mesh)
R = SFFM.MakeR( model, mesh)

println("")
println("The operator B looks like: ")
QBDidx = SFFM.MakeQBDidx(model, mesh)
display(Matrix(B.B[QBDidx,QBDidx]))

println("")
println("The subarrays of B look like this:")
println("for 11++")
display(Matrix(B[("+","+"),((1,1))]))
println("")
println("for 11+-")
display(Matrix(B[("+","-"),(1,1)]))
println("")
println("for 11-+")
display(Matrix(B[("-","+"),(1,1)]))
println("")
println("for 11--")
display(Matrix(B[("-","-"),(1,1)]))

println("")
println("for 12++")
display(Matrix(B[("+","+"),(1,2)]))
println("")
println("for 12+-")
display(Matrix(B[("+","-"),(1,2)]))
println("")
println("for 12-+")
display(Matrix(B[("-","+"),(1,2)]))
println("")
println("for 12--")
display(Matrix(B[("-","-"),(1,2)]))

println("")
println("for 21++")
display(Matrix(B[("+","+"),(2,1)]))
println("")
println("for 21+-")
display(Matrix(B[("+","-"),(2,1)]))
println("")
println("for 21-+")
display(Matrix(B[("-","+"),(2,1)]))
println("")
println("for 21--")
display(Matrix(B[("-","-"),(2,1)]))

println("")
println("for 22++")
display(Matrix(B[("+","+"),(2,2)]))
println("")
println("for 22+-")
display(Matrix(B[("+","-"),(2,2)]))
println("")
println("for 22-+")
display(Matrix(B[("-","+"),(2,2)]))
println("")
println("for 22--")
display(Matrix(B[("-","-"),(2,2)]))

plusIdx = QBDidx[[mesh.Fil["p+",:];repeat(mesh.Fil["+",:]',2)[:];mesh.Fil["q+",:]]]
minusIdx = QBDidx[[mesh.Fil["p-",:];repeat(mesh.Fil["-",:]',2)[:];mesh.Fil["q-",:]]]
println("")
println("for ++")
display(Matrix(B.B)[plusIdx,plusIdx])
println("")
println("for +-")
display(Matrix(B.B)[plusIdx,minusIdx])
println("")
println("for -+")
display(Matrix(B.B)[minusIdx,plusIdx])
println("")
println("for --")
display(Matrix(B.B)[minusIdx,minusIdx])

println("The operator R⁺ is ")
display(Matrix(R.RDict["+",:]))
println("")
println("The operator R⁻ is ")
display(Matrix(R.RDict["-",:]))

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
using FromFile
using Random: shuffle!
@from "Core.jl" import Options, Dataset, RecordType
@from "EquationUtils.jl" import stringTree
@from "PopMember.jl" import PopMember
@from "Population.jl" import Population, bestOfSample
@from "Mutate.jl" import nextGeneration
@from "Recorder.jl" import @recorder

# Pass through the population several times, replacing the oldest
# with the fittest of a small subsample
function regEvolCycle(dataset::Dataset{T},
                      baseline::T, pop::Population, temperature::T, curmaxsize::Int,
                      frequencyComplexity::AbstractVector{T},
                      options::Options,
                      record::RecordType)::Population where {T<:Real}
    # Batch over each subsample. Can give 15% improvement in speed; probably moreso for large pops.
    # but is ultimately a different algorithm than regularized evolution, and might not be
    # as good.
    if options.fast_cycle

        # These options are not implemented for fast_cycle:
        @recorder error("You cannot have the recorder and fast_cycle set to true at the same time!")
        @assert options.probPickFirst == 1.0

        shuffle!(pop.members)
        n_evol_cycles = round(Int, pop.n/options.ns)
        babies = Array{PopMember}(undef, n_evol_cycles)

        # Iterate each ns-member sub-sample
        @inbounds Threads.@threads for i=1:n_evol_cycles
            best_score = Inf
            best_idx = 1+(i-1)*options.ns
            # Calculate best member of the subsample:
            for sub_i=1+(i-1)*options.ns:i*options.ns
                if pop.members[sub_i].score < best_score
                    best_score = pop.members[sub_i].score
                    best_idx = sub_i
                end
            end
            allstar = pop.members[best_idx]
            mutation_recorder = RecordType()
            babies[i] = nextGeneration(dataset, baseline, allstar, temperature,
                                       curmaxsize, frequencyComplexity, options,
                                       tmp_recorder=mutation_recorder)
        end

        # Replace the n_evol_cycles-oldest members of each population
        @inbounds for i=1:n_evol_cycles
            oldest = argmin([pop.members[member].birth for member=1:pop.n])
            pop.members[oldest] = babies[i]
        end
    else
        for i=1:round(Int, pop.n/options.ns)
            allstar = bestOfSample(pop, options)
            mutation_recorder = RecordType()
            baby = nextGeneration(dataset, baseline, allstar, temperature,
                                  curmaxsize, frequencyComplexity, options,
                                  tmp_recorder=mutation_recorder)

            oldest = argmin([pop.members[member].birth for member=1:pop.n])

            @recorder begin
                if !haskey(record, "mutations")
                    record["mutations"] = RecordType()
                end
                for member in [allstar, baby, pop.members[oldest]]
                    if !haskey(record["mutations"], "$(member.ref)")
                        record["mutations"]["$(member.ref)"] = RecordType("events"=>Vector{RecordType}(),
                                                                          "tree"=>stringTree(member.tree, options),
                                                                          "score"=>member.score,
                                                                          "parent"=>member.parent)
                    end
                end
                mutate_event = RecordType("type"=>"mutate", "time"=>time(), "child"=>baby.ref, "mutation"=>mutation_recorder)
                death_event  = RecordType("type"=>"death",  "time"=>time())

                # Put in random key rather than vector; otherwise there are collisions!
                push!(record["mutations"]["$(allstar.ref)"]["events"], mutate_event)
                push!(record["mutations"]["$(pop.members[oldest].ref)"]["events"], death_event)
            end

            pop.members[oldest] = baby
        end
    end

    return pop
end


text:
function handle_step(obj::OrekitStepHandler, arg0::OrekitStepInterpolator, arg1::jboolean)
    return jcall(obj, "handleStep", void, (OrekitStepInterpolator, jboolean), arg0, arg1)
end

function init(obj::OrekitStepHandler, arg0::SpacecraftState, arg1::AbsoluteDate)
    return jcall(obj, "init", void, (SpacecraftState, AbsoluteDate), arg0, arg1)
end


text:
__precompile__()
module ReduceLinAlg
using Reduce, LinearAlgebra

#   This file is part of ReduceLinAlg.jl. It is licensed under the MIT license
#   Copyright (C) 2018 Michael Reed

const VectorAny = Union{Vector,Adjoint}

const lin = [
    :hessian,
    :mat_jacobian,
    :band_matrix,
    :block_matrix,
    :char_matrix,
    :char_poly,
    :diagonal,
    :extend,
    :gram_schmidt,
    :hilbert,
    :jordan_block,
    :kronecker_product,
]

const unr = [
    :cholesky,
    :coeff_matrix,
    :hermitian_tp,
    :symmetricp,
    :toeplitz,
    :triang_adjoint,
    :vandermonde,
]

:(export $([lin;unr]...)) |> eval

for fun in lin
    @eval begin
        $(Reduce.parsegen(fun,:args))
        $(Reduce.unfoldgen(fun,:args))
        function $fun(expr::String,s...;be=0)
            convert(String, $fun(RExpr(expr),s...;be=be))
        end
    end
end

for fun in unr
    @eval begin
        $(Reduce.parsegen(fun,:unary))
        $(Reduce.unfoldgen(fun,:unary))
        function $fun(expr::String;be=0)
            convert(String, $fun(RExpr(expr);be=be))
        end
    end
end

const MatExpr = Union{Array{Any,2},Array{Expr,2},Array{Symbol,2},Expr,Symbol}

band_matrix(r::Union{Vector,Adjoint,Expr,Symbol},v::Integer) = band_matrix(list(r),v) |> parse
band_matrix(r::T,v::Integer) where T <: Tuple = band_matrix(RExpr(r),v) |> parse
block_matrix(r::Integer,c::Integer,s::VectorAny) = block_matrix(RExpr(r),RExpr(c),list(s)) |> parse
block_matrix(r::Integer,c::Integer,s::T) where T <: Tuple = block_matrix(RExpr(r),RExpr(c),RExpr(s)) |> parse

char_matrix(r::MatExpr,v::Any) = mat(char_matrix(list(r),RExpr(v)) |> parse,r)
char_poly(r::MatExpr,v::Any) = mat(char_matrix(list(r),RExpr(v)) |> parse,r)

extend(a::MatExpr,r::Integer,c::Integer,s) = mat(extend(RExpr(a),r,c,RExpr(s)) |> parse,a)

hessian(r::Reduce.ExprSymbol,l::T) where T <: VectorAny  = hessian(r,list(l))
hessian(r::Reduce.ExprSymbol,l::T) where T <: Tuple  = hessian(r,RExpr(l))

@doc """
    hessian(expr,var_list::Vector)

Computes the Hessian matrix of `expr` with respect to the variables in `var_list`.

This is an n×n matrix where n is the number of variables and the (i,j)th entry is `df(expr,var_list[i],var_list[j])`.
""" hessian

hilbert(a::Integer,r) = hilbert(RExpr(a),RExpr(r)) |> parse

function mat_jacobian(r::T,v::S) where T <: VectorAny where S <: VectorAny
    mat_jacobian(list(r),list(v)) |> parse
end

function mat_jacobian(r::T,v::S) where T <: Tuple where S <: Tuple
    mat_jacobian(RExpr(r),RExpr(v)) |> parse
end

@doc """
    mat_jacobian(expr_list::Vector,var_list::Vector)

Computes the Jacobian matrix of `expr_list` with respect to `var_list`.

This is a matrix whose (i,j)th entry is `df(expr_list[i],var_list[j])`. The matrix is n×m where n is the number of variables and m is the number of expressions.
""" mat_jacobian

export jacobian

jacobian = mat_jacobian

@doc """
    jordan_block(expr,square_size::Integer)

Computes the square Jordan block matrix `J` of dimension `square_size`.

The entries of `J` are `J[i,i] = expr` for i = 1,...,n, `J[i,i+1] = 1` for i = 1,...,n-1, and all other entries are 0.
""" jordan_block

kronecker_product(a::MatExpr,b::MatExpr) = mat(kronecker(RExpr(a),RExpr(b)) |> parse,a,b)

cholesky(r::Array{T,2}) where T <: Number = cholesky(RExpr(r)) |> parse |> mat
coeff_matrix(r::VectorAny) = coeff_matrix(list(r)) |> parse
coeff_matrix(r::T) where T <: Tuple = coeff_matrix(RExpr(r)) |> parse
diagonal(r::VectorAny) = diagonal(list(r)) |> parse
diagonal(r::T) where T <: Tuple = diagonal(RExpr(r)) |> parse
gram_schmidt(r::Vector{<:Vector}) = gram_schmidt(list(r)) |> parse
gram_schmidt(r::T) where T <: Tuple = gram_schmidt(RExpr(r)) |> parse
gram_schmidt(r::Matrix) = gram_schmidt(list(r)) |> parse
hermitian_tp(r::MatExpr) = mat(hermitian_tp(RExpr(r)) |> parse, r)
symmetricp(r::MatExpr) = mat(symmetricp(RExpr(r)) |> parse, r)
toeplitz(r::VectorAny) = toeplitz(list(r)) |> parse
toeplitz(r::T) where T <: Tuple = toeplitz(RExpr(r)) |> parse
triang_adjoint(r::MatExpr) = mat(triang_adjoint(RExpr(r)) |> parse, r)
vandermonde(r::VectorAny) = vandermonde(list(r)) |> parse
vandermonde(r::T) where T <: Tuple = vandermonde(RExpr(r)) |> parse

function __init__()
    load_package(:linalg)
end

end # module


text:
using PlutoEducation, Test

@testset "PlutoEducation" begin
    @test 1 == 1
end


text:
# find distance between two vectors

using JuMP

function convertFloatToInt64(x,digits=2)
  # given a vec, normalize it, convert to an interger
  x = x/sum(x);
  M = 10^digits;
  x = round(x*10^digits);
  x = vec(map(tmp -> convert(Int64,round(tmp)),x));
  return x,M;
end

function getMinCostDistance(x,y,w,wext,verbose=true)
  # lables
  # x: source or supply
  # y: sink or demand
  n = length(x);
  if(n!=length(y))
    error("two vectors must of same size")
  end

  if((size(w)[1]!=size(w)[2]) | size(w)[1]!=n)
    error("cost matrix must be square")
  end
  # convert the vecs to integers
  x,M = convertFloatToInt64(x,2);
  y,M = convertFloatToInt64(y,2);

  # create exraneous source and sink nodes
  # in the case where supply is not equal to demand
  dif = sum(x)-sum(y);
  if(dif>0)
    # if supply is more, create extraneous demand node
    push!(x,0);
    push!(y,dif);
  else
    # if demand is more, create extraneous supply node
    push!(x,-dif);
    push!(y,0);
  end
  dif = sum(x)-sum(y);
  println("x: ",x)
  println("y: ",y)
  @assert dif == 0;

  N = n+1;
  M = sum(x);
  cost = Array{Int64,2}(zeros(N,N));
  # cost metric for given vec
  cost[1:n,1:n] = w;
  # cost to reach from ext.source to any sink node
  cost[N,1:n] = wext;
  # cost to reach ext.sink to any source node
  cost[1:n,N] = wext;

  mdl = Model();

  # flow is non-negative
  @defVar(mdl, flow[i=1:N, j=1:N] >= 0);
  # objective function:
  # minimize the overall flow from source to sink
  @setObjective(mdl, Min, sum{cost[i,j] * flow[i,j], i=1:N, j=1:N});
  # constraint: net out flow at source node should be equal its specified supply
  @addConstraint(mdl, xyconstr[i=1:1:N], sum{flow[i,j], j=1:N} == x[i]);
    # constraint: net out flow at source node should be equal its specified demand
  @addConstraint(mdl, xyconstr[j = 1:N], sum{flow[i,j], i=1:N} == y[j]);

  if(verbose)
    println("Solving original problem...")
  end
  status = solve(mdl);

  if status == :Optimal
    d = getObjectiveValue(mdl)/M/maximum(cost);
    flow = getValue(flow)
    if(verbose)
      @printf("Optimal!\n");
      @printf("Objective value: %d\n", d);
      @printf("Source:\n")
      println(x);
      @printf("Target:\n")
      println(y);
      @printf("Flow:\n");
      for j = 1:n
        @printf("\t%s", string("d",j));
      end
      @printf("\t%s", string("d.ext"));
      @printf("\n");
      for i = 1:n
        @printf("%s", string("s",i));
        for j = 1:N
          @printf("\t%d", (flow[i,j]));
        end
        @printf("\n");
      end
      @printf("%s", string("s.ext"));
      for j = 1:N
        @printf("\t%d", (flow[N,j]));
      end
      @printf("\n");
    end
  else
    d = [];
    if(verbose)
      @printf("No solution\n");
    end
  end
  return status, d, mdl, flow,cost;
end

x1 = [1 0 0 1 200]
x2 = [1 1 1 1 0]
x3 = [1 0 1 1 1]

w = [0 1 2 3 4;1 0 1 2 3;2 1 0 1 2;3 2 1 0 1;4 3 2 1 0]
wext = [1 1 1 1 1]
println("Source: ",x)
println("Target: ",y)
status, d1, mdl, flow, cost = getMinCostDistance(x1,x2,w,wext,false);
println("dist: ",d1)
status, d2, mdl, flow, cost = getMinCostDistance(x2,x3,w,wext,false);
println("dist: ",d2)
status, d3, mdl, flow, cost = getMinCostDistance(x1,x3,w,wext,false);
println("dist: ",d3)
println("triangle: ",[d1+d2,d3])

#

max_stars_repo_name: harshavardhanc/sunbird-analytics

max_issues_repo_name: harshavardhanc/sunbird-analytics

max_forks_repo_name: harshavardhanc/sunbird-analytics


text:
using HPack
using Test

include("codes.jl")
include("musts.jl")


text:
module LibExpat

import Base: getindex, show

@windows_only const libexpat = "libexpat-1"
@unix_only const libexpat = "libexpat"
include("lX_common_h.jl")
include("lX_defines_h.jl")
include("lX_expat_h.jl")
#include("lX_exports_h.jl")

@c Ptr{XML_LChar} XML_ErrorString (Cint,) libexpat

export ETree, xp_parse, find, xpath, @xpath_str
export ParsedData # deprecated

DEBUG = false

macro DBG_PRINT (s)
    quote
        if (DEBUG)
            println($s);
        end
    end
end

type ETree
    # XML Tag
    name::String
    # Dict of tag attributes as name-value pairs
    attr::Dict{String,String}
    # List of child elements.
    elements::Vector{Union(ETree,String)}
    parent::ETree

    ETree() = ETree("")
    function ETree(name)
        pd=new(
            name,
            Dict{String, String}(),
            Union(ETree,String)[])
        pd.parent=pd
        pd
    end
end
typealias ParsedData ETree

function show(io::IO, pd::ETree)
    print(io,'<',pd.name)
    for (name,value) in pd.attr
        print(io,' ',name,'=','"',replace(value,'"',"""),'"')
    end
    if length(pd.elements) == 0
        print(io,'/','>')
    else
        print(io,'>')
        for ele in pd.elements
            if isa(ele, ETree)
                show(io, ele)
            else
                print(io, replace(ele,'<',"<"))
            end
        end
        print(io,'<','/',pd.name,'>')
    end
end

string_value(pd::ETree) = takebuf_string(string_value(pd,IOBuffer()))
function string_value(pd::ETree, str::IOBuffer)
    for node in pd.elements
        if isa(node, String)
            write(str, node::String)
        elseif isa(node,ETree)
            string_value(node::ETree, str)
        end
    end
    str
end

type XPHandle
  parser::Union(XML_Parser,Nothing)
  pdata::ETree
  in_cdata::Bool

  XPHandle(p) = new(p, ETree(""), false)
end

function xp_make_parser(sep='\0')
    p::XML_Parser = (sep == '\0') ? XML_ParserCreate(C_NULL) : XML_ParserCreateNS(C_NULL, sep);
    if (p == C_NULL) error("XML_ParserCreate failed") end

    xph = XPHandle(p)
    p_xph = pointer_from_objref(xph)
    XML_SetUserData(p, p_xph);

    XML_SetCdataSectionHandler(p, cb_start_cdata, cb_end_cdata)
    XML_SetCharacterDataHandler(p, cb_cdata)
    XML_SetCommentHandler(p, cb_comment)
    XML_SetDefaultHandler(p, cb_default)
    XML_SetDefaultHandlerExpand(p, cb_default_expand)
    XML_SetElementHandler(p, cb_start_element, cb_end_element)
#    XML_SetExternalEntityRefHandler(p, f_ExternaEntity)
    XML_SetNamespaceDeclHandler(p, cb_start_namespace, cb_end_namespace)
#    XML_SetNotationDeclHandler(p, f_NotationDecl)
#    XML_SetNotStandaloneHandler(p, f_NotStandalone)
#    XML_SetProcessingInstructionHandler(p, f_ProcessingInstruction)
#    XML_SetUnparsedEntityDeclHandler(p, f_UnparsedEntityDecl)
#    XML_SetStartDoctypeDeclHandler(p, f_StartDoctypeDecl)

    return xph
end

function xp_geterror (p::Union(XML_Parser,Nothing))
    ec = XML_GetErrorCode(p)

    if ec != 0
        @DBG_PRINT (XML_GetErrorCode(p))
        @DBG_PRINT (bytestring(XML_ErrorString(XML_GetErrorCode(p))))

        return  ( bytestring(XML_ErrorString(XML_GetErrorCode(p))),
                XML_GetCurrentLineNumber(p),
                XML_GetCurrentColumnNumber(p) + 1,
                XML_GetCurrentByteIndex(p) + 1
            )
     else
        return  ( "", 0, 0, 0)
     end

end

function xp_geterror(xph::XPHandle)
    return xp_geterror(xph.parser)
end

function xp_close (xph::XPHandle)
  if (xph.parser != nothing)    XML_ParserFree(xph.parser) end
  xph.parser = nothing
end

function start_cdata (p_xph::Ptr{Void})
    xph = unsafe_pointer_to_objref(p_xph)::XPHandle
    @DBG_PRINT ("Found StartCdata")
    xph.in_cdata = true
    return
end
cb_start_cdata = cfunction(start_cdata, Void, (Ptr{Void},))

function end_cdata (p_xph::Ptr{Void})
    xph = unsafe_pointer_to_objref(p_xph)::XPHandle
    @DBG_PRINT ("Found EndCdata")
    xph.in_cdata = false
    return;
end
cb_end_cdata = cfunction(end_cdata, Void, (Ptr{Void},))

function cdata (p_xph::Ptr{Void}, s::Ptr{Uint8}, len::Cint)
    xph = unsafe_pointer_to_objref(p_xph)::XPHandle

    txt = bytestring(s, int(len))
    push!(xph.pdata.elements, txt)

    @DBG_PRINT ("Found CData : " * txt)
    return;
end
cb_cdata = cfunction(cdata, Void, (Ptr{Void},Ptr{Uint8}, Cint))

function comment (p_xph::Ptr{Void}, data::Ptr{Uint8})
    xph = unsafe_pointer_to_objref(p_xph)::XPHandle
    txt = bytestring(data)
    @DBG_PRINT ("Found comment : " * txt)
    return;
end
cb_comment = cfunction(comment, Void, (Ptr{Void},Ptr{Uint8}))

function default (p_xph::Ptr{Void}, data::Ptr{Uint8}, len::Cint)
    xph = unsafe_pointer_to_objref(p_xph)::XPHandle
    txt = bytestring(data)
    @DBG_PRINT ("Default : " * txt)
    return;
end
cb_default = cfunction(default, Void, (Ptr{Void},Ptr{Uint8}, Cint))

function default_expand (p_xph::Ptr{Void}, data::Ptr{Uint8}, len::Cint)
    xph = unsafe_pointer_to_objref(p_xph)::XPHandle
    txt = bytestring(data)
    @DBG_PRINT ("Default Expand : " * txt)
    return;
end
cb_default_expand = cfunction(default_expand, Void, (Ptr{Void},Ptr{Uint8}, Cint))

function attrs_in_to_dict(attrs_in::Ptr{Ptr{Uint8}})
    attrs = Dict{String,String}()

    if (attrs_in != C_NULL)
        i = 1
        attr = unsafe_load(attrs_in, i)
        while (attr != C_NULL)
            k = bytestring(attr)

            i=i+1
            attr = unsafe_load(attrs_in, i)

            if (attr == C_NULL) error("Attribute does not have a name!") end
            v = bytestring(attr)

            attrs[k] = v

            @DBG_PRINT ("$k, $v in $name")

            i=i+1
            attr = unsafe_load(attrs_in, i)
        end
    end

    return attrs
end

function start_element (p_xph::Ptr{Void}, name::Ptr{Uint8}, attrs_in::Ptr{Ptr{Uint8}})
    xph = unsafe_pointer_to_objref(p_xph)::XPHandle
    name = bytestring(name)
    @DBG_PRINT ("Start Elem name : $name,  current element: $(xph.pdata.name) ")

    new_elem = ETree(name)
    new_elem.parent = xph.pdata

    push!(xph.pdata.elements, new_elem)
    @DBG_PRINT ("Found  $name")

    merge!(new_elem.attr, attrs_in_to_dict(attrs_in))
    xph.pdata = new_elem

    return
end
cb_start_element = cfunction(start_element, Void, (Ptr{Void},Ptr{Uint8}, Ptr{Ptr{Uint8}}))

function end_element (p_xph::Ptr{Void}, name::Ptr{Uint8})
    xph = unsafe_pointer_to_objref(p_xph)::XPHandle
    txt = bytestring(name)
    @DBG_PRINT ("End element: $txt, current element: $(xph.pdata.name) ")

    xph.pdata = xph.pdata.parent

    return;
end
cb_end_element = cfunction(end_element, Void, (Ptr{Void},Ptr{Uint8}))

function start_namespace (p_xph::Ptr{Void}, prefix::Ptr{Uint8}, uri::Ptr{Uint8})
    xph = unsafe_pointer_to_objref(p_xph)::XPHandle
    prefix = bytestring(prefix)
    uri = bytestring(uri)
    @DBG_PRINT ("start namespace prefix : $prefix, uri: $uri")
    return;
end
cb_start_namespace = cfunction(start_namespace, Void, (Ptr{Void},Ptr{Uint8}, Ptr{Uint8}))

function end_namespace (p_xph::Ptr{Void}, prefix::Ptr{Uint8})
    xph = unsafe_pointer_to_objref(p_xph)::XPHandle
    prefix = bytestring(prefix)
    @DBG_PRINT ("end namespace prefix : $prefix")
    return;
end
cb_end_namespace = cfunction(end_namespace, Void, (Ptr{Void},Ptr{Uint8}))

# Unsupported callbacks: External Entity, NotationDecl, Not Stand Alone, Processing, UnparsedEntityDecl, StartDocType
# SetBase and GetBase

function xp_parse(txt::String)
    xph = nothing
    xph = xp_make_parser()

    try
        rc = XML_Parse(xph.parser, txt, length(txt.data), 1)
        if (rc != XML_STATUS_OK) error("Error parsing document : $rc") end

        # The root element will only have a single child element in a well formed XML
        return xph.pdata.elements[1]
    catch e
        stre = string(e)
        (err, line, column, pos) = xp_geterror(xph)
        @DBG_PRINT ("$e, $err, $line, $column, $pos")
        rethrow("$e, $err, $line, $column, $pos")

    finally
        if (xph != nothing) xp_close(xph) end
    end
end

function find{T<:String}(pd::ETree, path::T)
    # What are we looking for?
    what = :node
    attr = ""

    pathext = split(path, "#")
    if (length(pathext)) > 2 error("Invalid path syntax")
    elseif (length(pathext) == 2)
        if (pathext[2] == "string")
            what = :string
        else
            error("Unknown extension : [$(pathext[2])]")
        end
    end

    xp= Array((Symbol,Any),0)
    if path[1] == '/'
        # This will treat the incoming pd as the root of the tree
        push!(xp, (:root,:element))
        pathext[1] = pathext[1][2:end]
        #else - it will start searching the children....
    end

    nodes = split(pathext[1], "/")
    idx = false
    descendant = :child
    for n in nodes
        idx = false
        if length(n) == 0
            if descendant == :descendant
                error("too many / in a row")
            end
            descendant = :descendant
            continue
        end
        # Check to see if it is an index into an array has been requested, else default to 1
        m =  match(r"([\:\w]+)\s*(\[\s*(\d+)\s*\])?\s*(\{\s*(\w+)\s*\})?", n)

        if ((m == nothing) || (length(m.captures) != 5))
            error("Invalid name $n")
        else
            node = m.captures[1]
            push!(xp, (descendant,:element))
            descendant = :child
            push!(xp, (:name, convert(T,node)))

            if m.captures[5] != nothing
                if (n == nodes[end])
                    what = :attr
                end
                attr = m.captures[5]
                push!(xp, (:filter, (:attribute, convert(T,attr))))
            end

            if m.captures[3] != nothing
                push!(xp, (:filter, (:number,int(m.captures[3]))))
                idx = true
            end
        end
    end

    pd = xpath(pd, XPath{T,Vector{ETree}}((:xpath, xp)))
    if what == :node
        if idx
            if length(pd) == 1
                return pd[1]
            else
                return nothing
            end
        else
            # If caller did not specify an index, return a list of leaf nodes.
            return pd
        end
    elseif length(pd) == 0
        return nothing

    elseif length(pd) != 1
        error("More than one instance of $pd, please specify an index")

    else
        pd = pd[1]
        if what == :string
            return string_value(pd)

        elseif what == :attr
            return pd.attr[attr]

        else
            error("Unknown request type")
        end
    end

    return nothing
end

include("xpath.jl")

include("streaming.jl")
export XPCallbacks, XPStreamHandler
export parse, stop, pause, resume, free, parsefile

end


text:
using Plots
using Random
using Distributions
using LaTeXStrings

global e = Base.MathConstants.e

function estimation(counter, n, m, R)
    est = zero(BigInt)
    for i in 1:m
        est += counter(n, i, R)
    end
    est
end

function sd1(n, m, R)
    est = zero(BigInt)
    for r in 1:m-1
        est += 2*(e*n/r)^r*log(R)
    end
    est += (e*n/m)^m
    return est
end

function sd2(n, m, R)
    return big((n/m))^m * big((n/(n-m)))^(n-m) * log(e*n*R)
end

function sd3(n, m, R)
    est = zero(BigInt)
    for r in 1:m-1
        est += binomial(BigInt(n), floor(Int,r))*log(R)
    end
    est += binomial(BigInt(n), floor(Int,m))
    return est
end

function sd4(n, m, R)
    est = zero(BigInt)
    for r in 1:m-1
        for s in 1:r
            est += binomial(BigInt(n), floor(Int,s))*log(R)
        end
    end
    for s in 1:m-1
        est += binomial(BigInt(n), floor(Int,s))*log(R)
    end
    est += binomial(BigInt(n), floor(Int,m))
end

Plots.font("DejaVu")
Plots.PyPlotBackend()
function fullestimation(n, maxn)
    R = n^(5)

    yticks = ([n, binomial(BigInt(n), BigInt(sqrt(n))), binomial(BigInt(n), BigInt(n/2)) ],
    [L"n", L"\binom{\!\!n\!\!}{\!\!\sqrt{n}\!\!}" , L"\binom{\!\!n\!\!}{\!\!n/2\!\!}"])
    sd1data = [estimation(sd1, n, m, R) for m in 1:maxn]
    sd3data = [estimation(sd3, n, m, R) for m in 1:maxn]
    sd4data = [sd4(n, m, R) for m in 1:maxn]
    plt = plot(sd1data, yaxis=:log10, label=L"$T_1$  (SD-(1+1)EA)", size=(320,320), linestyle=:solid, dpi=500,
    legend=:bottomright, yticks=yticks, ylims=(1, Inf),
    xticks=([log(n), sqrt(n), n/4,n/2], [L"\ln n",L"\sqrt{n}", L"n/4" ,L"n/2"]),
    # xticks=([log(n), 1.5sqrt(n), n/4,n/2], [L"\ln n",L"\sqrt{n}", L"n/4" ,L"n/2"]),
    xlabel=L"m", ylabel=L"T_i(m)")
    plt = plot!(sd3data,  label=L"$T_2$  (SD-RLS$^p$)", linestyle=:dashdot)
    plt = plot!(sd4data, label=L"$T_3$  (SD-RLS$^r$)", linestyle=:dash)
    plt
end

plt = fullestimation(500, 250)
savefig(plt, "./estimation-large.png")

plt = fullestimation(500, 25)
savefig(plt, "./estimation-small.png")

max_stars_repo_name: DTUComputeTONIA/SDRLS

max_issues_repo_name: DTUComputeTONIA/SDRLS

max_forks_repo_name: DTUComputeTONIA/SDRLS


text:
## until these are defined in the Data.Source interface, use the reset! method
## to get non-random-access data sources ready to stream.
reset!(x) = x
reset!(x::TupleIterator) = reset!(x.source)

"""
    get_unique!(iter::TupleIterator, cols::Vector{Int})

Update the schema stored on `iter` with the unique values for the specified
columns.  The unique values are stored in a Dict in the metadata field of the
underlying schema.  `iter` is reset before and after this operation.

"""
function get_unique!{T<:Integer}(iter::TupleIterator, cols::Vector{T})
    isempty(cols) && return(iter)
    reset!(iter)
    types = Data.types(iter)
    uniqs = [Vector{types[i]}() for i in cols]
    seens = [Set{types[i]}() for i in cols]

    for t in iter
        for (i,col) in enumerate(cols)
            x = t[col]
            if !(x in seens[i])
                push!(seens[i], x)
                push!(uniqs[i], x)
            end
        end
    end

    sch = Data.schema(iter)

    uniq_dict = get!(sch.metadata, :unique, Dict())
    header = Data.header(sch)
    for (col, uniq) in zip(cols, uniqs)
        uniq_dict[header[col]] = uniq
    end

    reset!(iter)
    iter
end

get_unique!(iter::TupleIterator, cols::Vector) =
    (sch = Data.schema(iter); get_unique!(iter, [sch[string(c)] for c in cols]))


text:
"""
    CUDA.CUDNN

High level interface to cuDNN functions. See
https://github.com/JuliaGPU/CUDA.jl/blob/master/lib/cudnn/README.md
for a design overview.
"""
module CUDNN

using ..APIUtils

using ..CUDA
using ..CUDA: CUstream, libraryPropertyType
using ..CUDA: libcudnn, @retry_reclaim, isdebug, @context!

using CEnum

using Memoize

using DataStructures

# core library
include("libcudnn_common.jl")
include("error.jl")
include("libcudnn.jl")
include("libcudnn_deprecated.jl")

# low-level wrappers
include("util.jl")
include("base.jl")
include("descriptors.jl")
include("tensor.jl")
include("inplace.jl")
include("optensor.jl")
include("reduce.jl")
include("convolution.jl")
include("pooling.jl")
include("activation.jl")
include("softmax.jl")
include("dropout.jl")
include("rnn.jl")
include("multiheadattn.jl")
include("normalization.jl")

function math_mode(mode=CUDA.math_mode())
    if mode == CUDA.PEDANTIC_MATH
        # don't use tensor cores.
        # on A100, only use them for TF32
        CUDNN_DEFAULT_MATH
    elseif mode == CUDA.DEFAULT_MATH
        # allow tensor core usage
        CUDNN_TENSOR_OP_MATH
    elseif mode == CUDA.FAST_MATH
        # also downcast inputs
        CUDNN_TENSOR_OP_MATH_ALLOW_CONVERSION
    end
end

# cache for created, but unused handles
const idle_handles = HandleCache{CuContext,cudnnHandle_t}()

function handle()
    state = CUDA.active_state()
    handle, stream = get!(task_local_storage(), (:CUDNN, state.context)) do
        new_handle = pop!(idle_handles, state.context) do
            cudnnCreate()
        end

        finalizer(current_task()) do task
            push!(idle_handles, state.context, new_handle) do
                @context! skip_destroyed=true state.context cudnnDestroy(new_handle)
            end
        end

        cudnnSetStream(new_handle, state.stream)

        new_handle, state.stream
    end::Tuple{cudnnHandle_t,CuStream}

    if stream != state.stream
        cudnnSetStream(handle, state.stream)
        task_local_storage((:CUDNN, state.context), (handle, state.stream))
    end

    return handle
end

function log_message(sev, udata, dbg_ptr, ptr)
    # "Each line of this message is terminated by \0, and the end of the message is
    # terminated by \0\0"
    len = 0
    while true
        if unsafe_load(ptr, len+1) == '\0' && unsafe_load(ptr, len+2) == '\0'
            break
        end
        len += 1
    end
    str = unsafe_string(ptr, len)
    lines = split(str, '\0')
    msg = join(str, '\n')

    # TODO: inspect `sev` to generate an appropriate message (@debug, @info, etc)
    dbg = unsafe_load(dbg_ptr)
    println(msg)
    return
end

function __runtime_init__()
    # enable library logging when launched with JULIA_DEBUG=CUDNN
    # FIXME: this doesn't work, and the mask remains 0 (as observed with cudnnGetCallback)
    if isdebug(:init, CUDNN)
        callback = @cfunction(log_message, Nothing,
                              (cudnnSeverity_t, Ptr{Cvoid}, Ptr{cudnnDebug_t}, Ptr{UInt8}))
        cudnnSetCallback(typemax(UInt32), C_NULL, callback)
    end
end

end


text:
"""
    M_file(P, col::int64)

    Generate filepath and extract indices from spreadsheet log.xlsx

    Arguments
    - 'col::Int64' : column of file of interest in accompanying spread sheet
"""
function M_file(col::Int64)
    if Sys.iswindows()
        xf = XLSX.readxlsx("C:\\Users\\cwaha\\Dropbox\\My PC (DESKTOP-8JF2H49)\\Documents\\UCL\\raw_lab data\\log.xlsx")
    else
        xf = XLSX.readxlsx("/Users/christopherharbord/Dropbox/My PC (DESKTOP-8JF2H49)/Documents/UCL/raw_lab data/log.xlsx")
    end
    sh = xf[XLSX.sheetnames(xf)[1]]
    info = sh[:]
    exp_info = Dict()
    exp_info[:I]=[info[col,4],info[col,5],info[col,6]] # Indice at HP
    exp_info[:P_exp] = info[col,2]
    exp_info[:T] = info[col,7] # Temperature of experiment
    exp_info[:εr] = info[col,8] # Strain rate of test
    exp_info[:K_mm_kN] = info[col,9] # Machine compliance [mm kN^-1]
    exp_info[:L_mm] = info[col,10] # Sample length
    exp_info[:d_mm] = info[col,11] # Diameter of sample [m]
    fil = info[col,1]
    if Sys.iswindows()
        fid = "C:\\Users\\cwaha\\Dropbox\\My PC (DESKTOP-8JF2H49)\\Documents\\UCL\\raw_lab data\\"*fil*".tdms"
    else
        fid = "/Users/christopherharbord/Dropbox/My PC (DESKTOP-8JF2H49)/Documents/UCL/raw_lab data/"*fil*".tdms"
    end
    return fid, exp_info
end
"""
    M_read(P)

    Open and read contents of .tdms file
    #Arguments
    - fid : file path string
    return P
"""
function M_read(fid::String)
    tdmsIN = readtdms(fid)
    P = Dict()
    P[:t_s] = tdmsIN.groups["Numeric"]["TimeStamp"].data
    P[:F_kN] = tdmsIN.groups["Numeric"]["Load"].data
    P[:Ua_mm] = tdmsIN.groups["Numeric"]["Displacement"].data
    P[:U1_mm] = tdmsIN.groups["Numeric"]["LVDT 1"].data
    P[:U2_mm] = tdmsIN.groups["Numeric"]["LVDT 2"].data
    P[:Pc1_MPa] = tdmsIN.groups["Numeric"]["Pc 700MPa"].data
    P[:Pc2_MPa] = tdmsIN.groups["Numeric"]["PC 1400 MPa"].data
    P[:Pf_MPa] = tdmsIN.groups["Numeric"]["Pore Pressure"].data
    P[:PpVol_mm3] = tdmsIN.groups["Numeric"]["PP vol"].data
    return P
end

function t_conv(P)
    P[:t_date_time] = unix2datetime(datetime2unix(DateTime(1904,1,1,0,0,0)) .+Int.(Round.(P[:t_s]*1000)))
end

"""
    M_reduce!(P, exp_info)

    Reduce Murrell data
    #Arguments
    * P : dictionary containing raw data from the Murrell, must include indices of hit point and sample information
    * exp_info : dictionary containing experimental parameters
"""

function M_reduce!(P,exp_info)
    I1 = exp_info[:I][1]
    P[:t_s_c] = P[:t_s].-P[:t_s][I1]
    P[:F_kN_c] = P[:F_kN] .-P[:F_kN][I1]
    P[:U_mm_c] = ((P[:U1_mm].+P[:U2_mm]).-(P[:U1_mm][I1]+P[:U2_mm][I1]))./2
    P[:U_mm_fc] = P[:U_mm_c] .-(P[:F_kN_c]*exp_info[:K_mm_kN])
    P[:ε] = P[:U_mm_fc]./exp_info[:L_mm]
    P[:Jr] = JR!(P, exp_info)
    P[:F_kN_j] = P[:F_kN_c] .-P[:Jr]
    P[:σ_MPa] = P[:F_kN_c]./(0.25e-6π*exp_info[:d_mm]^2) .*1e-3
    P[:σ_MPa_j] = P[:F_kN_j]./(0.25e-6π*exp_info[:d_mm]^2) .*1e-3
end

"""
    H_mod(P,ds::Int64)

    Calculate tangent moduli and yield stress using robust numerical differentiation

    #Arguments
    - 'P' : dictionary information pertaining to the 'Murrell'
"""
function Youngs_mod!(P, exp_info)
    ε1 = Array{Float64,1}(undef,1)
    σ1 = Array{Float64,1}(undef,1)
    t1 = Array{Float64,1}(undef,1)
    I1 = exp_info[:I][1]
    I2 = exp_info[:I][2]
    ε = P[:ε][I1:I2]
    σ = P[:σ_MPa][I1:I2]
    t = P[:t_s][I1:I2]
    # Youngs modulus calculation
    while true
        fig = figure()
        plot(ε,σ)
        xlim(0,0.005)
        ylim(0,300)
        g_data = ginput(2)
        close(fig)
        IE = (findfirst(ε.>g_data[1][1]), findfirst(ε.>g_data[2][1]))
        P[:EE] = linfit(ε[IE[1]:IE[2]], σ[IE[1]:IE[2]])
        P[:E_GPa] = P[:EE][1]/1e3
        fig = figure()
        plot(ε[IE[1]:IE[2]],σ[IE[1]:IE[2]],ε[IE[1]:IE[2]],M(ε[IE[1]:IE[2]],P[:EE].param))
        title("Youngs' modulus ="*string(P[:E_GPa])*"GPa")
        println("Fit ok? (y/n)")
        if readline() == "y"
            close(fig)
            break
        else
            close(fig)
        end
    end
end

function YS!(P,exp_info)
    I1 = exp_info[:I][1]
    I2 = exp_info[:I][2]
    ε1 = Array{Float64,1}(undef,1)
    σ1 = Array{Float64,1}(undef,1)
    t1 = Array{Float64,1}(undef,1)
    ε = P[:ε][I1:I2-1000]
    σ = P[:σ_MPa][I1:I2-1000]
    t = P[:t_s][I1:I2-1000]
    while true
        println("Enter moving average:")
        ds = parse(Float64,readline())
        for i in 1:Int(floor(length(P[:ε][I1:I2-1000])/ds))-1
            push!(ε1,sum(ε[Int(i*ds-ds+1):Int(i*ds)])/ds)
            push!(σ1,sum(σ[Int(i*ds-ds+1):Int(i*ds)])/ds)
            push!(t1,t[Int(i*ds)+Int(ceil(ds/2))])
        end
        deleteat!(ε1,1)
        deleteat!(σ1,1)
        deleteat!(t1,1)
        dσ = TVRegDiff(σ1,200, 0.2, dx=1, ε=1e-6)
        dε = TVRegDiff(ε1,200, 0.2, dx=1, ε=1e-6)
        P[:h] = dσ./dε*1e-3 # Calculate the tangent modulus in GPa, excluding initial non-linearity
        P[:ε_h] = ε1
        h_max = maximum(P[:h])
        II = findfirst(P[:h] .== h_max)
        P[:H_GPa] = sum(P[:h][P[:h].<0.4P[:E_GPa]])/sum(P[:h].<0.4P[:E_GPa])
        P[:ε_E] = ε1[II]
        P[:σ_E] = σ1[II]
        P[:YSa_MPa] = σ1[findfirst((P[:h].<0.8P[:E_GPa]) .& (ε1 .> ε1[II]))] #80% yield
        P[:εa] = ε1[findfirst((P[:h].<0.8P[:E_GPa]) .& (ε1 .> ε1[II]))]
        P[:YSb_MPa] = σ1[findfirst((P[:h].<0.5P[:E_GPa]) .& (ε1 .> ε1[II]))] #50% yield
        P[:εb] = ε1[findfirst((P[:h].<0.5P[:E_GPa]) .& (ε1 .> ε1[II]))]
        M(ε, p) = p[1]*ε .+p[2]
        fig = figure()
        ax=[]
        push!(ax,subplot(2,1,1))
        ax[1].plot(ε,σ,"0.5")
        ax[1].plot(ε,M(ε,P[:EE].param),"--k")
        ax[1].scatter(P[:εa],P[:YSa_MPa],s=100)
        ax[1].scatter(P[:εb],P[:YSb_MPa],s=100)
        ax[1].set_xlim(0,P[:εb]+0.001)
        ax[1].set_ylim(0,P[:YSb_MPa]+20)
        push!(ax,subplot(2,1,2))
        ax[2].plot(P[:ε_h],P[:h])
        ax[2].plot(P[:ε_h],M(P[:ε_h],(0,P[:E_GPa])),"--k")
        ax[2].plot(P[:ε_h],M(P[:ε_h],(0,0.8P[:E_GPa])),"--k")
        ax[2].plot(P[:ε_h],M(P[:ε_h],(0,0.5P[:E_GPa])),"--k")
        ax[2].scatter(P[:εa],0.8P[:E_GPa])
        ax[2].scatter(P[:εb],0.5P[:E_GPa])
        ax[2].set_xlim(0,P[:εb]+0.001)
        ax[2].set_ylim(0,2P[:E_GPa])
        println("Fit ok (y/n)")
        if readline() == "y"
            close(fig)
            break
        else
            close(fig)
        end
    end
end

"""
    SFP_read()

    Read and return data from program SFP_furnace_control_V2.vi

    #Arguments
    - 'SFP' : empty dictionary to store extracted data
    - 'fil::string' : path to file
"""
function SFP_read(fid)
    path="/Users/christopherharbord/Dropbox/My PC (DESKTOP-8JF2H49)/Documents/UCL/Furnace_calibration/SFP_logging/"*fil
    dat = readdlm(path)
    headers = dat[:,1]
    I = findall(x->x==headers[1],dat[:,1])
    SFP = Dict()
    SFP[:t] = dat[:,1]
    SFP[:PF] = dat[:,3]
    SFP[:S1] = dat[:,4]
    SFP[:S2] = dat[:,5]
    SFP[:S3] = dat[:,6]
    SFP[:TC1] = dat[:,7]
    SFP[:TC2] = dat[:,8]
    close(dat)
    deleteat!(SFP[:t],I)
    deleteat!(SFP[:PF],I)
    deleteat!(SFP[:S1],I)
    deleteat!(SFP[:S2],I)
    deleteat!(SFP[:S3],I)
    deletat!(SFP[:TC1],I)
    deletat!(SFP[:TC2],I)
end
"""
    JR(P, exp_info)

    Calculate strength contribution of copper jacket for corrections

    #Arguments
    * 'P' : dictionary containing processed data from the 'Murrell'
    * exp_info : details of experiment
"""
function JR!(P, exp_info)
    EaJ = 197000 # Copper activation enthalpy
    n1J = 4.8 # Empirical factor 1
    n2J = 0.22 # Empirical factor 2
    ε̇j = P[:εr]*(exp_info[:d_mm]/exp_info[:L_mm]) # Jacket strain rate
    Jr = 2*10^((log10(ε̇j*exp(EaJ/(8.3145*(exp_info[:T]+278)))))/n1J-n2J) # Copper flow stress at experiment conditions
    Ja = π*exp_info[:d_mm]*exp_info[:L_mm]*1e-6 # Jacket area
    JR = Jr.*Ja.*(P[:ε]*(exp_info[:d_mm]/exp_info[:L_mm])).*1e-3 # Force due to jacket assuming linear increase due to incremental strain
end

function M_interp!(P, t_UT)
    P[:F_kN_i] = lininterp(P[:t_s],P[:F_kN_j], t_UT)
    P[:U_mm_i] = lininterp(P[:t_s],P[:U_mm_c], t_UT)
    P[:σ_MPa_i] = lininterp(P[:t_s],P[:σ_MPa_j], t_UT)
    P[:σ3_MPa_i] = lininterp(P[:t_s],P[:PC2_MPa], t_UT)
    P[:ε_i] = lininterp(P[:t_s],P[:ε], t_UT)
end

max_stars_repo_name: chrisharb/Murrell

max_issues_repo_name: chrisharb/Murrell

max_forks_repo_name: chrisharb/Murrell


text:
# All code is attached to its underlying database source
struct SourceCode{Source}
    source::Source
    code::Expr
end

"""
    struct BySQL{Source}

If you would like a statement to be evaluated by SQL, not Julia, and
none of the arguments are SQL code, you can use BySQL to hack dispatch.

```jldoctest
julia> using QuerySQLite

julia> using Query: @map

julia> using DataValues: DataValue

julia> database = Database(joinpath(pathof(QuerySQLite) |> dirname |> dirname, "test", "Chinook_Sqlite.sqlite"));

julia> result = database.Track |> @map({a = rand(BySQL(_), Int)});

julia> collect(result)[1].a isa DataValue{Int}
true
```
"""
struct BySQL{Source}
    source::Source
end

export BySQL

# Every time `SourceCode` objects are combined, check to see whether they all come from the same source
function pop_source!(sources, something)
    something
end
function pop_source!(sources, source_code::SourceCode)
    push!(sources, source_code.source)
    source_code.code
end
function pop_source!(sources, by_sql::BySQL)
    push!(sources, by_sql.source)
    by_sql
end
function key_pop_source!(sources, (key, source_code))
    code = pop_source!(sources, source_code)
    Expr(:kw, key, code)
end

function combine_sources(a_function, source_codes...; key_source_codes...)
    sources = Set(Any[])
    codes = [
        pop_source!(sources, source_code)
        for source_code in source_codes
    ]
    filter!(codes) do code
        !(code isa BySQL)
    end
    key_codes = [
        key_pop_source!(sources, key_source_code)
        for key_source_code in key_source_codes
    ]
    if length(sources) != 1
        error("Expected exactly one source; got ($(sources...))")
    else
        SourceCode(
            first(sources),
            Expr(:call, a_function, Expr(:parameters, key_codes...), codes...)
        )
    end
end

function numbered_argument(number)
    Symbol(string("argument", number))
end

function assert_type(argument, type)
    Expr(:(::), argument, type)
end

# Splat Varargs
function maybe_splat(argument, a_type)
    if @capture a_type Vararg{AType_}
        Expr(:(...), argument)
    else
        argument
    end
end

# `@code_instead` hijacks call to create Julia expressions instead of evaluating functions
function code_instead(location, a_function, types...)
    arguments = ntuple(numbered_argument, length(types))
    keywords = Expr(:parameters, Expr(:..., :keywords))
    Expr(:function,
        Expr(:call, a_function, keywords, map(assert_type, arguments, types)...),
        Expr(:block, location, Expr(:call,
            combine_sources,
            keywords,
            a_function,
            map(maybe_splat, arguments, types)...
        ))
    )
end

macro code_instead(a_function, types...)
    code_instead(__source__, a_function, types...) |> esc
end

# Currently, query doesn't do anything
function QueryOperators.query(source_code::SourceCode)
    source_code
end

@code_instead (==) SourceCode Any
@code_instead (==) Any SourceCode
@code_instead (==) SourceCode SourceCode

@code_instead (!=) SourceCode Any
@code_instead (!=) Any SourceCode
@code_instead (!=) SourceCode SourceCode

@code_instead (!) SourceCode

@code_instead (&) SourceCode Any
@code_instead (&) Any SourceCode
@code_instead (&) SourceCode SourceCode

@code_instead (|) SourceCode Any
@code_instead (|) Any SourceCode
@code_instead (|) SourceCode SourceCode

@code_instead (*) SourceCode Any
@code_instead (*) Any SourceCode
@code_instead (*) SourceCode SourceCode

@code_instead (-) SourceCode Any
@code_instead (-) Any SourceCode
@code_instead (-) SourceCode SourceCode

@code_instead (+) SourceCode Any
@code_instead (+) Any SourceCode
@code_instead (+) SourceCode SourceCode

@code_instead (/) SourceCode Any
@code_instead (/) Any SourceCode
@code_instead (/) SourceCode SourceCode

@code_instead (%) SourceCode Any
@code_instead (%) Any SourceCode
@code_instead (%) SourceCode SourceCode

@code_instead abs SourceCode

# TODO: add more methods
@code_instead char SourceCode Vararg{Any}

@code_instead coalesce SourceCode Vararg{Any}

@code_instead convert Type{Int} SourceCode

@code_instead QueryOperators.drop SourceCode Integer

@code_instead QueryOperators.filter SourceCode Any Expr

@code_instead format SourceCode AbstractString
@code_instead format AbstractString SourceCode
@code_instead format SourceCode SourceCode

@code_instead QueryOperators.groupby SourceCode Any Expr Any Expr

@code_instead hex SourceCode

@code_instead if_else SourceCode Any Any
@code_instead if_else Any SourceCode Any
@code_instead if_else Any Any SourceCode
@code_instead if_else Any SourceCode SourceCode
@code_instead if_else SourceCode Any SourceCode
@code_instead if_else SourceCode SourceCode Any
@code_instead if_else SourceCode SourceCode SourceCode

# TODO: add more methods
@code_instead in SourceCode Any

@code_instead instr SourceCode Any
@code_instead instr Any SourceCode
@code_instead instr SourceCode SourceCode

@code_instead isequal SourceCode Any
@code_instead isequal Any SourceCode
@code_instead isequal SourceCode SourceCode

@code_instead isless SourceCode Any
@code_instead isless Any SourceCode
@code_instead isless SourceCode SourceCode

@code_instead ismissing SourceCode

@code_instead join SourceCode

@code_instead QueryOperators.join SourceCode SourceCode Any Expr Any Expr Any Expr

@code_instead max SourceCode Vararg{Any}
@code_instead maximum SourceCode

@code_instead mean SourceCode

@code_instead min SourceCode Vararg{Any}
@code_instead minimum SourceCode

@code_instead length SourceCode

@code_instead lowercase SourceCode

@code_instead QueryOperators.map SourceCode Any Expr

@code_instead occursin AbstractString SourceCode
@code_instead occursin SourceCode AbstractString
@code_instead occursin SourceCode SourceCode

@code_instead QueryOperators.orderby SourceCode Any Expr

@code_instead QueryOperators.orderby_descending SourceCode Any Expr

@code_instead rand BySQL Type{Int}

@code_instead randstring BySQL Integer
@code_instead randstring SourceCode

@code_instead replace SourceCode Pair

@code_instead repr SourceCode

@code_instead round SourceCode

@code_instead secondary SourceCode

# TODO: add more methods
@code_instead string SourceCode Vararg{Any}

@code_instead strip SourceCode
@code_instead strip SourceCode Char

@code_instead SubString SourceCode Number Number
@code_instead SubString SourceCode Number

@code_instead sum SourceCode

@code_instead QueryOperators.take SourceCode Any

@code_instead QueryOperators.thenby SourceCode Any Expr

@code_instead QueryOperators.thenby_descending SourceCode Any Expr

@code_instead type_of SourceCode

@code_instead QueryOperators.unique SourceCode Any Expr

@code_instead uppercase SourceCode

# TODO: add
# printf
# zeroblob


text:
for (name,T) in [("byte", UInt8),("half",UInt16),("word",UInt32)]
    @eval @inline $(Symbol("load_$name"))(ctx, off::Int64) =
        unsafe_trunc($T,
                     ccall($("llvm.bpf.load.$name"), llvmcall, UInt64,
                           (Core.LLVMPtr{$T, 0}, Int64),
                           reinterpret(Core.LLVMPtr{$T, 0}, ctx), off))
    @eval @inline $(Symbol("load_$name"))(ctx, off::Union{Int8,Int16,Int32}) =
        $(Symbol("load_$name"))(ctx, Core.sext_int(Int64, off))
    @eval @inline $(Symbol("load_$name"))(ctx, off::Union{UInt8,UInt16,UInt32}) =
        $(Symbol("load_$name"))(ctx, Core.zext_int(Int64, off))
    # TODO: We should early exit with typemax(UInt64) if off > typemax(Int64)
    @eval @inline $(Symbol("load_$name"))(ctx, off::UInt64) =
        $(Symbol("load_$name"))(ctx, reinterpret(Int64, off))
end


text:
using Base.Filesystem

function parse_output(output)
    lines = split(output, '\n')
    return filter(line-> occursin("error", line) || occursin("error", lowercase.(line)), lines);
end

function main()
    filenames = map(x->"error" * string(x) * ".c", 1:15)
    filenames = map(x->joinpath("test_files", x), filenames)

    for name in filenames
        println("-- compiling: ", name)
        output = read(`./project3_16 $name`, String)
        parsed = parse_output(output)

        for i in parsed
            println(parsed)
        end
        print('\n')
    end
end

main()

max_stars_repo_name: Red-Portal/compiler-assignment

max_issues_repo_name: Red-Portal/compiler-assignment

max_forks_repo_name: Red-Portal/compiler-assignment


text:
module TestHandlers

using HTTP
using Test
import Base.==

==(a::HTTP.Response,b::HTTP.Response) = (a.status  == b.status)    &&
                                        (a.version == b.version)   &&
                                        (a.headers == b.headers)   &&
                                        (a.body    == b.body)

@testset "Handlers" begin
    @testset "FooRouter" begin
        r = HTTP.Router()
        f = HTTP.Handlers.RequestHandlerFunction((req) -> HTTP.Response(200))
        HTTP.@register(r, "/test", f)
        req = HTTP.Request()
        req.target = "/test"
        @test HTTP.handle(r, req) == HTTP.Response(200)
    end

    @testset "HTTP Router" begin
        f = HTTP.Handlers.RequestHandlerFunction((req) -> HTTP.Response(200))
        @test HTTP.handle(f, HTTP.Request()) == HTTP.Response(200)

        r = HTTP.Router()
        @test isempty(r.routes)
        @test HTTP.handle(r, HTTP.Request()) == HTTP.Response(404)

        HTTP.@register(r, "/path/to/greatness", f)
        req = HTTP.Request()
        req.target = "/path/to/greatness"
        @test HTTP.handle(r, req) == HTTP.Response(200)

        f2 = HTTP.Handlers.RequestHandlerFunction((req) -> HTTP.Response(201))
        HTTP.@register(r, "/next/path/to/greatness", f2)
        req = HTTP.Request()
        req.target = "/next/path/to/greatness"
        @test HTTP.handle(r, req) == HTTP.Response(201)

        r = HTTP.Router()
        HTTP.@register(r, "GET", "/sget", f)
        HTTP.@register(r, "POST", "/spost", f)
        HTTP.@register(r, "POST", "/tpost", f)
        req = HTTP.Request("GET", "/sget")
        @test HTTP.handle(r, req) == HTTP.Response(200)
        req = HTTP.Request("POST", "/sget")
        @test HTTP.handle(r, req) == HTTP.Response(404)
        req = HTTP.Request("GET", "/spost")
        @test HTTP.handle(r, req) == HTTP.Response(404)
        req = HTTP.Request("POST", "/spost")
        @test HTTP.handle(r, req) == HTTP.Response(200)
        req = HTTP.Request("GET", "/tpost")
        @test HTTP.handle(r, req) == HTTP.Response(404)
        req = HTTP.Request("POST", "/tpost")
        @test HTTP.handle(r, req) == HTTP.Response(200)

        r = HTTP.Router()
        HTTP.@register(r, "/test", f)
        HTTP.@register(r, "/test/*", f2)
        f3 = HTTP.Handlers.RequestHandlerFunction((req) -> HTTP.Response(202))
        HTTP.@register(r, "/test/sarv/ghotra", f3)
        f4 = HTTP.Handlers.RequestHandlerFunction((req) -> HTTP.Response(203))
        HTTP.@register(r, "/test/*/ghotra/seven", f4)

        req = HTTP.Request()
        req.target = "/test"
        @test HTTP.handle(r, req) == HTTP.Response(200)

        req.target = "/test/sarv"
        @test HTTP.handle(r, req) == HTTP.Response(201)

        req.target = "/test/sarv/ghotra"
        @test HTTP.handle(r, req) == HTTP.Response(202)

        req.target = "/test/sar/ghotra/seven"
        @test HTTP.handle(r, req) == HTTP.Response(203)
    end

end

end # module


text:
# Mld.jl

###### mld #####
"""
    mld(v)

Compute the Mean log deviation of a vector `v`.

# Examples
```julia
julia> using Inequality
julia> mld([8, 5, 1, 3, 5, 6, 7, 6, 3])
0.1397460530936332
```
"""
mld(v::AbstractVector{<:Real})::Float64 = log(Statistics.mean(v[v .!= 0])) - Statistics.mean(log.(v[v .!= 0]))

###### weighted mld #####
"""
    mld(v, w)

Compute the weighted Mean log deviation of a vector `v` using weights given by a weight vector `w`.

Weights must not be negative, missing or NaN. The weights and data vectors must have the same length.

# Examples
```julia
julia> using Inequality
julia> mld([8, 5, 1, 3, 5, 6, 7, 6, 3], collect(0.1:0.1:0.9))
0.10375545537468206
```
"""
function mld(v::AbstractVector{<:Real}, w::AbstractVector{<:Real})::Float64
    checks_weights(v, w)
    w = w[v .!= 0]
    v = v[v .!= 0]

    return -sum(w .* log.(v/mean(v, weights(w))) ) / sum(w)
end

function mld(v::AbstractVector{<:Real}, w::AbstractWeights)::Float64
    checks_weights(v, w)
    w = w[v .!= 0]
    v = v[v .!= 0]

    return -sum(w .* log.(v/mean(v, weights(w))))/ sum(w)
end

"""
    wmld(v, w)

Compute the Mean log deviationof `v` with weights `w`. See also [`mld`](@mld)
"""
wmld(v::AbstractVector{<:Real}, w::AbstractVector{<:Real}) = mld(v, w)

wmld(v::AbstractVector{<:Real}, w::AbstractWeights) = mld(v, w)


text:
using BioSequences.FASTA, DataFrames, CSV, nucpos, Statistics

sib_pos = "/bench/PhD/danpos_results/pooled/sib.Fnor.smooth.positions.xls"
rys_pos = "/bench/PhD/danpos_results/pooled/rys.Fnor.smooth.positions.xls"
sib_refpos = "/bench/PhD/danpos_results/pooled/sib.Fnor.smooth.positions.ref_adjust.xls"
rys_refpos = "/bench/PhD/danpos_results/pooled/rys.Fnor.smooth.positions.ref_adjust.xls"

sib_df=CSV.read(sib_pos)
rys_df=CSV.read(rys_pos)
sib_refdf=CSV.read(sib_refpos)
rys_refdf=CSV.read(rys_refpos)

Threads.@threads for t in 1:4
    t==1 && nucpos.map_positions!(sib_df, rys_df)
    t==2 && nucpos.map_positions!(rys_df, sib_df)
    t==3 && nucpos.map_positions!(sib_refdf, rys_refdf)
    t==4 && nucpos.map_positions!(rys_refdf, sib_refdf)
end

println("Unmapped positions: sib $(length(findall(iszero,sib_df.mapped_pos))), sib_ref $(length(findall(iszero,sib_refdf.mapped_pos)))")
println("Unmapped positions: rys $(length(findall(iszero,rys_df.mapped_pos))), rys_ref $(length(findall(iszero,rys_refdf.mapped_pos)))")
println("Avg overlap: sib $(mean(sib_df.rel_overlap)), sib_ref $(mean(sib_refdf.rel_overlap))")
println("Avg overlap: rys $(mean(rys_df.rel_overlap)), rys_ref $(mean(rys_refdf.rel_overlap))")


text:
export interval_overlaps

function intervals(dataset::HONData)
    simplices, nverts, times = dataset.simplices, dataset.nverts, dataset.times
    I, J, T = Int64[], Int64[], Int64[]
    curr_ind = 1
    for (nv, time) in zip(nverts, times)
        simplex = simplices[curr_ind:(curr_ind + nv - 1)]
        curr_ind += nv
        for (i, j) in combinations(simplex, 2)
            push!(I, i)
            push!(J, j)
            push!(T, time)
        end
    end
    sp = sortperm(I)
    I, J, T = I[sp], J[sp], T[sp]

    minmax_I = Int64[]
    minmax_J = Int64[]
    min_V = Int64[]
    max_V = Int64[]
    curr_ind = 1
    while curr_ind <= length(T)
        curr_i = I[curr_ind]
        curr_j = J[curr_ind]
        ijtimes = Int64[]
        # Gather timestamps until next static edge
        while curr_ind <= length(T) && I[curr_ind] == curr_i && J[curr_ind] == curr_j
            push!(ijtimes, T[curr_ind])
            curr_ind += 1
        end
        push!(minmax_I, curr_i)
        push!(minmax_J, curr_j)
        push!(min_V, minimum(ijtimes))
        push!(max_V, maximum(ijtimes))
    end

    n = maximum(simplices)
    min_M = convert(SpIntMat, sparse(minmax_I, minmax_J, min_V, n, n))
    max_M = convert(SpIntMat, sparse(minmax_I, minmax_J, max_V, n, n))
    return (min_M, max_M)
end

"""
interval_overlaps
--------------------

interval_overlaps(dataset::HONData)

Compute the number of active interval overlaps in open triangles.

dataset::HONData
    The dataset.
"""
function interval_overlaps(dataset::HONData)
    min_M, max_M = intervals(dataset)
    A, At, B = basic_matrices(dataset)
    simplex_order = simplex_degree_order(At)
    triangle_order = proj_graph_degree_order(B)

    n = size(B)[2]
    overlaps = zeros(Int64, 4, n)
    Threads.@threads for i = 1:n
        for (j, k) in neighbor_pairs(B, triangle_order, i)
            if B[j, k] > 0 && !triangle_closed(A, At, simplex_order, i, j, k)
                tij1, tij2 = min(min_M[i, j], min_M[j, i]), max(max_M[i, j], max_M[j, i])
                tik1, tik2 = min(min_M[i, k], min_M[k, i]), max(max_M[i, k], max_M[k, i])
                tjk1, tjk2 = min(min_M[j, k], min_M[k, j]), max(max_M[j, k], max_M[k, j])

                # Whether or not each interval overlaps
                ij_ik_ol = tij1 <= tik2 && tik1 <= tij2
                ij_jk_ol = tij1 <= tjk2 && tjk1 <= tij2
                ik_jk_ol = tik1 <= tjk2 && tjk1 <= tik2
                num_overlap = ij_ik_ol + ij_jk_ol + ik_jk_ol

                overlaps[num_overlap + 1, i] += 1
            end
        end
    end

    overlaps = sum(overlaps, dims=2)
    tot = sum(overlaps)
    frac_overlaps = overlaps / tot
    println("dataset & # open triangles & 0 overlaps & 1 overlap & 2 overlaps & 3 overlaps")
    println(@sprintf("%s & %d & %0.3f & %0.3f & %0.3f & %0.3f", dataset.name, tot, frac_overlaps...))
end


text:
using Test
using DiffEqOperators, OrdinaryDiffEq, LinearAlgebra

@testset "KdV equation (Single Solition)" begin
    N = 21
    Δx = 1/(N-1)
    r = 0.5

    # x = 10:Δx:30;
    x = -10:Δx:10;
    # ϕ(x,t) = (r/2)*sech.((sqrt(r)*(x-r*t)/2)-7).^2 # solution of the single forward moving wave
    ϕ(x,t) = (1/2)*sech.((x .- t)/2).^2 # solution of the single forward moving wave

    u0 = ϕ(x,0);
    oriu = zeros(size(x));

    #const du3 = zeros(size(x));
    du = zeros(size(x));
    # const temp = zeros(size(x));

    # A = CenteredDifference(1,2,Δx,length(x),:Dirichlet0,:Dirichlet0);
    A = UpwindDifference{Float64}(1,3,Δx,length(x),-1);
    # C = CenteredDifference(3,2,Δx,length(x),:Dirichlet0,:Dirichlet0);
    #C = UpwindDifference{Float64}(3,3,Δx,length(x),-1);

    function KdV(du, u, p, t)
        # bc = DirichletBC(ϕ(-10-Δx,t),ϕ(10+Δx,t))
        bc = GeneralBC([0,1,-6*ϕ(-10,t),0,-1],[0,1,-6*ϕ(10,t),0,-1],Δx,3)
        #mul!(du3,C,bc*u)
        mul!(du,A,bc*u)
        # @. temp = -0.5*u*du - 0.25*du3
        # copyto!(du,temp)

    end

    single_solition = ODEProblem(KdV, u0, (0.,5.));
    soln = solve(single_solition,Tsit5(),abstol=1e-6,reltol=1e-6);

    #=
    for t in 0:0.5:5
        @show maximum(soln(t)-ϕ(x,t))
    end
    =#
    for t in 0:0.5:5
        @test soln(t) ≈ ϕ(x,t) atol = 0.01;
    end

    # Using Biased Upwinds with 1 offside point
    A2 = UpwindDifference{Float64}(1,3,Δx,length(x),-1,offside=1);
    function KdV(du, u, p, t)
        bc = GeneralBC([0,1,-6*ϕ(-10,t),0,-1],[0,1,-6*ϕ(10,t),0,-1],Δx,3)
        mul!(du,A2,bc*u)
    end
    single_solition = ODEProblem(KdV, u0, (0.,5.));
    soln = solve(single_solition,Tsit5(),abstol=1e-6,reltol=1e-6);
    for t in 0:0.5:5
        @test soln(t) ≈ ϕ(x,t) atol = 0.01;
    end

    A3 = UpwindDifference{Float64}(1,3,Δx*ones(length(x)+1),length(x),-1,offside=1);
    function KdV(du, u, p, t)
        bc = GeneralBC([0,1,-6*ϕ(-10,t),0,-1],[0,1,-6*ϕ(10,t),0,-1],Δx,3)
        mul!(du,A3,bc*u)
    end
    single_solition = ODEProblem(KdV, u0, (0.,5.));
    soln = solve(single_solition,Tsit5(),abstol=1e-6,reltol=1e-6);
    for t in 0:0.5:5
        @test soln(t) ≈ ϕ(x,t) atol = 0.01;
    end

    # Using Biased Upwinds with 2 offside points
    A4 = UpwindDifference{Float64}(1,4,Δx,length(x),-1,offside=2);
    function KdV(du, u, p, t)
        bc = GeneralBC([0,1,-6*ϕ(-10,t),0,-1],[0,1,-6*ϕ(10,t),0,-1],Δx,3)
        mul!(du,A4,bc*u)
    end
    single_solition = ODEProblem(KdV, u0, (0.,5.));
    soln = solve(single_solition,Tsit5(),abstol=1e-6,reltol=1e-6);
    for t in 0:0.5:5
        @test soln(t) ≈ ϕ(x,t) atol = 0.01;
    end

    A5 = UpwindDifference{Float64}(1,4,Δx,length(x),1,offside=2);
    function KdV(du, u, p, t)
        bc = GeneralBC([0,1,-6*ϕ(-10,t),0,-1],[0,1,-6*ϕ(10,t),0,-1],Δx,3)
        mul!(du,-1*A5,bc*u)
    end
    single_solition = ODEProblem(KdV, u0, (0.,5.));
    soln = solve(single_solition,Tsit5(),abstol=1e-6,reltol=1e-6);
    for t in 0:0.5:5
        @test soln(t) ≈ ϕ(x,t) atol = 0.01;
    end
end

# Conduct interesting experiments by referring to
# http://lie.math.brocku.ca/~sanco/solitons/kdv_solitons.php
#  TODO The true solution for this case seems unstable, look for a workaround
#= @testset "KdV equation (Double Solition)" begin
    N = 10
    Δx = 1/(N-1)

    x = -50:Δx:50;
    c1,c2 = 25,16
    # ϕ1(x,t) = 3*p1^2*sech(.5*(p1*(x+2+t))).^2+3*p2^2*sech(.5*(p2*(x+1+t))).^2;

    function ϕ(x,t)
        # t = t-10
        num1 = 2(c1-c2)*(c1*cosh.(√c2*(x .- c2*t)/2).^2 + c2*sinh.(√c1*(x .- c1*t)/2).^2)
        den11 = (√c1-√c2)*cosh.((√c1*(x .- c1*t) + √c2*(x .- c2*t))/2)
        den12 = (√c1+√c2)*cosh.((√c1*(x .- c1*t) - √c2*(x .- c2*t))/2)
        den1 = (den11+den12).^2
        return num1./den1
    end

    # const du3 = zeros(size(x));
    du = zeros(size(x));
    # const temp = zeros(size(x));

    A = UpwindDifference{Float64}(1,3,Δx,length(x),-1);

    # C = UpwindDifference{Float64}(3,1,Δx,length(x),-1);

    function KdV(du, u, p, t)
        bc = GeneralBC([0,1,-6*ϕ(-50,t),0,-1],[0,1,-6*ϕ(50,t),0,-1],Δx,3)
        # mul!(du3,C,bc*u)
        mul!(du,A,bc*u)
        # @. temp = -6*u*du - du3
        # copyto!(du,temp)
    end

    u0 = ϕ(x,0);
    double_solition = ODEProblem(KdV, u0, (0.,5.));
    soln = solve(double_solition,Tsit5());

    # The solution is a forward moving soliton wave with speed = 1
    for t in 0:0.1:9
        @test_skip soln(t) ≈ ϕ(x,t) atol = 0.01
    end
end =#


text:
-- prints -3
x = 2 * (1 + 3) / 4 - 5
print(x)

max_stars_repo_name: JCGrant/JLang

max_issues_repo_name: JCGrant/JLang

max_forks_repo_name: JCGrant/JLang


text:
function draw_raster_n(rhd::RHD2000,han::Gui_Handles,num_chan::Int64)

    num_plot=num_chan

    k_in=num_chan*(han.num16)-num_chan+1

    if (k_in+num_chan-1)>size(rhd.v,2)
        num_chan=(size(rhd.v,2)-k_in+1)
    end

    maxid=find_max_id(rhd,han,k_in,num_chan)
    ctx=Gtk.getgc(han.c)

    xwidth=width(ctx)
    myheight=height(ctx)
    if num_plot<64
        yheight=300
        myoff=myheight-300.0
    else
        yheight=myheight
        myoff=0.0
    end

    ctx_step=yheight/num_chan

    @inbounds for thisid=1:maxid
        count=1
        for i=k_in:(k_in+num_chan-1)
            if han.enabled[i]
                for g=1:rhd.nums[i] #plotting every spike, do we need to do that?
                    if rhd.buf[g,i].id==thisid
                        offset=myoff+ctx_step*(count-1)
                        move_to(ctx,han.draws,offset)
                        line_to(ctx,han.draws,offset+ctx_step-3.0)
                    end
                end
            end
            count+=1
        end

        set_line_width(ctx,0.5);
        @inbounds select_color(ctx,thisid)
        stroke(ctx)
    end

    nothing
end

draw_raster16(rhd::RHD2000,han::Gui_Handles)=draw_raster_n(rhd,han,16)

draw_raster32(rhd::RHD2000,han::Gui_Handles)=draw_raster_n(rhd,han,32)

draw_raster64(rhd::RHD2000,han::Gui_Handles)=draw_raster_n(rhd,han,64)

function update_isi(rhd::RHD2000,han::Gui_Handles,i)

    spike_num=han.sc.spike
    clus=rhd.buf[i,spike_num].id
    mytime=rhd.time[rhd.buf[i,spike_num].inds.start,1]

    my_isi = mytime - han.isi_last_time[clus]

    han.isi_last_time[clus] = mytime

    han.isi[han.isi_ind] = my_isi

    han.isi_clus_ID[han.isi_ind] = clus

    han.isi_ind+=1
    if han.isi_ind>500
        han.isi_ind=1
    end

    han.isi_count+=1
    if han.isi_count>500
        han.isi_count=500
    end

    nothing
end

function draw_isi(rhd::RHD2000,han::Gui_Handles)

    ctx=Gtk.getgc(han.sc.c3)

    myheight=height(ctx)
    mywidth=width(ctx)

    total_clus = max(han.sc.total_clus+1,5)

    for i=1:han.sc.total_clus

        mycount=0
        myviolation=0
        for j=1:han.isi_count
            if han.isi_clus_ID[j] == i+1
                if han.isi[j]/rhd.sr < .0025
                   myviolation += 1
                end
                mycount+=1

                ind=1
                this_isi = han.isi[j]/rhd.sr
                for k=.001:.001:.049
                    if (this_isi > k - .001)&(this_isi < k)
                        han.isi_hist[ind]+=1
                    end
                    ind+=1
                end
            end
        end

        isi_f=round(myviolation/mycount*100,2)
        if isi_f>5.0
            set_source_rgb(ctx,1.0,0.0,0.0)
        else
            set_source_rgb(ctx,1.0,1.0,1.0)
        end
        startx=(i-1)*(han.sc.wave_points)+1
        Cairo.scale(ctx,mywidth/(han.sc.wave_points*total_clus),1)
        move_to(ctx,startx,10)
        show_text(ctx,string(isi_f))

        #Find Scale for histogram
        isi_scale = (maximum(han.isi_hist)>0) ? 50/maximum(han.isi_hist) : 1

        select_color(ctx,i+1)
        move_to(ctx,startx,130)
        line_to(ctx,startx,130-han.isi_hist[1]*isi_scale)
        han.isi_hist[1]=0

        for j=2:length(han.isi_hist)
            if j < han.wave_points
                move_to(ctx,startx+j,130)
                line_to(ctx,startx+j,130-han.isi_hist[j]*isi_scale)
            end
            han.isi_hist[j]=0
        end
        stroke(ctx)

        identity_matrix(ctx)
    end

    nothing
end

#=
Reset Canvas
=#

function clear_c(han::Gui_Handles)

    #paint it black
    ctx = Gtk.getgc(han.c)
    set_source_rgb(ctx,0.0,0.0,0.0)
    paint(ctx)

    if han.c_right_top==1
        prepare_16(ctx,han)
    elseif han.c_right_top==2
        prepare_32(ctx,han)
    elseif han.c_right_top==3
        prepare_64(ctx,han)
    else
    end

    if han.c_right_bottom==1
        prepare_events(ctx,han)
    elseif han.c_right_bottom==2
    elseif han.c_right_bottom==3
    elseif han.c_right_bottom==4
        prepare_scope(ctx,han)
    elseif han.c_right_bottom==5
    elseif han.c_right_bottom==7
        prepare_spectrogram(ctx,han)
    else

    end
    nothing
end

get_bounds(mydim,divs)=linspace(0.0,mydim,divs)

function prepare_raster16(ctx)
end

function prepare_raster32(ctx)
end

function prepare_events(ctx,han)

    myheight=height(ctx)

    for y in collect((myheight-250.0):50.0:(myheight-50.0))
        line(ctx,1.0,500.0,y,y)
    end
    set_source_rgb(ctx,1.0,1.0,1.0)
    set_line_width(ctx,1.0)
    stroke(ctx)

    for i=1:6
	if han.events[i]>-1
            move_to(ctx,10,(myheight-340)+i*50)
	    if han.events[i]<8 #analog
		show_text(ctx,string("A",han.events[i]))
	    else
		show_text(ctx,string("TTL",han.events[i]-8))
	    end
	end
    end

    nothing
end

function prepare_scope(ctx,han)

    myheight=height(ctx)

    dashes = [10.0,10.0,10.0]
    set_dash(ctx, dashes, 0.0)

    for y in collect((myheight-250.0):50.0:(myheight-50.0))
        line(ctx,1.0,500.0,y,y)
    end

    for x in collect(50:50:450)
        line(ctx,x,x,myheight,myheight-300)
    end

    set_source_rgba(ctx,1.0,1.0,1.0,0.3)
    stroke(ctx)

    han.soft.ctx=copy(ctx)

    set_dash(ctx,Float64[])
end

function draw_c3(rhd::RHD2000,han::Gui_Handles)

    ctx=Gtk.getgc(han.sc.c3)
    mywidth=width(ctx)

    spike_num=han.sc.spike
    reads=han.draws

    Cairo.scale(ctx,mywidth/500.0,1.0)

    for i=1:rhd.nums[spike_num]

        @inbounds move_to(ctx,reads,(rhd.buf[i,spike_num].id-1)*10.0+150.0)
        @inbounds line_to(ctx,reads,(rhd.buf[i,spike_num].id-1)*10.0+160.0)
        set_line_width(ctx,0.5);
	    @inbounds select_color(ctx,rhd.buf[i,spike_num].id)
	    stroke(ctx)
    end

    identity_matrix(ctx)

    if reads==1
        prepare_c3(han.sc)
        draw_templates_c3(han.sc)
        draw_isi(rhd,han)
    end

    nothing
end

function prepare_c3(sc::SpikeSorting.Single_Channel)

    ctx=Gtk.getgc(sc.c3)

    myheight=height(ctx)
    mywidth=width(ctx)

    total_clus = max(sc.total_clus+1,5)

    line(ctx,0,mywidth,130,130)
    set_source_rgb(ctx,1.0,1.0,1.0)
    stroke(ctx)

    Cairo.scale(ctx,mywidth/(sc.wave_points*total_clus),1)

    for i=1:total_clus

        line(ctx,i*sc.wave_points,i*sc.wave_points,0,130)
        set_source_rgb(ctx,1.0,1.0,1.0)
        stroke(ctx)

    end

    identity_matrix(ctx)

    if sc.buf.selected_clus>0
        (x1_f,x2_f,y1_f,y2_f)=SpikeSorting.get_template_dims(sc,sc.buf.selected_clus)
        draw_box(x1_f,y1_f,x2_f,y2_f,(1.0,0.0,1.0),1.0,ctx)
    end

    nothing
end