Merge pull request #35 from timcop/tests

Tests and bug fixes - Nice work Tim!

.gitignore

@@ -4,3 +4,7 @@
 *.DS_Store
 docs/build/
 Manifest.toml
+lcov.info
+.vscode/settings.json
+.vscode/*
+lcov.info

.vscode/settings.json

@@ -0,0 +1,2 @@
+{
+}

Project.toml

@@ -4,18 +4,21 @@ authors = ["Ashton Bradley <ashton.bradley@gmail.com>"]
 version = "0.2.7"
 
 [deps]
+Coverage = "a2441757-f6aa-5fb2-8edb-039e3f45d037"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 SimpleWeightedGraphs = "47aef6b3-ad0c-573a-a1e2-d07658019622"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24"
+UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [compat]
 FFTW = "1, 1.4"
@@ -29,3 +32,10 @@ SimpleWeightedGraphs = "1, 1.1"
 SpecialFunctions = "1"
 ToeplitzMatrices = "0.7"
 julia = "1.3"
+
+[extras]
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Test", "SafeTestsets"]

coverage.jl

@@ -0,0 +1,6 @@
+using Coverage, Pkg 
+# process '*.cov' files
+Pkg.test(; coverage=true)
+coverage = process_folder() # defaults to src/; alternatively, supply the folder name as argument
+LCOV.writefile("lcov.info", coverage)
+clean_folder("src")

src/creation.jl

@@ -236,53 +236,53 @@ function gpecore_exact(K,L=2,N=100,R = K)
     #Κ = 1
     #R = 2 # Stretching coordinate.  R ~ kappa seems to be a good ballpark
 
-# Convert Chebyshev grid to [0,L] from [1 -1];
-z,Dz = getChebDMatrix(N)
-blank,D2z = getChebD2Matrix(N)
-z = vec(reverse(z,dims=2))
-z = (z)*L./2
-Dz = -2*Dz./L
-D2z = 4*D2z./L.^2
-
-# y is the physical coordinate, range [0 infinity] (i.e y = r)
-y =  @. R*(1+z)/(1-z)
-
-#initial guess based on ansatz
-ψ = @. y/hypot(1/Λ,y)
-ψ[1] = 0
-ψ[end] = 1
-ψ0 = vec(ψ)
-
-Q = vec(z.^2 .- 2*z .+ 1)
-Qmat = repeat(Q,1,N)
-Zmat = repeat(2*(z .-1),1,N)
-Ymat = repeat(1 ./y,1,N);
-
-# Second Derivative
-residuals = -0.5*( (Q.*(D2z*ψ) + 2*(z .-1).*(Dz*ψ) ).*Q/(4*R^2)
-        + (Q/(2*R) ./y).*(Dz*ψ) )+ 0.5*K^2 *ψ ./y.^2 + ψ.^3 .- ψ
-residuals[1] = 0; residuals[end] =0
-
-while sum(abs.(residuals).^2) > 1e-12
+    # Convert Chebyshev grid to [0,L] from [1 -1];
+    z,Dz = getChebDMatrix(N)
+    blank,D2z = getChebD2Matrix(N)
+    z = vec(reverse(z,dims=2))
+    z = (z)*L./2
+    Dz = -2*Dz./L
+    D2z = 4*D2z./L.^2
+
+    # y is the physical coordinate, range [0 infinity] (i.e y = r)
+    y =  @. R*(1+z)/(1-z)
+
+    #initial guess based on ansatz
+    ψ = @. y/hypot(1/Λ,y)
+    ψ[1] = 0
+    ψ[end] = 1
+    ψ0 = vec(ψ)
+
+    Q = vec(z.^2 .- 2*z .+ 1)
+    Qmat = repeat(Q,1,N)
+    Zmat = repeat(2*(z .-1),1,N)
+    Ymat = repeat(1 ./y,1,N);
+
     # Second Derivative
     residuals = -0.5*( (Q.*(D2z*ψ) + 2*(z .-1).*(Dz*ψ) ).*Q/(4*R^2)
             + (Q/(2*R) ./y).*(Dz*ψ) )+ 0.5*K^2 *ψ ./y.^2 + ψ.^3 .- ψ
     residuals[1] = 0; residuals[end] =0
 
-    Jacobi = -0.5*( (Qmat.*(D2z) + Zmat.*(Dz) ).*(Qmat/(4*R^2))
-    + (Qmat/(2*R).*Ymat).*(Dz) ) + diagm(0 => 0.5*K^2 ./y.^2)+ diagm(0 => 3*ψ.^2 .- 1)
+    while sum(abs.(residuals).^2) > 1e-12
+        # Second Derivative
+        residuals = -0.5*( (Q.*(D2z*ψ) + 2*(z .-1).*(Dz*ψ) ).*Q/(4*R^2)
+                + (Q/(2*R) ./y).*(Dz*ψ) )+ 0.5*K^2 *ψ ./y.^2 + ψ.^3 .- ψ
+        residuals[1] = 0; residuals[end] =0
 
+        Jacobi = -0.5*( (Qmat.*(D2z) + Zmat.*(Dz) ).*(Qmat/(4*R^2))
+        + (Qmat/(2*R).*Ymat).*(Dz) ) + diagm(0 => 0.5*K^2 ./y.^2)+ diagm(0 => 3*ψ.^2 .- 1)
 
-    Jacobi[1,:] = [1 zeros(1,N-1)]
-    Jacobi[N,:] = [zeros(1,N-1) 1]
 
+        Jacobi[1,:] = [1 zeros(1,N-1)]
+        Jacobi[N,:] = [zeros(1,N-1) 1]
 
-    Δ = vec(-4/7*(Jacobi\residuals))
 
-    ψ = ψ + Δ
-    ψ[1] = 0
-    ψ[end] =1
-end
+        Δ = vec(-4/7*(Jacobi\residuals))
+
+        ψ = ψ + Δ
+        ψ[1] = 0
+        ψ[end] =1
+    end
     res = norm(residuals)^2
     return y,ψ,res
 end
@@ -324,7 +324,7 @@ function chebdif(N, M)
 
    n1,n2 = (floor(N/2)|> Int,ceil(N/2)|> Int)  # Indices used for flipping trick.
 
-   k = collect(0:N-1)' # Compute theta vector.
+   k = collect(0.0:N-1)' # Compute theta vector.
    th = k*π/(N-1)
 
    x = sin.(π*collect(N-1:-2:1-N)'./(2*(N-1))) # Compute Chebyshev points.
@@ -335,8 +335,10 @@ function chebdif(N, M)
    DX[diagind(DX)] .= 1 # Put 1's on the main diagonal of DX.
 
    C = Matrix(Toeplitz((-1).^k',(-1).^k'))  # C is the matrix with
-   C[1,:] = C[1,:]*2; C[N,:] = C[N,:]*2  # entries c(k)/c(j)
-   C[:,1] = C[:,1]/2; C[:,N] = C[:,N]/2
+   C[1,:] = C[1,:]*2; 
+   C[N,:] = C[N,:]*2  # entries c(k)/c(j)
+   C[:,1] = C[:,1]/2; 
+   C[:,N] = C[:,N]/2;
 
    Z = 1 ./DX  # Z contains entries 1/(x(k)-x(j))
    Z[diagind(Z)] .= 0  # with zeros on the diagonal.
@@ -352,7 +354,7 @@ function chebdif(N, M)
    end
 
    return x, DM
-   end
+end
 
 function getChebDMatrix(n)
    z, M = chebdif(n, 1)

src/detection.jl

@@ -64,7 +64,7 @@ function findvortices_grid(psi::Sphere;shift=true)
     w2 = sum(windvals[end,:])
     (sign(w2) != 0) && (rawvort = [rawvort; pi 0.0 w2])
 
-    vort = sortslices(vort,dims=1)
+    vort = sortslices(rawvort,dims=1)
     return PointVortex(vort)
 end
 

src/get_dipoles.jl

@@ -1,50 +1,50 @@
 function get_dipoles(xi,yi,n,nplus,Lx,Ly)
 
-if isempty(xi) || isempty(yi)
-  @info "Input position vector is empty"
-  xd = []
-  yd = []
-  is_dipole = []
-  num_found_in_dipoles = 0
-  return xd,yd,is_dipole,num_found_in_dipoles
-end
+  if isempty(xi) || isempty(yi)
+    @info "Input position vector is empty"
+    xd = []
+    yd = []
+    is_dipole = []
+    num_found_in_dipoles = 0
+    return xd,yd,is_dipole,num_found_in_dipoles
+  end
 
-κi = [ones(nplus); -ones(n-nplus)]
+  κi = [ones(nplus); -ones(n-nplus)]
 
-# x and y distances
-xij,yij = distances(xi,yi)
+  # x and y distances
+  xij,yij = distances(xi,yi)
 
-# incorporate periodicity
-periodic_distances!(xij,yij,Lx,Ly)
+  # incorporate periodicity
+  periodic_distances!(xij,yij,Lx,Ly)
 
-# inter-vortex distance matrix, set self distance to much larger than box size
-rij = hypot.(xij,yij)
-rij += diagm(0=>1e6*hypot(Lx,Ly)*ones(n))
+  # inter-vortex distance matrix, set self distance to much larger than box size
+  rij = hypot.(xij,yij)
+  rij += diagm(0=>1e6*hypot(Lx,Ly)*ones(n))
 
-# get index of minimum separations
-_, nn_index = findmin(rij,dims=2)
+  # get index of minimum separations
+  _, nn_index = findmin(rij,dims=2)
 
-nn_mat = ones(n)*collect(1:n)' .|> Int
-nn_vec = nn_mat[nn_index][:,1]
+  nn_mat = ones(n)*collect(1:n)' .|> Int
+  nn_vec = nn_mat[nn_index][:,1]
 
-mutuals = (collect(1:n) .== nn_vec[nn_vec])
+  mutuals = (collect(1:n) .== nn_vec[nn_vec])
 
-# find dipoles
-is_dipole = zeros(n)
-for i in eachindex(mutuals)
-is_dipole[i] = (mutuals[i] == 1 && κi[i]*κi[nn_vec[i]] == -1.0)
-end
+  # find dipoles
+  is_dipole = zeros(n)
+  for i in eachindex(mutuals)
+    is_dipole[i] = (mutuals[i] == 1 && κi[i]*κi[nn_vec[i]] == -1.0)
+  end
 
-dipoles_index = findall(is_dipole .> 0.0)
-dipoles_index2 = nn_vec[dipoles_index]
-num_found_in_dipoles = length(dipoles_index)
+  dipoles_index = findall(is_dipole .> 0.0)
+  dipoles_index2 = nn_vec[dipoles_index]
+  num_found_in_dipoles = length(dipoles_index)
 
-xd = xi[dipoles_index]
-yd = yi[dipoles_index]
-xd2 = xi[dipoles_index2]
-yd2 = yi[dipoles_index2]
-xd = [xd[1:Int(end/2)]; xd2[1:Int(end/2)]]
-yd = [yd[1:Int(end/2)]; yd2[1:Int(end/2)]]
+  xd = xi[dipoles_index]
+  yd = yi[dipoles_index]
+  xd2 = xi[dipoles_index2]
+  yd2 = yi[dipoles_index2]
+  xd = [xd[1:Int(end/2)]; xd2[1:Int(end/2)]]
+  yd = [yd[1:Int(end/2)]; yd2[1:Int(end/2)]]
 
-return xd,yd,is_dipole,num_found_in_dipoles
+  return xd,yd,is_dipole,num_found_in_dipoles
 end

src/utils.jl

@@ -21,6 +21,26 @@ function found_near(n)
     return near
 end
 
+function found_near(n, nvorts)
+    near = true
+    for j in 1:n
+        psi,vort = rand_vortexfield(nvorts)
+        vortfound = findvortices(psi)
+        vfdata = vortex_array(vortfound)
+        vdata = vortex_array(vort)
+        if length(vfdata[:, 1]) != nvorts
+            return false
+        end
+        dx = Δ(psi.x)
+        sort!(vdata, dims=1)
+        sort!(vfdata, dims=1)
+        for i in 1:nvorts
+            near *= isapprox(vdata[i, :], vfdata[i, :], rtol = dx/4)
+        end
+    end
+    return near
+end
+
 """
     vortices = remove_vortices_edge(vort::Array{PointVortex,1},x,y,edge=1)
 

test/creations_test.jl

@@ -0,0 +1,12 @@
+using VortexDistributions
+
+## Basic unit testing for coverage 
+@test scalar_ansatz(1) == scalar_ansatz(-1)
+f = Ansatz()
+x = LinRange(0, 10, 100)
+y = f.(x)
+@test length(y) == length(x)
+@test f.(3,4) == f(5)
+
+y, ψ, res = VortexDistributions.gpecore_exact(1,2,100)
+@test length(y) == 100

test/keep_vortices_test.jl

@@ -0,0 +1,23 @@
+using VortexDistributions
+
+# simple test field
+Nx = 400; Ny = 400
+Lx = 200; Ly = Lx
+x = LinRange(-Lx/2,Lx/2, Nx+1)[1:end-1]; y = LinRange(-Ly/2,Ly/2, Ny+1)[1:end-1];
+psi0 = one.(x*y') |> complex
+
+psi = Torus(copy(psi0), x, y);
+
+
+xp = 0.0
+yp = 0.0
+vp = PointVortex(xp,yp,1)
+sp = ScalarVortex(vp)
+
+vortex!(psi, sp)
+
+vfound = findvortices(psi)
+
+vortm = keep_vortices(vfound)
+
+@test length(vfound) == length(vortm)

test/multivort_test.jl

@@ -0,0 +1,6 @@
+using VortexDistributions, Test
+
+
+# @testset "Multivortex creation and detection" begin
+@test found_near(30)
+# end

test/multivort_test2.jl

@@ -0,0 +1,6 @@
+using VortexDistributions, Test
+
+
+# @testset "Multivortex creation and detection" begin
+@test found_near(30, 10)
+# end