Merge pull request #4 from slimgroup/f32

prevent f32 overflow

Author: PetersBas

.github/workflows/ci-test.yml

@@ -19,7 +19,7 @@ jobs:
       fail-fast: false
 
       matrix:
-        version: ['1.5', '1.6', '1.7', 'nightly']
+        version: ['1.6', '1.7']
         os: [ubuntu-latest]
 
         include:

src/PARSDMM.jl

@@ -247,7 +247,7 @@ for i=1:maxit #main loop
 end #end Q-update timer
 
   if i==maxit
-    println("PARSDMM reached maxit")
+    constr_log("PARSDMM reached maxit")
     (TD_OP,AtA,log_PARSDMM) = output_check_PARSDMM(x,TD_OP,AtA,log_PARSDMM,i,counter)
   end
 

src/PARSDMM_initialize.jl

@@ -99,14 +99,14 @@ function PARSDMM_initialize(
     end
   end
   if maximum(feasibility_initial)<options.feas_tol #accept input as feasible and return
-      println("input to PARSDMM is feasible, returning")
+      constr_log("input to PARSDMM is feasible, returning")
       stop = true
   end
 
   # if one of the sets is non-convex, use different lambda and rho update frequency, don't update gamma and set a different fixed gamma
   for ii=1:pp
       if set_Prop.ncvx[ii] == true
-          println("non-convex set(s) involved, using special settings")
+          constr_log("non-convex set(s) involved, using special settings")
           rho_update_frequency  = 3;
           adjust_gamma          = false
           gamma_ini             = TF(0.75)

src/SetIntersectionProjection.jl

@@ -21,6 +21,9 @@ using TimerOutputs
 
 export log_type_PARSDMM, set_properties, PARSDMM_options, set_definitions
 
+const _verbose = false
+constr_log(msg...) = _verbose ? nothing : println(msg...)
+
 #main scripts
 include("PARSDMM.jl")
 include("PARSDMM_multi_level.jl")

src/cg.jl

@@ -62,8 +62,8 @@ function cg(A::Function,b::Vector{TF}; tol::Real=1e-2,maxIter::Integer=100,M::Fu
 
 
 	if out==2
-		println("=== cg ===")
-		println(@sprintf("%4s\t%7s","iter","relres"))
+		constr_log("=== cg ===")
+		constr_log(@sprintf("%4s\t%7s","iter","relres"))
 	end
 
 	resvec = zeros(TF,maxIter)
@@ -99,7 +99,7 @@ function cg(A::Function,b::Vector{TF}; tol::Real=1e-2,maxIter::Integer=100,M::Fu
 		#resvec[iter] = BLAS.nrm2(n, r, 1) / nr0#
 		resvec[iter]  = norm(r)/nr0
 		if out==2
-			println(iter,resvec[iter])
+			constr_log(iter,resvec[iter])
 		end
 		if resvec[iter] <= tol
 			flag = 0; break
@@ -116,12 +116,12 @@ function cg(A::Function,b::Vector{TF}; tol::Real=1e-2,maxIter::Integer=100,M::Fu
 
 	if out>=0
 		if flag==-1
-			println("cg iterated maxIter (=%d) times but reached only residual norm %1.2e instead of tol=%1.2e.",
+			constr_log("cg iterated maxIter (=%d) times but reached only residual norm %1.2e instead of tol=%1.2e.",
 																								maxIter,resvec[lastIter],tol)
 		elseif flag==-2
-			println("Matrix A in cg has to be positive definite.")
+			constr_log("Matrix A in cg has to be positive definite.")
 		elseif flag==0 && out>=1
-			println("cg achieved desired tolerance at iteration %d. Residual norm is %1.2e.",lastIter,resvec[lastIter])
+			constr_log("cg achieved desired tolerance at iteration %d. Residual norm is %1.2e.",lastIter,resvec[lastIter])
 		end
 	end
     return x,flag,resvec[lastIter],lastIter,resvec[1:lastIter]
@@ -194,12 +194,12 @@ end
 #
 # 	if out>=0
 # 		if flag==-1
-# 			println("cg iterated maxIter (=%d) times but reached only residual norm %1.2e instead of tol=%1.2e.",
+# 			constr_log("cg iterated maxIter (=%d) times but reached only residual norm %1.2e instead of tol=%1.2e.",
 # 																								maxIter,resvec[lastIter],tol)
 # 		elseif flag==-2
-# 			println("Matrix A in cg has to be positive definite.")
+# 			constr_log("Matrix A in cg has to be positive definite.")
 # 		elseif flag==0 && out>=1
-# 			println("cg achieved desired tolerance at iteration %d. Residual norm is %1.2e.",lastIter,resvec[lastIter])
+# 			constr_log("cg achieved desired tolerance at iteration %d. Residual norm is %1.2e.",lastIter,resvec[lastIter])
 # 		end
 # 	end
 

src/default_PARSDMM_options.jl

@@ -3,7 +3,7 @@ export default_PARSDMM_options
 """
 Returns a set of default options for the PARSDMM solver
 """
-function default_PARSDMM_options(options,TF)
+function default_PARSDMM_options(options,TF; verbose=false)
 
   if     TF == Float64
     TI = Int64
@@ -28,5 +28,7 @@ function default_PARSDMM_options(options,TF)
   options.parallel              = false       #comput proximal mappings, multiplier updates, rho and gamma updates in parallel
   options.zero_ini_guess        = true        #zero initial guess for primal, auxilliary, and multipliers
   Minkowski                     = false       #the intersection of sets includes a Minkowski set
+
+  _verbose = verbose
   return options
 end

src/projectors/project_l1_Duchi!.jl

@@ -26,7 +26,7 @@ function project_l1_Duchi!(v::Union{Vector{TF},Vector{Complex{TF}}}, b::TF) wher
   u  = similar(v)
   sv = Vector{TF}(undef, lv)
 
-  #use RadixSort for Float32 (short keywords)
+  # use RadixSort for Float32 (short keywords)
   copyto!(u, v)
   u .= abs.(u)
   u  = convert(Vector{TF},u)
@@ -35,19 +35,14 @@ function project_l1_Duchi!(v::Union{Vector{TF},Vector{Complex{TF}}}, b::TF) wher
   else
     u = sort!(u, rev=true, alg=QuickSort)
   end
-  
-  
-  # if TF==Float32
-  #   u = sort!(abs.(u), rev=true, alg=RadixSort)
-  # else
-  #   u = sort!(abs.(u), rev=true, alg=QuickSort)
-  # end
-
   cumsum!(sv, u)
 
   # Thresholding level
-  temp  = TF(1.0):TF(1.0):TF(lv)
-  rho   = max(1, min(lv, findlast(u .> ((sv.-b) ./ temp ) ) ))::Int
+  rho = 0
+  while u[rho+1] > ((sv[rho+1] - b)/(rho+1)) && (rho+1) < lv
+    rho += 1
+  end
+  rho = max(1, rho)
   theta = max.(TF(0) , (sv[rho] .- b) ./ rho)::TF
 
   # Projection as soft thresholding

src/setup_multi_level_PARSDMM.jl

@@ -87,7 +87,7 @@ for i=2:n_levels
   constraint_level = constraint2coarse(constraint_level,comp_grid_levels[i],coarsening_factor)
 
   #set up constraints on new level
-  println(TF)
+  constr_log(TF)
   (P_sub_l,TD_OP_l,set_Prop_l)  = setup_constraints(constraint_level,comp_grid_levels[i],TF)
   (TD_OP_l,AtA_l,dummy1,dummy2) = PARSDMM_precompute_distribute(TD_OP_l,set_Prop_l,comp_grid_levels[i],options)
 

src/stop_PARSDMM.jl

@@ -21,19 +21,19 @@ function stop_PARSDMM(
 
     #stop if objective value does not change and x is sufficiently feasible for all sets
     if i>6 && maximum(log_PARSDMM.set_feasibility[counter-1,:])<feas_tol && maximum(abs.( (log_PARSDMM.obj[i-5:i]-log_PARSDMM.obj[i-1-5:i-1])./log_PARSDMM.obj[i-1-5:i-1] )) < obj_tol
-        println("stationary objective and reached feasibility, exiting PARSDMM (iteration ",i,")")
+        constr_log("stationary objective and reached feasibility, exiting PARSDMM (iteration ",i,")")
         stop=true;
     end
 
     #stop if x doesn't change significantly anyjore
     if i>5 && maximum(log_PARSDMM.evol_x[i-5:i])<evol_rel_tol
-      println("relative evolution to small, exiting PARSDMM (iteration ",i,")")
+      constr_log("relative evolution to small, exiting PARSDMM (iteration ",i,")")
       stop=true;
     end
 
     # fix rho to ensure regular ADMM convergence if primal residual does not decrease over a 20 iteration window
     if i>20 && adjust_rho==true && log_PARSDMM.r_pri_total[i]>maximum(log_PARSDMM.r_pri_total[(i-1):-1:max((i-50),1)])
-      println("no primal residual reduction, fixing PARSDMM rho & gamma (iteration ",i,")")
+      constr_log("no primal residual reduction, fixing PARSDMM rho & gamma (iteration ",i,")")
       adjust_rho = false;
       adjust_feasibility_rho = false;
       adjust_gamma = false;
@@ -47,7 +47,7 @@ function stop_PARSDMM(
 
     #if rho is fixed and still no decrease in primal residual is observed over a window, we give up
     if adjust_rho==false && i>(ind_ref+25) && log_PARSDMM.r_pri_total[i]>maximum(log_PARSDMM.r_pri_total[(i-1):-1:max(ind_ref,max((i-50),1))])
-      println("no primal residual reduction, exiting PARSDMM (iteration ",i,")")
+      constr_log("no primal residual reduction, exiting PARSDMM (iteration ",i,")")
       stop = true;
     end
     return stop,adjust_rho,adjust_gamma,adjust_feasibility_rho,ind_ref