Merge pull request #24 from bblankrot/dev

Added adjoint-state optimization and power optimization

src/PS_types.jl

@@ -191,7 +191,7 @@ end
 """
 	LineSource(x, y)
 
-Constructor for the `LineSource` type, where `(x0,y0)` is the coordinate of the
+Constructor for the `LineSource` type, where `(x,y)` is the coordinate of the
 current filament.
 """
 struct LineSource <: Einc

src/ParticleScattering.jl

@@ -51,16 +51,22 @@ export optimize_φ
 #methods, optimize_rs.jl
 export optimize_radius
 #methods, incident.jl
-export uinc
-
-# For advanced plotting with pgfplots
-import DataFrames, CSV, PGFPlotsX; const pgf = PGFPlotsX
-include("visualization_pgf.jl")
-#methods, visualization_pgf.jl
-export plot_near_field_pgf
+export uinc, hxinc, hyinc
+#consts, incident.jl
+export eta0
 
 #temp
 export divideSpace, FMMtruncation, particleExpansion, FMMbuildMatrices
 export FMMbuffer, FMMmatrices
 
+include("optimize_rs_mf.jl")
+export optimize_rs_mf
+
+include("poynting.jl")
+export calc_power
+
+include("optimize_rs_pwr.jl")
+export PowerBuffer, OptimProblemBuffer, optMatrixPwr, optimize_pwr_rs_f, optimize_pwr_rs_g!
+include("optimize_phis_pwr.jl")
+export optimize_pwr_φ_f, optimize_pwr_φ_g!
 end

src/fmm_mvp.jl

@@ -137,3 +137,72 @@ function FMM_mainMVP_pre2!(output, beta, scatteringMatrices, φs::Vector{Float64
     end #inbounds
     return output
 end
+
+function FMM_mainMVP_transpose!(output, beta, scatteringMatrices, φs::Vector{Float64}, ids::Vector{Int}, P, mFMM, pre_agg, translated_sum)
+    #calculate matrix^T-vector product - devectorized with pre-preconditioning
+
+    @inbounds begin
+    #here we first multiply by Xᵀ
+
+    #Xβ storage can be moved outward, but for now this preallocation prevents most
+    #dynamic mem alloc by this MVP
+    Xβ = Array{Complex{Float64}}(length(beta))
+    temp = Array{Complex{Float64}}(2*P+1)
+    for ic = 1:length(ids)
+        rng = (ic-1)*(2*P+1) + (1:2*P+1)
+        v = view(Xβ, rng)
+        if φs[ic] == 0.0
+            At_mul_B!(v, scatteringMatrices[ids[ic]], view(beta, rng))
+        else
+            #rotate without matrix - transposed
+            rotateMultipole!(temp, view(beta, rng), φs[ic], P)
+            At_mul_B!(v, scatteringMatrices[ids[ic]], temp)
+            rotateMultipole!(v, -φs[ic], P)
+        end
+    end
+
+    At_mul_B!(output, mFMM.Znear, Xβ)
+    G = length(mFMM.groups)
+    fill!(pre_agg,0.0)
+    #first preagg all
+    for ig2 = 1:G
+        for is = 1:mFMM.groups[ig2].size
+            indices = (mFMM.groups[ig2].point_ids[is]-1)*(2*P+1)
+            for ii = 1:2*P+1
+                for iq = 1:mFMM.Q
+                    pre_agg[iq,ig2] += conj(mFMM.Agg[ig2][iq,(is-1)*(2*P+1) + ii])*Xβ[indices + ii]
+                end
+            end
+        end
+    end
+
+    for ig1 = 1:G
+        #translate plane waves from ig2 to ig1
+        fill!(translated_sum,0.0)
+        for ig2 = 1:G
+            if isempty(mFMM.Trans[(ig2-1)*G + ig1]) #transposed
+                continue
+            else
+                for iQ = 1:mFMM.Q
+                    translated_sum[iQ] -= mFMM.Trans[(ig2-1)*G + ig1][iQ]*pre_agg[iQ,ig2] #transposed
+                end
+                #minus sign because the real equation is (I-TᵀXᵀ)β=b
+            end
+        end
+        #disaggregate from ig1 center to ig1's scatterers
+        for is = 1:mFMM.groups[ig1].size
+            for ip = 1:2*P+1
+                disagged = 0.0im
+                for iq = 1:mFMM.Q
+                    disagged += mFMM.Agg[ig1][iq,(is-1)*(2*P+1) + ip]*translated_sum[iq]
+                end
+                output[(mFMM.groups[ig1].point_ids[is]-1)*(2*P+1) + ip] += disagged
+            end
+        end
+    end
+
+    #add identity matrix (Iβ)
+    output .+= beta
+    end #inbounds
+    return output
+end

src/incident.jl

@@ -1,3 +1,5 @@
+const eta0 = 4π*299792458e-7
+
 #plane wave
 function u2α(k, u::PlaneWave, centers::Array{T,2}, P) where T <: Real
 	#build incoming coefficients for plane wave incident field
@@ -57,7 +59,7 @@ function uinc(k, points::Array{T,1}, u::LineSource) where T <: Real
 	0.25im*besselh(0, 1, k*r)
 end
 
-#current source
+#current source TODO: simple avg->trap rule
 function u2α(k, u::CurrentSource, centers::Array{T,2}, P) where T <: Real
 	α = zeros(Complex{Float64}, size(centers,1)*(2P + 1))
 	c = 0.25im*u.len/length(u.σ)
@@ -123,3 +125,62 @@ function err_α(k, ui::Einc, P, shape, center)
 	u_α = sum(α[p + P + 1]*besselj(p, k*shape.R)*exp.(1im*p*θ) for p = -P:P)
 	norm(u_α - u_exact)/norm(u_exact)
 end
+
+
+function hxinc(k, p, u::PlaneWave)
+    (sin(u.θi)/eta0)*exp.(1.0im*k*(cos(u.θi)*p[:,1] + sin(u.θi)*p[:,2]))
+end
+
+function hyinc(k, p, u::PlaneWave)
+    (-cos(u.θi)/eta0)*exp.(1.0im*k*(cos(u.θi)*p[:,1] + sin(u.θi)*p[:,2]))
+end
+
+function hxinc(k, p, u::LineSource)
+    R = sqrt.((view(p,:,1) - u.x).^2 + (view(p,:,2) - u.y).^2)
+    if any(R .== 0)
+		warn("Hinc: encountered singularity in incident field, returned NaN")
+		R[R.==0] = NaN
+	end
+    h = (0.25/eta0./R).*besselh.(1,k*R).*(u.y - view(p,:,2))
+end
+
+function hyinc(k, p, u::LineSource)
+    R = sqrt.((view(p,:,1) - u.x).^2 + (view(p,:,2) - u.y).^2)
+    if any(R .== 0)
+		warn("Hinc: encountered singularity in incident field, returned NaN")
+		R[R.==0] = NaN
+	end
+    h = (0.25/eta0./R).*besselh.(1,k*R).*(view(p,:,1) - u.x)
+end
+
+#these untested
+function hxinc(k, points, u::CurrentSource)
+    res = Array{Complex{Float64}}(size(points, 1))
+	r = Array{Float64}(size(u.p, 1))
+    y = Array{Float64}(size(u.p, 1))
+	for i = 1:size(points, 1)
+		y .= points[i,2] - u.p[:,2]
+        r .= hypot.(points[i,1] - u.p[:,1], y)
+		if any(r .== 0) #point is on source
+			warn("Hinc: encountered singularity in incident field, returned NaN")
+			res[i] = NaN
+		end
+		res[i] = (-0.25*u.len/length(u.σ)/eta0)*sum(besselh.(1, 1, k*r).*u.σ.*y./r)
+	end
+	res
+end
+function hyinc(k, points, u::CurrentSource)
+    res = Array{Complex{Float64}}(size(points, 1))
+	r = Array{Float64}(size(u.p, 1))
+    x = Array{Float64}(size(u.p, 1))
+	for i = 1:size(points, 1)
+		x .= points[i,1] - u.p[:,1]
+        r .= hypot.(x, points[i,2] - u.p[:,2])
+		if any(r .== 0) #point is on source
+			warn("Hinc: encountered singularity in incident field, returned NaN")
+			res[i] = NaN
+		end
+		res[i] = (0.25*u.len/length(u.σ)/eta0)*sum(besselh.(1, 1, k*r).*u.σ.*x./r)
+	end
+	res
+end

src/multipole.jl

@@ -96,13 +96,17 @@ function M2Lmatrix!(T, k, P, d)
 	end
 end
 
-function scattered_field_multipole(k0, beta, centers::Array{Float64,2}, points::Array{Float64,2})
+function scattered_field_multipole(k0, beta, centers::Array{Float64,2}, points::Array{Float64,2}; recurrence = false)
 	Ns = size(centers,1)
 	P = div(div(length(beta),Ns)-1,2)
 	len = size(points,1)
 	Esc = zeros(Complex{Float64}, len)
 
-	scattered_field_multipole!(Esc, k0, beta, P, centers, 1:Ns, points, 1:len)
+	if recurrence
+		scattered_field_multipole_recurrence!(Esc, k0, beta, P, centers, 1:Ns, points, 1:len)
+	else
+		scattered_field_multipole!(Esc, k0, beta, P, centers, 1:Ns, points, 1:len)
+	end
 	return Esc
 end
 
@@ -115,18 +119,41 @@ function scattered_field_multipole!(Esc::Array{Complex{Float64},1}, k0, beta, P,
 
 			rs_moved = hypot(points_moved1, points_moved2)
 			ts_moved = atan2(points_moved2, points_moved1)
-			try
-				Esc[ip] += beta[ind]*besselh(0, 1, k0*rs_moved)
-			catch
+			if rs_moved == 0
 				error("rs_moved == 0, center=$(centers[ic,:]), point=$(points[ip,:]), k0 = $k0, ic=$ic, ip=$ip")
 			end
+			Esc[ip] += beta[ind]*besselh(0, 1, k0*rs_moved)
 			for p = 1:P
 				Esc[ip] += besselh(p, 1, k0*rs_moved)*(beta[p + ind]*exp(1im*p*ts_moved) + (-1)^p*beta[-p + ind]*exp(-1im*p*ts_moved))
 			end
 		end
 	end
 end
 
+function scattered_field_multipole_recurrence!(Esc::Array{Complex{Float64},1}, k0, beta, P, centers::Array{Float64,2}, ind_centers, points::Array{Float64,2}, ind_points)
+	for ic in ind_centers
+		ind = (ic-1)*(2*P+1) + P + 1
+		for ip in ind_points
+			points_moved1 = points[ip,1] - centers[ic,1]
+			points_moved2 = points[ip,2] - centers[ic,2]
+
+			rs_moved = hypot(points_moved1, points_moved2)
+			ts_moved = atan2(points_moved2, points_moved1)
+			if rs_moved == 0
+				error("rs_moved == 0, center=$(centers[ic,:]), point=$(points[ip,:]), k0 = $k0, ic=$ic, ip=$ip")
+			end
+			Hₚ₋₂ = besselh(-1, 1, k0*rs_moved)
+			Hₚ₋₁ = besselh(0, 1, k0*rs_moved)
+			Esc[ip] += beta[ind]*Hₚ₋₁
+			for p = 1:P
+				Hₚ = (2*(p-1)/(k0*rs_moved))*Hₚ₋₁ - Hₚ₋₂
+				Esc[ip] += Hₚ*(beta[p + ind]*exp(1im*p*ts_moved) + (-1)^p*beta[-p + ind]*exp(-1im*p*ts_moved))
+				Hₚ₋₂ = Hₚ₋₁; Hₚ₋₁ = Hₚ
+			end
+		end
+	end
+end
+
 function circleScatteringMatrix(kout, kin, R, P; gamma = false)
     #non-vectorized, reuses bessel
     S = Array{Complex{Float64}}(2*P+1)

src/optimize_phis.jl

@@ -6,10 +6,12 @@ Optimize the rotation angles of a particle collection for minimization or
 maximization (depending on `minimize`) of the field intensity at `points`.
 `optimopts` and `method` define the optimization method, convergence criteria,
 and other optimization parameters.
+`adjoint` dictates whether the gradient is calculated using the adjoint method
+(faster) or directly.
 Returns an object of type `Optim.MultivariateOptimizationResults`.
 """
 function optimize_φ(φs0, points, P, ui::Einc, k0, kin, shapes, centers, ids, fmmopts,
-                    optimopts::Optim.Options, method; minimize = true)
+                    optimopts::Optim.Options, method; minimize = true, adjoint = true)
 
     #stuff that is done once
     verify_min_distance(shapes, centers, ids, points) || error("Particles are too close.")
@@ -25,27 +27,22 @@ function optimize_φ(φs0, points, P, ui::Einc, k0, kin, shapes, centers, ids, f
     last_φs = similar(initial_φs)
     last_φs == initial_φs && (last_φs[1] += 1) #should never happen but has
 
-    if minimize
-        df = OnceDifferentiable(
-            φs -> optimize_φ_f(φs, shared_var, last_φs, α,
-                                    H, points, P, uinc_, Ns, k0, centers,
-                                    scatteringMatrices, ids, mFMM, fmmopts, buf),
-            (grad_stor, φs) -> optimize_φ_g!(grad_stor, φs, shared_var,
-                                    last_φs, α, H, points, P, uinc_,
-                                    Ns, k0, centers, scatteringMatrices,
-                                    scatteringLU, ids, mFMM, fmmopts, buf, minimize),
-                                    initial_φs)
+    fopt = φs -> ifelse(minimize,1,-1)*optimize_φ_f(φs, shared_var, last_φs, α,
+                                        H, points, P, uinc_, Ns, k0, centers,
+                                        scatteringMatrices, ids, mFMM, fmmopts,
+                                        buf)
+    if adjoint
+        gopt! = (grad_stor, φs) -> optimize_φ_adj_g!(grad_stor, φs, shared_var,
+                                    last_φs, α, H, points, P, uinc_, Ns, k0,
+                                    centers, scatteringMatrices, scatteringLU,
+                                    ids, mFMM, fmmopts, buf, minimize)
     else
-        df = OnceDifferentiable(
-            φs -> -optimize_φ_f(φs, shared_var, last_φs, α,
-                                    H, points, P, uinc_, Ns, k0, centers,
-                                    scatteringMatrices, ids, mFMM, fmmopts, buf),
-            (grad_stor, φs) -> optimize_φ_g!(grad_stor, φs, shared_var,
-                                    last_φs, α, H, points, P, uinc_,
-                                    Ns, k0, centers, scatteringMatrices,
-                                    scatteringLU, ids, mFMM, fmmopts, buf, minimize),
-                                    initial_φs)
+        gopt! = (grad_stor, φs) -> optimize_φ_g!(grad_stor, φs, shared_var,
+                                    last_φs, α, H, points, P, uinc_, Ns, k0,
+                                    centers, scatteringMatrices, scatteringLU,
+                                    ids, mFMM, fmmopts, buf, minimize)
     end
+    df = OnceDifferentiable(fopt, gopt!, initial_φs)
     optimize(df, initial_φs, method, optimopts)
 end
 
@@ -149,6 +146,49 @@ function optimize_φ_g!(grad_stor, φs, shared_var, last_φs, α, H, points, P,
                     real(shared_var.∂β.'*(H*conj(shared_var.f)))
 end
 
+
+function optimize_φ_adj_g!(grad_stor, φs, shared_var, last_φs, α, H, points, P, uinc_, Ns, k0, centers, scatteringMatrices, scatteringLU, ids, mFMM, opt, buf, minimize)
+    optimize_φ_common!(φs, last_φs, shared_var, α, H, points, P, uinc_, Ns,
+        k0, centers,scatteringMatrices, ids, mFMM, opt, buf)
+
+    MVP = LinearMap{eltype(buf.rhs)}(
+            (output_, x_) -> FMM_mainMVP_transpose!(output_, x_,
+                                scatteringMatrices, φs, ids, P, mFMM,
+                                buf.pre_agg, buf.trans),
+            Ns*(2*P+1), Ns*(2*P+1), ismutating = true)
+
+    #build rhs of adjoint problem
+    shared_var.rhs_grad[:] = -(H*conj(shared_var.f))
+    #solve adjoint problem
+    λadj, ch = gmres(MVP, shared_var.rhs_grad, restart = Ns*(2*P+1),
+                    tol = opt.tol, log = true)
+    # λadj, ch = gmres!(λadj, MVP, shared_var.rhs_grad, restart = Ns*(2*P+1),
+    #                 tol = opt.tol, log = true, initially_zero = true)
+    if ch.isconverged == false
+        display("FMM process did not converge for adjoint system.")
+        error("..")
+    end
+
+    #note: this assumes that there are exactly 2P+1 non zero elements in ∂X. If
+    #not, v must be (2P+1)Ns × 1.
+    D = -1.0im*collect(-P:P)
+    v = Array{Complex{Float64}}(2P+1) #TODO: minimize dynamic alloc
+    for n = 1:Ns
+        #compute n-th element of gradient
+        rng = (n-1)*(2*P+1) + (1:2*P+1)
+        rotateMultipole!(v, view(shared_var.β,rng), -φs[n], P)
+        v[:] = scatteringLU[ids[n]]\v #LU decomp with pivoting
+        v[:] .*= -D
+        v[:] = scatteringMatrices[ids[n]]*v
+        rotateMultipole!(v, φs[n], P)
+        v[:] += D.*shared_var.β[rng]
+
+        grad_stor[n] = ifelse(minimize, -2, 2)*real(view(λadj,rng).'*v)
+        #prepare for next one - #TODO: check why this is here
+        v[:] = 0.0im
+    end
+end
+
 function optimizationHmatrix(points, centers, Ns, P, k0)
     points_moved = Array{Float64}(2)
     H = Array{Complex{Float64}}(Ns*(2*P+1), size(points,1))