From 6d8cca376333d874ebae1619973fd8c274c61699 Mon Sep 17 00:00:00 2001 From: bblankrot <34473471+bblankrot@users.noreply.github.com> Date: Wed, 28 Mar 2018 16:10:14 +0200 Subject: [PATCH 1/2] minimumNP --- docs/src/minimalNP.md | 2 +- src/minimum_N_P.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/minimalNP.md b/docs/src/minimalNP.md index 22414ed..ad15ad0 100644 --- a/docs/src/minimalNP.md +++ b/docs/src/minimalNP.md @@ -25,7 +25,7 @@ Specifically, this error is measured on a circle of radius R_multipole*maximum(hypot.(s.ft[:,1],s.ft[:,2])) ``` -as this is the scattering disk on which the translation to cylindrical harmonics +as this is the scattering disc on which the translation to cylindrical harmonics (Bessel & Hankel functions) are performed (and beyond which any gain or loss of accuracy due to `N` is mostly irrelevant). diff --git a/src/minimum_N_P.jl b/src/minimum_N_P.jl index bfd7fdb..8308925 100644 --- a/src/minimum_N_P.jl +++ b/src/minimum_N_P.jl @@ -10,7 +10,7 @@ Error is calculated on `N_points` points on the scattering disk (`s.R`), by assuming a fictitious line source and comparing its field to that produced by the resulting potential densities. -Since for moderate wavelengths and errors, \$\varepsilon \\propto N^3\$, we +Since the error scales with \$N^{-3}\$ for moderate wavelengths and errors, we estimate `N` using the error of `N_start`, then binary search based on that guess between `N_min` and `N_max`. """ From 9ccc7bb271040e3a22c65261e8163743cde060f1 Mon Sep 17 00:00:00 2001 From: bblankrot <34473471+bblankrot@users.noreply.github.com> Date: Thu, 5 Apr 2018 13:14:07 +0200 Subject: [PATCH 2/2] Code cleanup + more tests --- src/ParticleScattering.jl | 8 +- src/multipole.jl | 3 +- src/optimize_phis.jl | 4 +- src/optimize_rs_old.jl | 157 -------------------------------------- src/shapes.jl | 2 +- src/visualization.jl | 2 +- test/fmm_test.jl | 4 +- test/multipole_test.jl | 56 ++++++++++++++ test/runtests.jl | 4 +- test/utilities_test.jl | 13 ++++ 10 files changed, 81 insertions(+), 172 deletions(-) delete mode 100644 src/optimize_rs_old.jl create mode 100644 test/multipole_test.jl create mode 100644 test/utilities_test.jl diff --git a/src/ParticleScattering.jl b/src/ParticleScattering.jl index 2e9e86e..7dfa316 100644 --- a/src/ParticleScattering.jl +++ b/src/ParticleScattering.jl @@ -11,7 +11,7 @@ using IterativeSolvers, LinearMaps, Optim import LineSearches #For plotting with PyPlot using PyPlot, PyCall -@pyimport matplotlib.patches as patch #circles,polygons +@pyimport matplotlib.patches as patch #circles, polygons include("PS_types.jl") include("shapes.jl") @@ -44,9 +44,6 @@ export ScatteringProblem, OptimBuffer, FMMoptions, R_multipole, ShapeParams, CircleParams, AbstractShapeParams #methods, utilities.jl export find_border, uniqueind - -### documented till here#### - #methods, optimize_phis.jl export optimize_φ #methods, optimize_rs.jl @@ -60,10 +57,7 @@ include("visualization_pgf.jl") export plotNearField_pgf, drawShapes_pgf #temp -include("optimize_rs_old.jl") - export divideSpace, FMMtruncation, particleExpansion, FMMbuildMatrices - export FMMbuffer, FMMmatrices end diff --git a/src/multipole.jl b/src/multipole.jl index 14ca5f8..d1bcdec 100644 --- a/src/multipole.jl +++ b/src/multipole.jl @@ -185,7 +185,8 @@ function innerFieldCircle(kin, gamma, center::Array{Float64,1}, points::Array{Fl bess = [besselj(p,kin*rs_moved[ii]) for ii=1:len, p=0:P] E = gamma[P + 1]*bess[:,1] for p = 1:P - E += bess[:,p+1].*(gamma[p + P + 1]*exp.(1.0im*p*ts_moved) + (-1)^p*gamma[-p + P + 1]*exp.(-1.0im*p*ts_moved)) + E += bess[:,p+1].*(gamma[p + P + 1]*exp.(1.0im*p*ts_moved) + + (-1)^p*gamma[-p + P + 1]*exp.(-1.0im*p*ts_moved)) end return E end diff --git a/src/optimize_phis.jl b/src/optimize_phis.jl index b3ab802..7a1ff30 100644 --- a/src/optimize_phis.jl +++ b/src/optimize_phis.jl @@ -3,8 +3,8 @@ optimopts::Optim.Options, minimize = true) Optimize the rotation angles of a particle collection for minimization or -maximization of the field intensity at `points`, depending on `minimize`. -`optimopts` and `method` define the optimization emthod, convergence criteria, +maximization (depending on `minimize`) of the field intensity at `points`. +`optimopts` and `method` define the optimization method, convergence criteria, and other optimization parameters. Returns an object of type `Optim.MultivariateOptimizationResults`. """ diff --git a/src/optimize_rs_old.jl b/src/optimize_rs_old.jl deleted file mode 100644 index 1c5056a..0000000 --- a/src/optimize_rs_old.jl +++ /dev/null @@ -1,157 +0,0 @@ -function optimize_radius_dep(rs0, r_min, r_max, points, P, θ_i, k0, kin, centers, fmmopts, optimopts, minimize, method, linesearch) - #assumes all are circles with their own IDs. - #TODO: separate ids from rs to allow symmetry. - #setup FMM reusal - Ns = size(centers,1) - (groups, boxSize) = divideSpace(centers, fmmopts) - - (P2, Q) = FMMtruncation(fmmopts.acc, boxSize, k0) - mFMM = FMMbuildMatrices(k0, P, P2, Q, groups, centers, boxSize, tri=true) - - ids = collect(1:Ns) #all different sizes - scatteringMatrices = Array{SparseMatrixCSC{Complex{Float64},Int}}(0) - #allocate derivative - dS_S = Array{SparseMatrixCSC{Complex{Float64},Int}}(0) - for ic = 1:Ns - push!(scatteringMatrices,speye(2*P+1)) - push!(dS_S,speye(2*P+1)) - end - - buf = FMMbuffer(Ns,P,Q,length(groups)) - - #stuff that is done once - points_moved = Array{Float64}(2) - H = Array{Complex{Float64}}(Ns*(2*P+1),size(points,1)) - for i = 1:size(points,1), ic = 1:Ns - points_moved[1] = points[i,1] - centers[ic,1] - points_moved[2] = points[i,2] - centers[ic,2] - r_angle = atan2(points_moved[2],points_moved[1]) - kr = k0*sqrt(points_moved[1]^2 + points_moved[2]^2) - for l = -P:P - H[(ic-1)*(2*P+1) + P + 1 + l,i] = besselh(l,kr)*exp(1.0im*l*r_angle) - end - end - α_inc = Complex{Float64}[exp(1.0im*p*(pi/2-θ_i)) for p=-P:P] - φs = zeros(Float64,Ns) - - # Allocate buffers - #TODO: accept vecotr and scalar r_min,r_max - shared_var = OptimBuffer(Ns,P,size(points,1)) - if length(rs0) == 1 - initial_rs = rs0*ones(Float64,Ns) - else - initial_rs = copy(rs0) - end - last_rs = similar(initial_rs) - - if minimize - df = OnceDifferentiable(rs -> optimize_radius_f_dep(rs, last_rs, shared_var, φs, α_inc, H, points, P, θ_i, Ns, k0, kin, centers,scatteringMatrices, dS_S, ids, mFMM, fmmopts, buf), - (grad_stor, rs) -> optimize_radius_g!_dep(grad_stor, rs, last_rs, shared_var, φs, α_inc, H, points, P, θ_i, Ns, k0, kin, centers, scatteringMatrices, dS_S, ids, mFMM, fmmopts, buf, minimize)) - else - df = OnceDifferentiable(rs -> -optimize_radius_f_dep(rs, last_rs, shared_var, φs, α_inc, H, points, P, θ_i, Ns, k0, kin, centers,scatteringMatrices, dS_S, ids, mFMM, fmmopts, buf), - (grad_stor, rs) -> optimize_radius_g!_dep(grad_stor, rs, last_rs, shared_var, φs, α_inc, H, points, P, θ_i, Ns, k0, kin, centers, scatteringMatrices, dS_S, ids, mFMM, fmmopts, buf, minimize)) - end - - outer_iterations = optimopts.iterations #? - inner_iterations = optimopts.iterations #? - - if method == "LBFGS" - optimize(df, initial_rs, r_min, r_max, Fminbox{LBFGS}(); - optimizer_o = optimopts, iterations = outer_iterations, - linesearch = linesearch, x_tol = optimopts.x_tol, - f_tol = optimopts.f_tol, g_tol = optimopts.g_tol) - elseif method == "BFGS" - optimize(df, initial_rs, r_min, r_max, Fminbox{BFGS}(); - optimizer_o = optimopts, iterations = outer_iterations, - linesearch = linesearch, x_tol = optimopts.x_tol, - f_tol = optimopts.f_tol, g_tol = optimopts.g_tol) - end -end - -function optimize_radius_common!_dep(rs, last_rs, shared_var, φs, α_inc, H, points, P, θ_i, Ns, k0, kin, centers, scatteringMatrices, dS_S, ids, mFMM, opt, buf) - if (rs != last_rs) - copy!(last_rs, rs) - #do whatever common calculations and save to shared_var - #construct rhs - for ic = 1:Ns - rng = (ic-1)*(2*P+1) + (1:2*P+1) - try - updateCircleScatteringDerivative!(scatteringMatrices[ids[ic]], dS_S[ids[ic]], k0, kin, rs[ic], P) - catch - warn("Could not calculate derivatives for ic=$ic,k0=$k0,kin=$kin,R=$(rs[ic])") - error() - end - buf.rhs[rng] = scatteringMatrices[ids[ic]]*α_inc - #phase shift added to move cylinder coords - phase = exp(1.0im*k0*(cos(θ_i)*centers[ic,1] + sin(θ_i)*centers[ic,2])) - buf.rhs[rng] *= phase - end - - if opt.method == "pre" - MVP = LinearMap{eltype(buf.rhs)}((output_, x_) -> FMM_mainMVP_pre!(output_, - x_, scatteringMatrices, φs, ids, P, mFMM, buf.pre_agg, buf.trans), - Ns*(2*P+1), Ns*(2*P+1), ismutating = true) - elseif opt.method == "pre2" - MVP = LinearMap{eltype(buf.rhs)}((output_, x_) -> FMM_mainMVP_pre2!(output_, - x_, scatteringMatrices, φs, ids, P, mFMM, buf.pre_agg, buf.trans), - Ns*(2*P+1), Ns*(2*P+1), ismutating = true) - end - - fill!(shared_var.β,0.0) - shared_var.β,ch = gmres!(shared_var.β, MVP, buf.rhs, restart = Ns*(2*P+1) + 1, maxiter = Ns*(2*P+1), tol = opt.tol, log = true) #no restart, preconditioning - if !ch.isconverged - error("FMM process did not converge") - end - shared_var.f[:] = H.'*shared_var.β - shared_var.f[:] += exp.(1.0im*k0*(cos(θ_i)*points[:,1] + sin(θ_i)*points[:,2])) #incident - end -end - -function optimize_radius_f_dep(rs, last_rs, shared_var, φs, α_inc, H, points, P, θ_i, Ns, k0, kin, centers, scatteringMatrices, dS_S, ids, mFMM, opt, buf) - optimize_radius_common!_dep(rs, last_rs, shared_var, φs, α_inc, H, points, P, θ_i, Ns, k0, kin, centers, scatteringMatrices, dS_S, ids, mFMM, opt, buf) - - func = sum(abs2,shared_var.f) -end - -function optimize_radius_g!_dep(grad_stor, rs, last_rs, shared_var, φs, α_inc, H, points, P, θ_i, Ns, k0, kin, centers, scatteringMatrices, dS_S, ids, mFMM, opt, buf, minimize) - optimize_radius_common!_dep(rs, last_rs, shared_var, φs, α_inc, H, points, P, θ_i, Ns, k0, kin, centers, scatteringMatrices, dS_S, ids, mFMM, opt, buf) - - if opt.method == "pre" - MVP = LinearMap{eltype(buf.rhs)}((output_, x_) -> FMM_mainMVP_pre!(output_, - x_, scatteringMatrices, φs, ids, P, mFMM, buf.pre_agg, buf.trans), - Ns*(2*P+1), Ns*(2*P+1), ismutating = true) - elseif opt.method == "pre2" - MVP = LinearMap{eltype(buf.rhs)}((output_, x_) -> FMM_mainMVP_pre2!(output_, - x_, scatteringMatrices, φs, ids, P, mFMM, buf.pre_agg, buf.trans), - Ns*(2*P+1), Ns*(2*P+1), ismutating = true) - end - - #time for gradient - shared_var.rhs_grad[:] = 0.0 - shared_var.∂β[:] = 0.0 - # tempn = Array{Complex{Float64}}(2*P+1) - for n = 1:Ns - #compute n-th gradient - rng = (n-1)*(2*P+1) + (1:2*P+1) - shared_var.rhs_grad[rng] = dS_S[ids[n]]*shared_var.β[rng] - shared_var.∂β[:,n], ch = gmres!(shared_var.∂β[:,n], MVP, shared_var.rhs_grad, restart = Ns*(2*P+1) + 1, maxiter = Ns*(2*P+1), tol = opt.tol, log = true) - #warn("using dbdn_tol = 10*opt.tol = $(10*opt.tol)") - - if ch.isconverged == false - display(ch) - display("rs:") - display(rs) - display("β:") - display(shared_var.β) - display("rhs_grad:") - display(shared_var.rhs_grad) - display("∂β:") - display(shared_var.∂β[:,n]) - error("FMM process did not converge for partial derivative $n/$Ns. ") - end - #prepare for next one - shared_var.rhs_grad[rng] = 0.0 - end - - grad_stor[:] = ifelse(minimize,2,-2)*real(shared_var.∂β.'*(H*conj(shared_var.f))) -end diff --git a/src/shapes.jl b/src/shapes.jl index 3ca8adc..a2ebd8b 100644 --- a/src/shapes.jl +++ b/src/shapes.jl @@ -126,7 +126,7 @@ function randpoints(M, dmin, width, height, points; failures = 100) dmin2 = dmin^2 for i = 2:size(points,1), j = 1:i-1 dist2 = sum(x -> x^2, points[i,:] - points[j,:]) - dist2 <= dmin2 && error("randpoints: given points have distance > dmin") + dist2 <= dmin2 && error("randpoints: given points have distance <= dmin") end x_res = [rand(Float64)*width] diff --git a/src/visualization.jl b/src/visualization.jl index d026940..60dceb6 100644 --- a/src/visualization.jl +++ b/src/visualization.jl @@ -9,7 +9,7 @@ angle `θ_i` scattering from the ScatteringProblem `sp`, using matplotlib's `pcolormesh`. Can accept number of sampling points in each direction plus bounding box or calculate automatically. -Uses the FMM options given by `opt` (default behavious is disabled FMM); +Uses the FMM options given by `opt` (FMM is disabled by default); `use_multipole` dictates whether electric field is calculated using the multipole/cylindrical harmonics (true) or falls back on potential densities (false). Either way, the multiple-scattering system is solved in the cylindrical diff --git a/test/fmm_test.jl b/test/fmm_test.jl index 8858ee6..b200b8c 100644 --- a/test/fmm_test.jl +++ b/test/fmm_test.jl @@ -14,10 +14,10 @@ φs = 2π*rand(M^2) #random rotation angles dist = 2*maximum(shapes[i].R for i=1:2) try - centers = randpoints(M^2, dist, 5λ0, 5λ0, [0.0 0.0], + centers = randpoints(M^2, dist, 5λ0, 5λ0, [0.0 0.0; 0.01+dist 0.01], failures = 1_000) sp = ScatteringProblem(shapes, ids, centers, φs) - @test verify_min_distance(sp, [0.0 0.0]) + @test verify_min_distance(sp, [0.0 0.0; 0.01+dist 0.01]) catch warn("Could not find random points (1)") end diff --git a/test/multipole_test.jl b/test/multipole_test.jl new file mode 100644 index 0000000..d087c0a --- /dev/null +++ b/test/multipole_test.jl @@ -0,0 +1,56 @@ +@testset "multipole" begin + #let's check that rotation works properly + k0 = 0.1 + kin = 100k0 + λ0 = 2π/k0 + θ_i = 0.0 + + N = 1000 + P = 5 + shapes = [squircle(2π/λ0, N)] + centers = [-1.5λ0 0; 1.5λ0 0] + Ns = size(centers,1) + φs = [0.0; 0.23] + ids = ones(Int, Ns) + sp = ScatteringProblem(shapes, ids, centers, φs) + points = λ0*[5 0; -2 3] + + verify_min_distance(shapes, centers, ids, points) + β1, σ1 = solve_particle_scattering(k0, kin, P, sp, θ_i; + get_inner = true, verbose = false) + Ez1 = scattered_field_multipole(k0, β1, centers, points) + + #now rotate everything by some angle and compare + θ_r = 1.23#2π*rand() + #notation is transposed due to structure of centers + shapes2 = [squircle(2π/λ0, N); CircleParams(15)] #for extra code coverage + centers2 = centers*([cos(θ_r) -sin(θ_r); sin(θ_r) cos(θ_r)].') + φs2 = θ_r + φs + sp2 = ScatteringProblem(shapes2, ids, centers2, φs2) + points2 = points*([cos(θ_r) -sin(θ_r); sin(θ_r) cos(θ_r)].') + + β2, σ2 = solve_particle_scattering(k0, kin, P, sp2, θ_i + θ_r; + get_inner = true, verbose = true) + Ez2 = scattered_field_multipole(k0, β2, centers2, points2) + + @test Ez1 ≈ Ez2 + β2_r = copy(β2) + for ic = 1:Ns, p = -P:P + β2_r[(ic-1)*(2P+1) + p + P + 1] *= exp(1.0im*θ_r*p) + end + @test β1 ≈ β2_r + @test σ1 ≈ σ2 +end + +@testset "circle" begin + k0 = 1 + kin = 0.1 + R = 0.2 + P = 10 + S, gamma = ParticleScattering.circleScatteringMatrix(k0, kin, R, P; gamma = true) + center = [0 2R] + points = center .+ [0 0.5R; 0.1R -0.2R] + E1 = ParticleScattering.innerFieldCircle(kin, gamma, center[1,:], points) + E2 = ParticleScattering.innerFieldCircle(kin, gamma, center[1,:], points[1,:]) + @test E1[1] == E2 +end diff --git a/test/runtests.jl b/test/runtests.jl index 9fa94d5..70803cd 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,7 +2,9 @@ using ParticleScattering using Base.Test include("scatteredfield_test.jl") +include("multipole_test.jl") +include("fmm_test.jl") include("optimize_radius_test.jl") include("optimize_angle_test.jl") include("minimum_N_P_test.jl") -include("fmm_test.jl") +include("utilities_test.jl") diff --git a/test/utilities_test.jl b/test/utilities_test.jl new file mode 100644 index 0000000..4794ca5 --- /dev/null +++ b/test/utilities_test.jl @@ -0,0 +1,13 @@ +@testset "utilities" begin + s = rounded_star(1, 0.1, 7, 100) + @test ParticleScattering.pInPolygon([0.0 0.0], s.ft) + @test !ParticleScattering.pInPolygon([0.84 0.35], s.ft) + @test ParticleScattering.pInPolygon([0.84 0.35], s.ft+0.1) + + sp = ScatteringProblem([s],[1],[0.1 0.1], [15π]) + border = find_border(sp) + border2 = find_border(sp,[0.1 0.1]) + @test border == border2 + border3 = find_border(sp,[0.1 0.9999]) + @test border != border3 +end