Merge pull request #13 from bblankrot/dev

Improved coverage and documentation

.travis.yml

@@ -5,6 +5,9 @@ os:
   - osx
 julia:
   - 0.6
+branches:
+  except:
+    - dev
 notifications:
   email: false
 env:

REQUIRE

@@ -1,8 +1,8 @@
 julia 0.6
 
-IterativeSolvers
+IterativeSolvers 0.5.0
 LinearMaps
-Optim
+Optim 0.14.0
 PyPlot
 DataFrames
 CSV

docs/make.jl

@@ -9,6 +9,7 @@ makedocs(format = :html,
         "Tutorials" => Any[
             "tutorial1.md",
             "tutorial2.md",
+            "tutorial_optim_angle.md",
             "tutorial_optim_radius.md"
         ],
         "Choosing Minimal N and P" => "minimalNP.md",

docs/src/index.md

@@ -27,6 +27,7 @@ PGFPlotsX
 ## Getting Started
 
 Users are encouraged to follow the tutorials, as they provide a gradual
-introduction to this package. Complex and complete examples involving
-optimization are available in the [examples](https://github.com/bblankrot/ParticleScattering.jl/tree/master/examples)
+introduction to this package yet cover most of its functionality.
+Complex and complete examples involving optimization are available in the
+ [examples](https://github.com/bblankrot/ParticleScattering.jl/tree/master/examples)
 folder.

docs/src/tutorial2.md

@@ -1,4 +1,4 @@
-# Tutorial 2: Accelerating Solutions with FMM
+# [Tutorial 2: Accelerating Solutions with FMM](@id tutorial2)
 
 In this tutorial, we examine scattering from several hundreds of particles, and
 use the built-in Fast Multipole Method (FMM) implementation to provide faster
@@ -58,7 +58,7 @@ divideSpace(centers, fmm_options; drawGroups = true)
 
 In this plot, the red markers denote the group centers while stars denote
 particle centers (the particles can be drawn on top of this plot with
-`draw_shapes`). At first, it might look strange that most each particle lies
+`draw_shapes`). At first, it might look strange that parts of many particles lie
 outside the FMM group; however, the FMM is used only after the particles are
 converted to line sources, and are thus fully contained in the FMM grid.
 

docs/src/tutorial_optim_angle.md

@@ -0,0 +1,133 @@
+#  Tutorial 3: Angle Optimization
+
+In this tutorial, we build upon the [previous tutorial]@(ref tutorial2) by
+optimizing the rotation angles of the particles (`φs`) to maximize the field
+intensity at a specific point.
+Depending on the scattering problem, wavelengths, and incident field,
+optimization can have a major or minor impact on the field.
+The simplest way to perform this optimization is by calling [`optimize_φ`](@ref), which
+in turn utilizes `Optim`, a Julia package for nonlinear optimization. The type of
+objective function handled here is given by
+
+```math
+f_{\mathrm{obj}} = \sum_{\mathbf{r} \in I} \| u(\mathbf{r})\|^2
+```
+
+where ``I`` is a set of points that lie outside all scattering discs
+and ``u`` is the ``z``-component of the electric field. This function can be
+minimized directly, or maximized by minimizing ``-f_{\mathrm{obj}}``.
+
+First we set up our scattering problem:
+
+```julia
+λ0 = 1 #doesn't matter since everything is normalized to λ0
+k0 = 2π/λ0
+kin = 0.5k0
+θ_i = 0 #incident wave e^{i k_0 (1/sqrt{2},1/sqrt{2}) \cdot \mathbf{r}}
+M = 20
+shapes = [rounded_star(0.35λ0, 0.1λ0, 4, 202)]
+P = 12
+centers =  rect_grid(2, div(M,2), λ0, λ0) #2xM/2 grid
+ids = ones(Int, M)
+φs0 = zeros(M)
+sp = ScatteringProblem(shapes, ids, centers, φs0)
+fmm_options = FMMoptions(true, acc = 6, dx = 2λ0)
+points = 0.05*λ0*[-1 0; 1 0; 0 1; 0 -1]
+```
+
+where `φs0` is the starting point for the optimization method, and `points` are
+the locations at which we intend to maximize or minimize the field intensity.
+In this case, we want to optimize intensity at a small area around the origin.
+We now select the optimization method and select its options.
+In most cases, this combination of BFGS with a backtracking line search will
+yield accurate results in fast time;
+other line searches that require re-evaluation of the gradient will be
+significantly slower but may converge more accurately. The possible convergence
+criteria are set by the bounds `f_tol`, `g_tol`, and `x_tol`, for
+a relative change in the function, gradient norm, or variables, respectively.
+In addition, we can set a maximum number of `iterations`. Verbosity of the output
+is set with `show_trace` and `extended_trace`.
+
+```julia
+optim_options = Optim.Options(f_tol = 1e-6, iterations = 100,
+                    store_trace = true, show_trace = true)
+optim_method = Optim.BFGS(;linesearch = LineSearches.BackTracking())
+```
+
+We now run both minimization and maximization:
+
+```julia
+res_min = optimize_φ(φs0, points, P, θ_i, k0, kin, shapes, centers, ids,
+            fmm_options, optim_options, optim_method; minimize = true)
+res_max = optimize_φ(φs0, points, P, θ_i, k0, kin, shapes, centers, ids,
+            fmm_options, optim_options, optim_method; minimize = false)
+sp_min = ScatteringProblem(shapes, ids, centers, res_min.minimizer)
+sp_max = ScatteringProblem(shapes, ids, centers, res_max.minimizer)
+```
+
+Once the optimization is done, we can visualize each `ScatteringProblem`
+separately with [`plot_near_field`](@ref) or compare them side by side with the
+following PyPlot code:
+
+```julia
+plts = Array{Any}(3)
+plts[1] = plot_near_field(k0, kin, P, sp, θ_i, x_points = 100, y_points = 300,
+        opt = fmm_options, border = find_border(sp, points))
+plts[2] = plot_near_field(k0, kin, P, sp_min, θ_i, x_points = 100, y_points = 300,
+        opt = fmm_options, border = find_border(sp, points))
+plts[3] = plot_near_field(k0, kin, P, sp_max, θ_i, x_points = 100, y_points = 300,
+        opt = fmm_options, border = find_border(sp, points))
+close("all")
+
+fig, axs = subplots(ncols=3); msh = 0
+for (i, spi) in enumerate([sp;sp_min;sp_max])
+    msh = axs[i][:pcolormesh](plts[i][2][1], plts[i][2][2], abs.(plts[i][2][3]),
+                        vmin = 0, vmax = 3.4, cmap="viridis")
+    draw_shapes(spi.shapes, spi.centers, spi.ids, spi.φs, axs[i])
+    axs[i][:set_aspect]("equal", adjustable = "box")
+    axs[i][:set_xlim]([border[1];border[2]])
+    axs[i][:set_ylim]([border[3];border[4]])
+    axs[i][:set_xticks]([-1,0,1])
+    i > 1 && axs[i][:set_yticks]([])
+end
+subplots_adjust(left=0.05, right=0.8, top=0.98, bottom = 0.05, wspace = 0.1)
+cbar_ax = fig[:add_axes]([0.85, 0.05, 0.05, 0.93])
+fig[:colorbar](msh, cax=cbar_ax)
+```
+```@raw html
+<div style="text-align:center">
+<img alt=optim_angle src="./assets/optim_angle.png" style="width:80%; height:auto; margin:1%; max-width: 600px">
+</div><p style="clear:both;">
+```
+
+From left to right, we see the electric field before optimization, after
+minimization, and after maximization. The field intensity at the origin is
+notably different in both optimization results, with minimization decreasing
+the intensity by 95%, and maximization increasing it by over 700%. The convergence
+of the optimization method for both examples can be plotted via:
+
+```julia
+iters = length(res_min.trace)
+fobj = [res_min.trace[i].value for i=1:iters]
+gobj = [res_min.trace[i].g_norm for i=1:iters]
+iters2 = length(res_max.trace)
+fobj2 = -[res_max.trace[i].value for i=1:iters2]
+gobj2 = [res_max.trace[i].g_norm for i=1:iters2]
+
+fig, axs = subplots(ncols=2, figsize=[7,5])
+axs[1][:semilogy](0:iters-1, fobj, linewidth=2)
+axs[2][:semilogy](0:iters-1, gobj, linewidth=2)
+axs[1][:semilogy](0:iters2-1, fobj2, linewidth=2, "--")
+axs[2][:semilogy](0:iters2-1, gobj2, linewidth=2, "--")
+axs[1][:legend](["\$f_\\mathrm{obj}\$ (min)";
+                "\$f_\\mathrm{obj}\$ (max)"], loc="right")
+axs[2][:legend](["\$\\|\\mathbf{g}_\\mathrm{obj}\\|\$ (min)";
+                "\$\\|\\mathbf{g}_\\mathrm{obj}\\|\$ (max)"], loc="best")
+axs[1][:set_xlabel]("Iteration")
+axs[2][:set_xlabel]("Iteration")
+axs[1][:set_ylim](ymax=40)
+```
+
+```@raw html
+<p style="text-align:center;"><img alt=optim_angle_conv src="./assets/optim_angle_conv.png" style="width:70%; height:auto; max-width:600px"></p>
+```

docs/src/tutorial_optim_radius.md

@@ -85,10 +85,12 @@ xlabel("x/a")
 ylabel("y/a")
 ```
 
+```@raw html
 <div style="text-align:center">
 <img alt=optim_radius_before src="./assets/optim_radius_before.png" style="width:40%; height:auto; margin:1%; max-width: 300px">
 <img alt=optim_radius_after src="./assets/optim_radius_after.png" style="width:40%; height:auto; margin:1%; max-width: 300px">
 </div><p style="clear:both;">
+```
 
 `res` also stores the objective value as well as the g
 radient norm in each iteration.
@@ -103,6 +105,9 @@ rng = iters .== 0
 ```
 where `rng` now contains the indices at which a new outer iteration has begun.
 Finally, plotting `fobj` and `gobj` for this example yields the following plot:
+
+```@raw html
 <p style="text-align:center;"><img alt=optim_radius_conv src="./assets/optim_radius_conv.png" style="width:60%; height:auto; max-width:400px"></p>
+```
 
 where markers denote the start of an outer iteration.

src/fmm_matrices.jl

@@ -136,12 +136,12 @@ function FMMnearMatrix(k, P, groups, centers, boxSize, num)
     Is = Array{Int}(num*W^2)
     Js = Array{Int}(num*W^2)
     Zs = Array{Complex{Float64}}(num*W^2)
+    bess = Array{Complex{Float64}}(2*P+1)
     mindist2 = 3*boxSize^2 #anywhere between 2 and 4
     ind = 0
     for ig1 = 1:G, ig2 = 1:G
         x = groups[ig1].center - groups[ig2].center
         sum(abs2,x) > mindist2 && continue
-
         #for each enclosed scatterer, build translation matrix
         for ic1 = 1:length(groups[ig1].point_ids)
             for ic2 = 1:length(groups[ig2].point_ids)

src/optimize_phis.jl

@@ -1,7 +1,6 @@
 """
     optimize_φ(φs0, points, P, θ_i, k0, kin, shapes, centers, ids, fmmopts,
-        optimopts::Optim.Options, method::Optim.AbstractOptimizer;
-        minimize = true)
+        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`.
@@ -10,10 +9,10 @@ and other optimization parameters.
 Returns an object of type `Optim.MultivariateOptimizationResults`.
 """
 function optimize_φ(φs0, points, P, θ_i, k0, kin, shapes, centers, ids, fmmopts,
-                    optimopts::Optim.Options, method::Optim.AbstractOptimizer;
-                    minimize = true)
+                    optimopts::Optim.Options, method; minimize = true)
 
     #stuff that is done once
+    verify_min_distance(shapes, centers, ids, points) || error("Particles are too close.")
     mFMM, scatteringMatrices, scatteringLU, buf =
         prepare_fmm_reusal_φs(k0, kin, P, shapes, centers, ids, fmmopts)
     Ns = size(centers,1)

src/optimize_rs.jl

@@ -11,9 +11,9 @@ Here, `ids` allows for grouping particles - for example, to maintain symmetry of
 the optimized device.
 `optimopts` defines the convergence criteria and other optimization parameters
 for both the inner and outer iterations. `method` can be either `"BFGS"` or
-`"LBFGS"`. See the `Fminbox` documentation for more details.
+`"LBFGS"`. See the `Optim.Fminbox` documentation for more details.
 
-Returns an object of type ??????
+Returns an object of type `Optim.MultivariateOptimizationResults`.
 """
 function optimize_radius(rs0, r_min, r_max, points, ids, P, θ_i, k0, kin,
                         centers, fmmopts, optimopts::Optim.Options;
@@ -32,6 +32,8 @@ function optimize_radius(rs0, r_min, r_max, points, ids, P, θ_i, k0, kin,
     else
         assert(J == length(r_max))
     end
+    verify_min_distance([CircleParams(r_max[i]) for i = 1:J], centers, ids,
+        points) || error("Particles are too close.")
 
     #setup FMM reusal
     groups, boxSize = divideSpace(centers, fmmopts)

src/shapes.jl

@@ -169,6 +169,7 @@ end
 
 """
     verify_min_distance(shapes, centers::Array{Float64,2}, ids)
+    verify_min_distance(sp::ScatteringProblem)
 Returns `true` if the shapes placed at `centers` are properly distanced
 (non-intersecting scattering disks).
 """
@@ -184,6 +185,7 @@ end
 
 """
     verify_min_distance(shapes, centers::Array{Float64,2}, ids, points::Array{Float64,2})
+    verify_min_distance(sp::ScatteringProblem, points)
 Returns `true` if the shapes placed at `centers` are properly distanced
 (non-intersecting scattering disks), and all `points` are outside the scattering
 disks.

test/fmm_test.jl

@@ -0,0 +1,48 @@
+@testset "FMM" begin
+    λ0 = 1
+    k0 = 2π/λ0
+    kin = 3k0
+    θ_i = π/4 #incident wave e^{i k_0 (1/sqrt{2},1/sqrt{2}) \cdot \mathbf{r}}
+    N_squircle = 200
+    N_star = 210
+    P = 10
+
+    M = 10
+    shapes = [rounded_star(0.1λ0, 0.03λ0, 5, N_star);
+                squircle(0.15λ0, N_squircle)]
+    ids = rand(1:2, M^2)
+    φ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],
+                failures = 1_000)
+    sp = ScatteringProblem(shapes, ids, centers, φs)
+    @test verify_min_distance(sp, [0.0 0.0])
+    catch
+    warn("Could not find random points (1)")
+    end
+    try
+    centers = randpoints(M^2, dist, 5λ0, 5λ0, failures = 1_000)
+    sp = ScatteringProblem(shapes, ids, centers, φs)
+    @test verify_min_distance(sp)
+    catch
+    warn("Could not find random points (2)")
+    end
+
+    centers =  rect_grid(M, M, 0.8λ0, λ0)
+    sp = ScatteringProblem(shapes, ids, centers, φs)
+    fmm_options = FMMoptions(true, acc = 6, nx = div(M,2))
+
+    groups,boxSize = divideSpace(centers, fmm_options; drawGroups = false)
+    P2,Q = FMMtruncation(6, boxSize, k0)
+    mFMM1 = FMMbuildMatrices(k0, P, P2, Q, groups, centers, boxSize; tri = false)
+    mFMM2 = FMMbuildMatrices(k0, P, P2, Q, groups, centers, boxSize; tri = true)
+    @test norm(mFMM1.Znear - mFMM2.Znear, Inf) == 0
+
+    points = [linspace(-λ0*M, λ0*M, 200) zeros(200)]
+    u1 = calc_near_field(k0, kin, P, sp, points, θ_i; opt = fmm_options,
+            verbose = false)
+    u2 = calc_near_field(k0, kin, P, sp, points, θ_i; opt = fmm_options,
+            verbose = true, use_multipole = false)
+    @test norm(u1 - u2)/norm(u1) < 1e-6
+end

test/optimize_angle_test.jl

@@ -0,0 +1,55 @@
+#based on the optimizing angle tutorial
+import Optim
+
+@testset "optimize angle" begin
+    λ0 = 1 #doesn't matter since everything is normalized to λ0
+    k0 = 2π/λ0
+    kin = 0.5k0
+    θ_i = 0.0
+    M = 20
+
+    sfun(N) = rounded_star(0.35λ0, 0.1λ0, 4, N)
+    # N, errN = minimumN(k0, kin, sfun; tol=1e-6, N_start = 200, N_min = 100,
+                # N_max = 400)
+    shapes = [sfun(202)]
+    P = 12#P, errP = minimumP(k0, kin, shapes[1]; P_max = 30, tol = 1e-7)
+
+    centers =  rect_grid(2, div(M,2), λ0, λ0) #2xM/2 grid with distance λ0
+    ids = ones(Int, M)
+    φs0 = zeros(M)
+    sp = ScatteringProblem(shapes, ids, centers, φs0)
+
+    points = [-0.05 0.0; 0.05 0.0; 0.0 0.05; 0.0 -0.05]
+    assert(verify_min_distance(sp, points))
+
+    fmm_options = FMMoptions(true, acc = 6, dx = 2λ0)
+
+    optim_options = Optim.Options(f_tol = 1e-6, iterations = 50,
+                        store_trace = true, show_trace = false)
+    optim_method = Optim.BFGS(;linesearch = LineSearches.BackTracking())
+    res_min = optimize_φ(φs0, points, P, θ_i, k0, kin, shapes, centers, ids, fmm_options,
+            optim_options, optim_method; minimize = true)
+    @test res_min.f_converged
+    res_max = optimize_φ(φs0, points, P, θ_i, k0, kin, shapes, centers, ids, fmm_options,
+            optim_options, optim_method; minimize = false)
+    @test res_max.f_converged
+    sp_min = ScatteringProblem(shapes, ids, centers, res_min.minimizer)
+    sp_max = ScatteringProblem(shapes, ids, centers, res_max.minimizer)
+
+    u_bef = calc_near_field(k0, kin, P, sp, points, θ_i; opt = fmm_options,
+            verbose = false)
+    u_min = calc_near_field(k0, kin, P, sp_min, points, θ_i; opt = fmm_options,
+            verbose = false)
+    u_max = calc_near_field(k0, kin, P, sp_max, points, θ_i; opt = fmm_options,
+            verbose = false)
+
+    fobj0 = sum(abs2.(u_bef))
+    fobj_min = sum(abs2.(u_min))
+    fobj_max = sum(abs2.(u_max))
+
+    @test isapprox(res_min.trace[1].value, fobj0)
+    @test isapprox(res_min.minimum, fobj_min)
+    @test isapprox(-res_max.minimum, fobj_max)
+
+    border = find_border(sp, points)
+end