function x = shrinkFuncComp(y,alpha,qNorm, method)
f = @(x) x+alpha.*qNorm.*sign(x).*abs(x).^(qNorm-1);
df = @(x) 1+alpha.*qNorm.*(qNorm-1).*abs(x).^(qNorm-2);
maxN = 10000;
tol = 1E-7;
Y = y(:);
tic;
if strcmp(method, 'newton')
x = Y;
diff = f(x)-Y;
dx = 1;
n = 0;
while n<=maxN && max(abs(dx))>tol
n = n+1;
tmp = x;
x(x>0) = max(tol*x(x>0), x(x>0) - diff(x>0)./df(x(x>0)));
x(x<0) = min(tol*x(x<0), x(x<0) - diff(x<0)./df(x(x<0)));
diff = f(x)-Y;
dx = x - tmp;
end
else
x = zeros(size(Y));
lowerBounds = min(x, Y);
upperBounds = max(x, Y);
x = (lowerBounds+upperBounds)/2;
diff = f(x)-Y;
dx = 1;
n = 0;
while n<=maxN && max(abs(dx))>tol
n = n+1;
tmp = x;
upperBounds(diff>0) = x(diff>0);
lowerBounds(diff<0) = x(diff<0);
x = (lowerBounds+upperBounds)/2;
diff = f(x)-Y;
dx = x - tmp;
end
end
t = toc;
if n>maxN
fprintf('WARNING: shrinkFuncComp.m : maximum number of iterations for inversion reached!\n');
fprintf(' n = %i, ||diff|| = %4.2e, ||dx|| = %4.2e, t = %4.4fs\n', n, max(abs(diff)), max(abs(dx)),t);
end
x = reshape(x,size(y));
end