jonathf/matlab2cpp

Lambda function create

Closed this issue · 11 comments

Creating lambda functions is not supported in grammar nor translation

Example:
func = @(x) x^2

This could be translated as e.g:

std::function<double(double)> func = [] (double x) { return x*x; };

Jeg er usikker på om det er en god idé.

A)
Har vi lyst til binde oss til C++11?

B)
Dette fungerer kun hvis typene til paramterene er fastsatt til double. Hva vhsidet ikke er tilfellet?
Mer rigoriøst er å la func defineres som en full fuknsjon med alt det innebærer. Da kan brukeren (eller programvaren) definere typene som de andre funksjonene.
Dette er mye enklere oversettelses-messig.

A) Ja. Detta är mycket grundläggande C++11, alla vettiga kompilatorer har stöd.

B) En std::function är en typ precis som "int" eller "mat". Så i mitt exempel har användaren specificerat "func" till typen "std::function<double(double)>".

A) Ikke mitt tema, men det er noe motstand mot C++11 i forskningsmiljøet
jeg er en del av. Men det blir bare politikk, så jeg har i grunnen ikke
problemer med det.

B) Selve funksjonen ja. Jeg tenkte mer på parameterene

2015-03-17 15:13 GMT+01:00 aronandersson notifications@github.com:

A) Ja. Detta är mycket grundläggande C++11, alla vettiga kompilatorer har
stöd.

B) En std::function är en typ precis som "int" eller "mat". Så i mitt
exempel har användaren specificerat "func" till typen "std::function".

Reply to this email directly or view it on GitHub
#18 (comment).

Jag förstår inte riktigt varför någonting skulle bli annorlunda jämfört med en vanlig funktion - men ok. Tänk på att det kan bli ganska bökigt utan C++11 om man vill göra något icketriviellt med funktionen i fråga e.g:

P = rand(3,1)
func = @(x)(P(3)*x.^2 + P(2)*x + P(1));
x = fzero(func,0);

Det normala sättet att göra detta i "vanlig" C++ är att definiera en klass med en virtuell funktion (som fzero tar som argument), och sen låta en annan klass ärva från den klassen och ha P som en medlemsvariabel etc... Med lambda
kan man skriva:

vec P = rand(3,1);
std::function<double(double)> func = [P] (double x) { return P(2)*x*x+P(1)*x+P(0); };
double x = fzero(afun,0);

Jeg ser poenget ditt med c++11, så vi går for det.

Når det gjelder lambda vs. full func handler det litt om hvor komplisert
kode for håndtering av funksjoner og funksjons-scope er. Men jeg skal nok
klare å lage en løsning som ikke krever at jeg gjenoppfinner hjulet som
lager lambdafunksjoner i output.

2015-03-17 16:17 GMT+01:00 aronandersson notifications@github.com:

Jag förstår inte riktigt varför någonting skulle bli annorlunda jämfört
med en vanlig funktion - men ok. Tänk på att det kan bli ganska bökigt utan
C++11 om man vill göra något icketriviellt med funktionen i fråga e.g:

P = rand(3,1)
func = @(x)(P(3)_x.^2 + P(2)_x + P(1));
x = fzero(func,0);

Det normala sättet att göra detta i "vanlig" C++ är att definiera en klass
med en virtuell funktion (som fzero tar som argument), och sen låta en
annan klass ärva från den klassen och ha P som en medlemsvariabel etc...
Med lambda
kan man skriva:

vec P = rand(3,1);
std::function<double(double)> func = [P](double x) { return P(2)_x_x+P(1)*x+P(0); };
double x = fzero(afun,0);

Reply to this email directly or view it on GitHub
#18 (comment).

Skal være støttet som C++11-kode.

Verkar som om detta kraschar igen, e.g.

1758 19 Assign 'AAt = @(z) (A(At(z)));'
1758 19 Var 'AAt'
Traceback (most recent call last):
File ".\mconvert.py", line 26, in
tree = matlab2cpp.main(opt, args)
File "C:\Anaconda\lib\site-packages\matlab2cpp__init__.py", line 45, in main
unassigned = builder.load(filename)
File "C:\Anaconda\lib\site-packages\matlab2cpp\treebuilder.py", line 64, in load
self.create_program(filename)
File "C:\Anaconda\lib\site-packages\matlab2cpp\treebuilder.py", line 238, in create_program
cur, line = self.create_function(program, cur, line)
File "C:\Anaconda\lib\site-packages\matlab2cpp\treebuilder.py", line 433, in create_function
cur, line = self.fill_codeblock(block, cur, line)
File "C:\Anaconda\lib\site-packages\matlab2cpp\treebuilder.py", line 527, in fill_codeblock
cur, line = self.create_if(block, cur, line)
File "C:\Anaconda\lib\site-packages\matlab2cpp\treebuilder.py", line 1002, in create_if
end, line = self.fill_codeblock(block, cur, line)
File "C:\Anaconda\lib\site-packages\matlab2cpp\treebuilder.py", line 527, in fill_codeblock
cur, line = self.create_if(block, cur, line)
File "C:\Anaconda\lib\site-packages\matlab2cpp\treebuilder.py", line 1002, in create_if
end, line = self.fill_codeblock(block, cur, line)
File "C:\Anaconda\lib\site-packages\matlab2cpp\treebuilder.py", line 555, in fill_codeblock
cur, line = self.create_lambda(block, cur, line, eq_loc)
File "C:\Anaconda\lib\site-packages\matlab2cpp\treebuilder.py", line 1626, in create_lambda
assign[0].set_global_type("func_lambda")
AttributeError: 'Var' object has no attribute 'set_global_type'

Fixed

Ny krasch:

generate translation
iterating 1125 nodes
Traceback (most recent call last):
File ".\mconvert.py", line 26, in
tree = matlab2cpp.main(opt, args)
File "C:\Anaconda\lib\site-packages\matlab2cpp__init__.py", line 79, in main
program.translate_tree(opt)
File "C:\Anaconda\lib\site-packages\matlab2cpp\node.py", line 78, in translate_tree
node.translate_node(opt)
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\func_lambda.py", line 62, in Lambda
node_.declare()
TypeError: 'Get' object is not callable

Fil:

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

largescale = isa(A,'function_handle');

if (nargin < 5), pdtol = 1e-3;  end
if (nargin < 6), pdmaxiter = 50;  end
if (nargin < 7), cgtol = 1e-8;  end
if (nargin < 8), cgmaxiter = 200;  end

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 (largescale)
  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
else
  if (norm(A*x0-b)/norm(b) > cgtol)
    disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
    opts.POSDEF = true; opts.SYM = true;
    [w, hcond] = linsolve(A*A', b, opts);
    if (hcond < 1e-14)
      disp('A*At is ill-conditioned: cannot find starting point');
      xp = x0;
      return;
    end
    x0 = A'*w;
  end  
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;
if (largescale)
  v = -A(lamu1-lamu2);
  Atv = At(v);
  rpri = A(x) - b;
else
  v = -A*(lamu1-lamu2);
  Atv = A'*v;
  rpri = A*x - b;
end

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;

  if (largescale)
    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);
  else
    w1p = -(w3 - A*(w1./sigx - w2.*sig2./(sigx.*sig1)));
    H11p = A*(sparse(diag(1./sigx))*A');
    opts.POSDEF = true; opts.SYM = true;
    [dv,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;
      return
    end
    dx = (w1 - w2.*sig2./sig1 - A'*dv)./sigx;
    Adx = A*dx;
    Atdv = A'*dv;
  end

  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)));
  if (largescale)
    disp(sprintf('                  CG Res = %8.3e, CG Iter = %d', cgres, cgiter));
  else
    disp(sprintf('                  H11p condition number = %8.3e', hcond));
  end

end

Resolved