Commit 6f31e72d authored by Russell Luke's avatar Russell Luke
Browse files

Tried naive implementation of MagProj...doesn't work!

parent 62221e42
......@@ -47,28 +47,40 @@ class ProxOperator:
def magproj(u, constr):
"""
Projection operator onto a magnitude constraint
Parameters
----------
u : array_like - The function to be projected onto constr (can be complex)
constr : array_like - A nonnegative array that is the magnitude constraint
Returns
-------
array_like - The projection
"""
eps = 1e-23
"""
Projection operator onto a magnitude constraint
modsq_u = conj(u) * u
denom = modsq_u+eps
denom2 = sqrt(denom)
r_eps = (modsq_u/denom2) - constr
dr_eps = (denom+eps)/(denom*denom2)
return (1 - (dr_eps*r_eps)) * u
Parameters
----------
u : array_like - The function to be projected onto constr (can be complex)
constr : array_like - A nonnegative array that is the magnitude constraint
Returns
-------
array_like - The projection
"""
"""
# naive implementation:
mod_u = sqrt(conj(u) * u)
return constr*(u/mod_u)
"""
"""
Inexact, but stable implementation of magnitude projection.
See LukeBurkeLyon, SIREV 2002
"""
eps = 1e-23
modsq_u = conj(u) * u
denom = modsq_u+eps
denom2 = sqrt(denom)
r_eps = (modsq_u/denom2) - constr
dr_eps = (denom+eps)/(denom*denom2)
return (1 - (dr_eps*r_eps)) * u
class P_diag(ProxOperator):
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment