Krasch

[aronandersson]

configuring tvqc_newton.m
configure tree
Traceback (most recent call last):
  File ".\mconvert.py", line 35, in <module>
    tree = matlab2cpp.main(args)
  File "C:\Anaconda\lib\site-packages\matlab2cpp\__init__.py", line 71, in main
    builder.configure()
  File "C:\Anaconda\lib\site-packages\matlab2cpp\treebuilder.py", line 171, in configure
    node.translate_node()
  File "C:\Anaconda\lib\site-packages\matlab2cpp\node.py", line 107, in translate_node
    value = value(node)
  File "C:\Anaconda\lib\site-packages\matlab2cpp\targets\matrix.py", line 98, in Matrix
    return str(node.auxiliary())
  File "C:\Anaconda\lib\site-packages\matlab2cpp\node.py", line 209, in auxiliary
    rhs.translate_node()
  File "C:\Anaconda\lib\site-packages\matlab2cpp\node.py", line 107, in translate_node
    value = value(node)
  File "C:\Anaconda\lib\site-packages\matlab2cpp\targets\matrix.py", line 161, in Assign
    assert rhs.cls in ("Matrix", "Cell")
AssertionError

fil:

function [xp, tp, niter] = tvqc_newton(x0, t0, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter) 

%largescale = isa(A,'function_handle'); 

alpha = 0.01;
beta = 0.5;  

N = length(x0);
n = round(sqrt(N));

% create (sparse) differencing matrices for TV
Dv = spdiags([reshape([-ones(n-1,n); zeros(1,n)],N,1) ...
  reshape([zeros(1,n); ones(n-1,n)],N,1)], [0 1], N, N);
Dh = spdiags([reshape([-ones(n,n-1) zeros(n,1)],N,1) ...
  reshape([zeros(n,1) ones(n,n-1)],N,1)], [0 n], N, N);

AtA = A'*A;

% initial point
x = x0;
t = t0;

r = A*x - b;
Dhx = Dh*x;  Dvx = Dv*x;
ft = 1/2*(Dhx.^2 + Dvx.^2 - t.^2);
fe = 1/2*(r'*r - epsilon^2);
f = sum(t) - (1/tau)*(sum(log(-ft)) + log(-fe));

niter = 0;
done = 0;
while (~done)

  Atr = A'*r;
  ntgx = Dh'*((1./ft).*Dhx) + Dv'*((1./ft).*Dvx) + 1/fe*Atr;
  ntgt = -tau - t./ft;
  gradf = -(1/tau)*[ntgx; ntgt];

  sig22 = 1./ft + (t.^2)./(ft.^2);
  sig12 = -t./ft.^2;
  sigb = 1./ft.^2 - (sig12.^2)./sig22;

  w1p = ntgx - Dh'*(Dhx.*(sig12./sig22).*ntgt) - Dv'*(Dvx.*(sig12./sig22).*ntgt);

    H11p =  Dh'*sparse(diag(-1./ft + sigb.*Dhx.^2))*Dh + ...
      Dv'*sparse(diag(-1./ft + sigb.*Dvx.^2))*Dv + ...
      Dh'*sparse(diag(sigb.*Dhx.*Dvx))*Dv + ...
      Dv'*sparse(diag(sigb.*Dhx.*Dvx))*Dh - ...
      (1/fe)*AtA + (1/fe^2)*Atr*Atr';
    opts.POSDEF = true; opts.SYM = true;
    [dx,hcond] = linsolve(H11p, w1p, opts);
    if (hcond < 1e-14)
      disp('Matrix ill-conditioned.  Returning previous iterate.  (See Section 4 of notes for more information.)');
      xp = x;  tp = t;
      return
    end
    Adx = A*dx;

  Dhdx = Dh*dx;  Dvdx = Dv*dx;
  dt = (1./sig22).*(ntgt - sig12.*(Dhx.*Dhdx + Dvx.*Dvdx));

  % minimum step size that stays in the interior
  aqt = Dhdx.^2 + Dvdx.^2 - dt.^2;   
  bqt = 2*(Dhdx.*Dhx + Dvdx.*Dvx - t.*dt);  
  cqt = Dhx.^2 + Dvx.^2 - t.^2;
  tsols = [(-bqt+sqrt(bqt.^2-4*aqt.*cqt))./(2*aqt); ...
    (-bqt-sqrt(bqt.^2-4*aqt.*cqt))./(2*aqt) ];
  indt = find([(bqt.^2 > 4*aqt.*cqt); (bqt.^2 > 4*aqt.*cqt)] & (tsols > 0));
  aqe = Adx'*Adx;   bqe = 2*r'*Adx;   cqe = r'*r - epsilon^2;
  smax = min(1,min([...
    tsols(indt); ...
    (-bqe+sqrt(bqe^2-4*aqe*cqe))/(2*aqe)
    ]));
  s = (0.99)*smax;

  % backtracking line search
  suffdec = 0;
  backiter = 0;
  while (~suffdec)
    xp = x + s*dx;  tp = t + s*dt;
    rp = r + s*Adx;  Dhxp = Dhx + s*Dhdx;  Dvxp = Dvx + s*Dvdx;
    ftp = 1/2*(Dhxp.^2 + Dvxp.^2 - tp.^2);
    fep = 1/2*(rp'*rp - epsilon^2);
    fp = sum(tp) - (1/tau)*(sum(log(-ftp)) + log(-fep));
    flin = f + alpha*s*(gradf'*[dx; dt]);
    suffdec = (fp <= flin);
    s = beta*s;
    backiter = backiter + 1;
    if (backiter > 32)
      disp('Stuck on backtracking line search, returning previous iterate.  (See Section 4 of notes for more information.)');
      xp = x;  tp = t;
      return
    end
  end

  % set up for next iteration
  x = xp; t = tp;
  r = rp;  Dvx = Dvxp;  Dhx = Dhxp; 
  ft = ftp; fe = fep; f = fp;

  lambda2 = -(gradf'*[dx; dt]);
  stepsize = s*norm([dx; dt]);
  niter = niter + 1;
  done = (lambda2/2 < newtontol) | (niter >= newtonmaxiter);

  disp(sprintf('Newton iter = %d, Functional = %8.3f, Newton decrement = %8.3f, Stepsize = %8.3e', ...
    niter, f, lambda2/2, stepsize));
    disp(sprintf('                  H11p condition number = %8.3e', hcond));


end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% H11p auxiliary function
function y = H11p(v, A, At, Dh, Dv, Dhx, Dvx, sigb, ft, fe, atr)

Dhv = Dh*v;
Dvv = Dv*v;

y = Dh'*((-1./ft + sigb.*Dhx.^2).*Dhv + sigb.*Dhx.*Dvx.*Dvv) + ...
  Dv'*((-1./ft + sigb.*Dvx.^2).*Dvv +A sigb.*Dhx.*Dvx.*Dhv) - ...
  1/fe*At(A(v)) + 1/fe^2*(atr'*v)*atr;  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%