function x = shrinkFuncComp(y,alpha,qNorm, method)
% solve f(x)=y for x yields shrinkage function

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;

    % use dx as condition due to high slope at x=0
    while n<=maxN && max(abs(dx))>tol
        n = n+1;
        tmp = x;

        % newton step, but do not let x change its sign
        % (f is strictly increasing + f(0) = 0)
        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 % bisection
    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