jonathf/matlab2cpp

Crash when automatically linking multiple files

Closed this issue · 1 comments

Jag får en krash när jag försöker konvertera "l1eq_pd":

Traceback (most recent call last):
  File ".\mconvert.py", line 58, in <module>
    matlab2cpp.main(args)
  File "C:\Anaconda\lib\site-packages\matlab2cpp\__init__.py", line 96, in main
    builder.configure(suggest=2*args.suggest)
  File "C:\Anaconda\lib\site-packages\matlab2cpp\treebuilder.py", line 151, in configure
    node.name+".m")][1][0]

function xp = l1eq_pd(x0, A, At, b, pdtol, pdmaxiter, cgtol, cgmaxiter)

N = length(x0);

alpha = 0.01;
beta = 0.5;
mu = 10;

gradf0 = [zeros(N,1); ones(N,1)];

% starting point --- make sure that it is feasible
if (norm(A(x0)-b)/norm(b) > cgtol)
    disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
    AAt = @(z) A*At*z;
    [w, cgres, cgiter] = cgsolve(AAt, b, cgtol, cgmaxiter, 0);
    if (cgres > 1/2)
        disp('A*At is ill-conditioned: cannot find starting point');
        xp = x0;
        return;
    end
    x0 = At(w);
end

x = x0;
u = (0.95)*abs(x0) + (0.10)*max(abs(x0));

% set up for the first iteration
fu1 = x - u;
fu2 = -x - u;
lamu1 = -1./fu1;
lamu2 = -1./fu2;

v = -A(lamu1-lamu2);
Atv = At(v);
rpri = A(x) - b;

sdg = -(fu1'*lamu1 + fu2'*lamu2);
tau = mu*2*N/sdg;

rcent = [-lamu1.*fu1; -lamu2.*fu2] - (1/tau);
rdual = gradf0 + [lamu1-lamu2; -lamu1-lamu2] + [Atv; zeros(N,1)];
resnorm = norm([rdual; rcent; rpri]);

pditer = 0;
done = (sdg < pdtol) | (pditer >= pdmaxiter);
while (~done)

    pditer = pditer + 1;

    w1 = -1/tau*(-1./fu1 + 1./fu2) - Atv;
    w2 = -1 - 1/tau*(1./fu1 + 1./fu2);
    w3 = -rpri;

    sig1 = -lamu1./fu1 - lamu2./fu2;
    sig2 = lamu1./fu1 - lamu2./fu2;
    sigx = sig1 - sig2.^2./sig1;

    w1p = w3 - A(w1./sigx - w2.*sig2./(sigx.*sig1));
    h11pfun = @(z) -A(1./sigx.*At(z));
    [dv, cgres, cgiter] = cgsolve(h11pfun, w1p, cgtol, cgmaxiter, 0);
    if (cgres > 1/2)
        disp('Cannot solve system.  Returning previous iterate.  (See Section 4 of notes for more information.)');
        xp = x;
        return
    end
    dx = (w1 - w2.*sig2./sig1 - At(dv))./sigx;
    Adx = A(dx);
    Atdv = At(dv);

    du = (w2 - sig2.*dx)./sig1;

    dlamu1 = (lamu1./fu1).*(-dx+du) - lamu1 - (1/tau)*1./fu1;
    dlamu2 = (lamu2./fu2).*(dx+du) - lamu2 - 1/tau*1./fu2;

    % make sure that the step is feasible: keeps lamu1,lamu2 > 0, fu1,fu2 < 0
    indp = find(dlamu1 < 0);  indn = find(dlamu2 < 0);
    s = min([1; -lamu1(indp)./dlamu1(indp); -lamu2(indn)./dlamu2(indn)]);
    indp = find((dx-du) > 0);  indn = find((-dx-du) > 0);
    s = (0.99)*min([s; -fu1(indp)./(dx(indp)-du(indp)); -fu2(indn)./(-dx(indn)-du(indn))]);

    % backtracking line search
    suffdec = 0;
    backiter = 0;
    while (~suffdec)
        xp = x + s*dx;  up = u + s*du;
        vp = v + s*dv;  Atvp = Atv + s*Atdv;
        lamu1p = lamu1 + s*dlamu1;  lamu2p = lamu2 + s*dlamu2;
        fu1p = xp - up;  fu2p = -xp - up;
        rdp = gradf0 + [lamu1p-lamu2p; -lamu1p-lamu2p] + [Atvp; zeros(N,1)];
        rcp = [-lamu1p.*fu1p; -lamu2p.*fu2p] - (1/tau);
        rpp = rpri + s*Adx;
        suffdec = (norm([rdp; rcp; rpp]) <= (1-alpha*s)*resnorm);
        s = beta*s;
        backiter = backiter + 1;
        if (backiter > 32)
            disp('Stuck backtracking, returning last iterate.  (See Section 4 of notes for more information.)')
            xp = x;
            return
        end
    end


    % next iteration
    x = xp;  u = up;
    v = vp;  Atv = Atvp;
    lamu1 = lamu1p;  lamu2 = lamu2p;
    fu1 = fu1p;  fu2 = fu2p;

    % surrogate duality gap
    sdg = -(fu1'*lamu1 + fu2'*lamu2);
    tau = mu*2*N/sdg;
    rpri = rpp;
    rcent = [-lamu1.*fu1; -lamu2.*fu2] - (1/tau);
    rdual = gradf0 + [lamu1-lamu2; -lamu1-lamu2] + [Atv; zeros(N,1)];
    resnorm = norm([rdual; rcent; rpri]);

    done = (sdg < pdtol) | (pditer >= pdmaxiter);

    disp(sprintf('Iteration = %d, tau = %8.3e, Primal = %8.3e, PDGap = %8.3e, Dual res = %8.3e, Primal res = %8.3e',...
        pditer, tau, sum(u), sdg, norm(rdual), norm(rpri)));

    disp(sprintf('                  CG Res = %8.3e, CG Iter = %d', cgres, cgiter));

end

function [x, res, iter] = cgsolve(A, b, tol, maxiter, verbose)

if (nargin < 5), verbose = 1; end

implicit = isa(A,'function_handle');

x = zeros(length(b),1);
r = b;
d = r;
delta = r'*r;
delta0 = b'*b;
numiter = 0;
bestx = x;
bestres = sqrt(delta/delta0); 
while ((numiter < maxiter) && (delta > tol^2*delta0))

  % q = A*d
  if (implicit), q = A(d);  else  q = A*d;  end

  alpha = delta/(d'*q);
  x = x + alpha*d;

  if (mod(numiter+1,50) == 0)
    % r = b - Aux*x
    if (implicit), r = b - A(x);  else  r = b - A*x;  end
  else
    r = r - alpha*q;
  end

  deltaold = delta;
  delta = r'*r;
  beta = delta/deltaold;
  d = r + beta*d;
  numiter = numiter + 1;
  if (sqrt(delta/delta0) < bestres)
    bestx = x;
    bestres = sqrt(delta/delta0);
  end    

  if ((verbose) && (mod(numiter,verbose)==0))
    disp(sprintf('cg: Iter = %d, Best residual = %8.3e, Current residual = %8.3e', ...
      numiter, bestres, sqrt(delta/delta0)));
  end

end

if (verbose)
  disp(sprintf('cg: Iterations = %d, best residual = %14.8e', numiter, bestres));
end
x = bestx;
res = bestres;
iter = numiter;

Prøv nå