function [im,opts,flog,plog,sclog,rlog] = L2BTV_superres(W,Y,avim,alp,nu,P,irold,opts)
%
% function [im,opts,flog,plog,sclog] = superres_huber(W,Y,La,Lb,avim,alp,nu,opts);
%
% W: W matrix (N*KM)
% Y: Y vector (KM*1)
% La: multiplicative photometric parameter.
% Lb: additive photometric parameter.
% avim: average image (biv*bih)
% alp: Huber alpha parameter.
% nu: (nua) -- ratio of Huber prior strength to image noise precision.
% opts: optional options vector to pass to scg.
%

GRADTEST = 0; % Set this to 1 in order to run local_gradcheck below.

if (nargin<6)
    fprintf('Using default options vector.\n');
    opts = zeros(1,18);
    opts(1) = 1;
    opts(2) = 1e-10;
    opts(3) = 1e-10;
    opts(14) = 40;
end
biv = size(avim,1);
bih = size(avim,2);

KM = size(Y,1);
%W = mex_amub(W,sparse(1:KM,1:KM,La));
W = W*sparse(1:KM,1:KM,ones(KM,1));
% The line above would be "W = W*sparse(1:KM,1:KM,La);", but matlab doesn't
% handle memory very well in this case, so it generates a completely
% unnecessary "out of memory" error. I got mex_amub from the file exchange,
% though now "ssmult" seems to be preferred.

x = avim(:)';

if (GRADTEST)
    [D,G] = local_gradcheck(x,W,Y,alp,nu,biv,bih);
    im = D; opts = G; % overwrite usual outputs when testing grad.
    return;
end


if (opts(1)>=0), fprintf('hmem_ml iterating (%i iters):\n',opts(14)); end;
%if (nargout <4)
    [X,opts,flog,plog,sclog,rlog] = hibscg(@local_func,x,opts,@local_grad,W,Y,alp,nu,biv,bih,P);
% else
%     [X,opts,flog,plog,sclog] = scg(@local_func,x,opts,@local_grad,W,Y,alp,nu,biv,bih);
% end
im = reshape(X,biv,bih);
if (opts(1)>=0), fprintf('Done!\n'); end;
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [e rnew stopsign]= local_func(x,W,Y,alp,nu,biv,bih,P,irold,stopsign)
ee=Y'-x*W;
l=length(ee);
e1=squre_sum(ee,l);
ml=median(ee);
sigml=mean(abs(ee-ml));
sigmg=std(ee,1);
rnew=sigml/sigmg;

deltr=rnew-irold;
if abs(deltr)<.00005;
    stopsign=stopsign+1;
else
    stopsign=0;
end

% l=length(ee);
% if rl>gam
%     %e = sum(e.^2) + elle_eval_huber(x,alp,nu,biv,bih);
%     e1 = squre_sum(ee,l) ;
% else
%     e1 = sum(abs(ee)) ;
% end

%e1 = sum(abs(Y'-x*W));
x1=reshape(x,biv,bih);
G=zeros(biv,bih);
% Create an inflated version of Xn so shifting operation is simpler
Xpad = padarray(x1, [P P], 'symmetric');

% Compute a grid of l=-P:P and m=0:P such that l+m>=0
for l=-P:P
  for m=-P:P

    % Shift HR by l and m
    Xshift = Xpad(1+P-l:end-P-l, 1+P-m:end-P-m);
    % G = G + nu.*alp.^(abs(l)+abs(m)).*(.5-abs( x1-Xshift));
    G = G + nu.*alp.^(abs(l)+abs(m)).*abs( x1-Xshift);

  end
end
e2=sum(abs(G(:)));
e = e1+e2;
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function g = local_grad(x,W,Y,alp,nu,biv,bih,P,gam)
r = Y'-x*W;
%g = W*sign(r)';

%  ml=median(r);
% sigml=mean(abs(r-ml));
% sigmg=std(r,1);
% rl=sigml/sigmg;
%  if rl>gam
%     g = W*r';
%     %g=metric_product(W,r);
% else
%     g = W*sign(r'); 
%     %g=metric_product(W,sign(r'));
%  end
 
g=W*r';
g1 = -2.*g'; 

x1=reshape(x,biv,bih);
G=zeros(biv,bih);
% Create an inflated version of Xn so shifting operation is simpler
Xpad = padarray(x1, [P P], 'symmetric');
% Compute a grid of l=-P:P and m=0:P such that l+m>=0
for l=-P:P
  for m=-P:P

    % Shift HR by l and m
    Xshift = Xpad(1+P-l:end-P-l, 1+P-m:end-P-m);

    % Subtract from HR image and compute sign
    Xsign = sign(x1-Xshift);
    
    % Shift Xsign back by -l and -m
    Xsignpad = padarray(Xsign, [P P], 0);

    Xshift = Xsignpad(1+P+l:end-P+l, 1+P+m:end-P+m);

    G = G + nu.*alp.^(abs(l)+abs(m)).*(Xsign-Xshift);

  end
end
g2=G(:)';
g=g1+g2;
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [D,G] = local_gradcheck(x,varargin)
% This function checks the gradient numerically, for a random subsection of
% the image pixels (since it takes rather too long to do all of them as in
% scg's automatic gradient test option, selected by setting opts(9)=1).

ep = 1e-4; % numerical grdient step size.
numsamp = 30; % approximate number of pixels I'd like to check.

x = x + reshape(gsamp(0,0.15^2,numel(x)),size(x)); % Add in a bit of randome noise.
G = local_grad(x,varargin{:});
L = numel(x);
if (L~=numel(G)), error('lengths of parameter vector and gradient vector are inconsistent.'); end;
step = zeros(1,L);
D = step;
t = numsamp/L;
for i = 1:L
    if (mod(i,10000)<1), fprintf('Element %i of %i.\n',i,L); end;
    if(rand<t)
        step(i) = ep;
        l1 = local_func(x+step,varargin{:});
        l2 = local_func(x-step,varargin{:});
        D(i) = (0.5/ep)*(l1-l2);
        step(i) = 0;
    else
        D(i) = NaN;
    end
end
fprintf('Done all elements!\n');
m = logical(~isnan(D));
G = G(m);
D = D(m);
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%