Add options and plan to Afft and Asense (#124)

* Add return x,y

* for in

* Use plan in Afft

* Test new Afft

* Add fft_forward option

* Use plan

* Test plan

* Plain Afft

* Refine Afft

* Test Afft

* Add unitary option

* Test unitary Asense

* v0.17

* Test fft_forward

* Check 1-based

Author: JeffFessler

Project.toml

@@ -1,7 +1,7 @@
 name = "MIRT"
 uuid = "7035ae7a-3787-11e9-139a-5545ed3dc201"
 authors = ["fessler <[email protected]>"]
-version = "0.16.0"
+version = "0.17.0"
 
 [deps]
 AVSfldIO = "b6189060-daf9-4c28-845a-cc0984b81781"

src/algorithm/general/ncg.jl

@@ -27,9 +27,9 @@ where ``C_j(u)`` is diagonal matrix of curvatures.
 This CG method uses a majorize-minimize (MM) line search.
 
 in
-- `B`		array of ``J`` blocks ``B_1,...,B_J``
-- `gradf`	array of ``J`` functions return gradients of ``f_1,...,f_J``
-- `curvf`	array of ``J`` functions `z -> curv(z)` that return a scalar
+- `B`		vector of ``J`` blocks ``B_1,...,B_J``
+- `gradf`	vector of ``J`` functions return gradients of ``f_1,...,f_J``
+- `curvf`	vector of ``J`` functions `z -> curv(z)` that return a scalar
 or a vector of curvature values for each element of ``z``
 - `x0`	initial guess; need `length(x) == size(B[j],2)` for ``j=1...J``
 
@@ -58,10 +58,12 @@ function ncg(
     niter::Int = 50,
     ninner::Int = 5,
     P = I, # trick: this is an overloaded I (by LinearMapsAA)
-    betahow::Symbol=:dai_yuan,
+    betahow::Symbol = :dai_yuan,
     fun::Function = (x,iter) -> 0,
 )
 
+    Base.require_one_based_indexing(B, gradf, curvf)
+
 	out = Array{Any}(undef, niter+1)
 	out[1] = fun(x0, 0)
 
@@ -72,10 +74,10 @@ function ncg(
 	grad_old = []
 	grad_new = []
 
-	Bx = [B[j] * x for j=1:J] # u_j in course notes
-	grad = (Bx) -> sum([B[j]' * gradf[j](Bx[j]) for j=1:J])
+	Bx = [B[j] * x for j in 1:J] # u_j in course notes
+	grad = (Bx) -> sum([B[j]' * gradf[j](Bx[j]) for j in 1:J])
 
-for iter = 1:niter
+for iter in 1:niter
 	grad_new = grad(Bx) # gradient
 	npgrad = -(P * grad_new)
 	if iter == 1
@@ -97,14 +99,14 @@ for iter = 1:niter
 
 	# MM-based line search for step size alpha
 	# using h(a) = sum_j f_j(uj + a vj)
-	Bd = [B[j] * dir for j=1:J] # v_j in course notes
+	Bd = [B[j] * dir for j in 1:J] # v_j in course notes
 
 	alf = 0
-	for ii=1:ninner
-	#	derh = alf -> sum([Bd[j]' * gradf[j](Bx[j] + alf * Bd[j]) for j=1:J])
+	for ii in 1:ninner
+	#	derh = alf -> sum([Bd[j]' * gradf[j](Bx[j] + alf * Bd[j]) for j in 1:J])
 		derh = 0 # derivative of h(a)
 		curv = 0
-		for j=1:J
+		for j in 1:J
 			tmp = Bx[j] + alf * Bd[j]
 			derh += real(dot(Bd[j], gradf[j](tmp)))
 			curv += sum(curvf[j](tmp) .* abs2.(Bd[j]))
@@ -119,7 +121,7 @@ for iter = 1:niter
 	end
 
 	x += alf * dir
-	for j=1:J # update Bj * x
+	for j in 1:J # update Bj * x
 		Bx[j] += alf * Bd[j]
 	end
 	out[iter+1] = fun(x, iter)

src/algorithm/general/ogm_ls.jl

@@ -24,23 +24,23 @@ where ``C_j(u)`` is diagonal matrix of curvatures.
 
 This OGM method uses a majorize-minimize (MM) line search.
 
-in
-- `B`		array of ``J`` blocks ``B_1,...,B_J``
-- `gradf`	array of ``J`` functions return gradients of ``f_1,...,f_J``
-- `curvf`	array of ``J`` functions `z -> curv(z)` that return a scalar
-or a vector of curvature values for each element of ``z``
-- `x0`	initial guess; need `length(x) == size(B[j],2)` for ``j=1...J``
-
-option
-- `niter`	# number of outer iterations; default 50
-- `ninner`	# number of inner iterations of MM line search; default 5
-- `fun`		User-defined function to be evaluated with two arguments (x,iter).
+# in
+- `B` vector of ``J`` blocks ``B_1,...,B_J``
+- `gradf` vector of ``J`` functions return gradients of ``f_1,...,f_J``
+- `curvf` vector of ``J`` functions `z -> curv(z)` that return a scalar
+  or a vector of curvature values for each element of ``z``
+- `x0` initial guess; need `length(x) == size(B[j],2)` for ``j=1...J``
+
+# option
+- `niter` # number of outer iterations; default 50
+- `ninner` # number of inner iterations of MM line search; default 5
+- `fun`	 User-defined function to be evaluated with two arguments (x,iter).
  * It is evaluated at (x0,0) and then after each iteration.
 
-output
-- `x`		final iterate
-- `out`		`[niter+1] (fun(x0,0), fun(x1,1), ..., fun(x_niter,niter))`
-  * (all 0 by default). This is an array of length `niter+1`
+# output
+- `x` final iterate
+- `out	(niter+1) (fun(x0,0), fun(x1,1), ..., fun(x_niter,niter))`
+  * (all 0 by default). This is a vector of length `niter+1`.
 """
 function ogm_ls(
     B::AbstractVector{<:Any},
@@ -52,6 +52,8 @@ function ogm_ls(
     fun::Function = (x,iter) -> 0,
 )
 
+    Base.require_one_based_indexing(B, gradf, curvf)
+
 out = Array{Any}(undef, niter+1)
 out[1] = fun(x0, 0)
 
@@ -65,12 +67,12 @@ grad_sum = zeros(size(x0))
 ti = 1
 thetai = 1
 
-B0 = [B[j] * x for j=1:J]
+B0 = [B[j] * x for j in 1:J]
 Bx = copy(B0)
 By = copy(B0)
-grad = (Bx) -> sum([B[j]' * gradf[j](Bx[j]) for j=1:J])
+grad = (Bx) -> sum([B[j]' * gradf[j](Bx[j]) for j in 1:J])
 
-for iter = 1:niter
+for iter in 1:niter
 	grad_new = grad(Bx) # gradient of x_{iter-1}
 	grad_sum += ti * grad_new # sum_{j=0}^{iter-1} t_j * gradient_j
 
@@ -80,21 +82,21 @@ for iter = 1:niter
 	tt = (iter < niter) ? ti : thetai # use theta_i factor for last iteration
 	yi = (1 - 1/tt) * x + (1/tt) * x0
 
-	for j=1:J # update Bj * yi
+	for j in 1:J # update Bj * yi
 		By[j] = (1 - 1/tt) * Bx[j] + (1/tt) * B0[j]
 	end
 
 	dir = -(1 - 1/tt) * grad_new - (2/tt) * grad_sum # -d_i
 
 	# MM-based line search for step size alpha
 	# using h(a) = sum_j f_j(By_j + a * Bd_j)
-	Bd = [B[j] * dir for j=1:J]
+	Bd = [B[j] * dir for j in 1:J]
 
 	alf = 0
-	for ii=1:ninner
+	for ii in 1:ninner
 		derh = 0 # derivative of h(a)
 		curv = 0
-		for j=1:J
+		for j in 1:J
 			tmp = By[j] + alf * Bd[j]
 			derh += real(dot(Bd[j], gradf[j](tmp)))
 			curv += sum(curvf[j](tmp) .* abs2.(Bd[j]))
@@ -109,20 +111,20 @@ for iter = 1:niter
 	end
 
 #	# derivative of h(a) = cost(x + a * dir) where \alpha is real
-#	dh = alf -> real(sum([Bd[j]' * gradf[j](By[j] + alf * Bd[j]) for j=1:J]))
-#	Ldh = sum([Lgf[j] * norm(Bd[j])^2 for j=1:J]) # Lipschitz constant for dh
+#	dh = alf -> real(sum([Bd[j]' * gradf[j](By[j] + alf * Bd[j]) for j in 1:J]))
+#	Ldh = sum([Lgf[j] * norm(Bd[j])^2 for j in 1:J]) # Lipschitz constant for dh
 #	(alf, ) = gd(dh, Ldh, 0, niter=ninner) # GD-based line search
 # todo
 
 	x = yi + alf * dir
 
 	if iter < niter
-		for j=1:J # update Bj * x
+		for j in 1:J # update Bj * x
 			Bx[j] = By[j] + alf * Bd[j]
 		end
 	end
 
-#	for j=1:J # recursive update Bj * yi ???
+#	for j in 1:J # recursive update Bj * yi ???
 #		By[j] = (1 - 1/ti) * (By[j] + alf * Bd[j]) + (1/ti) * B0[j]
 #	end
 
@@ -136,10 +138,11 @@ end
 """
     (x,out) = ogm_ls(grad, curv, x0, ...)
 
-special case of `ogm_ls` (OGM with line search) for minimizing a cost function
+Special case of `ogm_ls` (OGM with line search)
+for minimizing a cost function
 whose gradient is `grad(x)`
 and that has a quadratic majorizer with diagonal Hessian given by `curv(x)`.
-Typically `curv = (x) -> L` where `L` is the Lipschitz constant of `grad`
+Typically `curv = (x) -> L` where `L` is the Lipschitz constant of `grad`.
 """
 function ogm_ls(
     grad::Function,

src/algorithm/general/pogm_restart.jl

@@ -148,7 +148,7 @@ function pogm_restart(
 	ynew = []
 
 # iterations
-for iter=1:niter
+for iter in 1:niter
 
 	# proximal gradient (PGM) update
 	if mom === :pgm && mu != 0

src/algorithm/general/poweriter.jl

@@ -38,7 +38,7 @@ function poweriter(
 
 	x = copy(x0)
 	ratio_old = Inf
-	for iter=1:niter
+	for iter in 1:niter
 		Ax = A * x
 		ratio = norm(Ax) / norm(x)
 		if abs(ratio - ratio_old) / ratio < tol