function [a,w,b]=twolayersnorm(n,m,step,itmax)
  pkg load statistics % only in octave
% try with [a,w,b]=twolayersnorm(100,20,1,10000);
% goal: approach any function by sum_j a_j relu(w_j x + b_j)
% n: nb of relu functions
% m: nb of interpolation points
% a: contains the a_j
% w: contains the w_j
% b: contains the b_j
  x=sparse(m,1);
  y=sparse(m,1);
  for i=1:m
    x(i)=-1+2*i/(m+1);
    y(i)=exp(-10*x(i)^2);
  end
  a=normrnd(0,1,n,1);
  w=normrnd(0,1,n,1);
  b=normrnd(0,1,n,1);
  mata=sparse(m,n);
  matw=sparse(m,n);
  matb=sparse(m,n);
% gradient descent
  for iter=1:itmax
   for i=1:m
    for j=1:n
     mata(i,j)=???;
     matb(i,j)=???;
     matw(i,j)=???;
    end
   end
   mata=mata/n;
   matb=matb/n;
   matw=matw/n;
   res=mata*a-y;
   grada=mata'*res;
   gradw=matw'*res;
   gradb=matb'*res;
   a=a-step*grada;
   w=w-step*gradw;
   b=b-step*gradb;
   if (mod(iter,10)==0)
     fprintf(" iter %d norm(grada) %e norm(gradw) %e norm(gradb) %e \n",iter,norm(grada),norm(gradw),norm(gradb));
     plot(x,y,'+',x,mata*a,'x',-b./w,zeros(n,1),'o');   
     xlim([-5 5]);
     ylim([-0.5 1]);
     refresh();
   endif  % in matlab end
  end
end
function relu=relu(x)
  if (x<0)
    relu=0.;
  else
    relu=x;
  endif % in matlab end
end
function relup=relup(x)
  if (x<0)
   relup=0.;
  else 
   relup=1;
  endif % in matlab end
end


