%  Multirresolución mediante la recosntrucción AC con SR, pero utilizando la medida E_2 para la tecnica ENO. 
%  Además se utiliza la idea de aproximar TANTO EN PARES COMO IMPARES y
%  posteriormente SE INTRODUCE CONSISTENCIA.
%    capeliminacapiamincmcenosruideapacoe2(xx,ff,fsin,nl,nr,grau,nivells)
%    Entradas:
%           xx       = nodos de la funcion inicial
%           ff       = f(xx)
%           fsin     = f sin ruido
%           nl       = numero de nodos a la izquierda de f^{k-1}_{i} 
%           nl       = numero de nodos a la derecha de f^{k-1}_{i}         
%           grau     = grado para F (grau-1 para F)
%           nivells  = niveles de multiresolucion 
%    Salidas:
%          Plot de la reconstruccion(-), funcion de entrada con  ruido (.), funcion original sin ruido (--).
%
% Autor: José Jaime Noguera.   +INFO: www.uv.es/anims/noguera

function capeliminacapiamincmcenosruideapacoe2(xx,ff,fsin,nl,nr,grau,nivells)

nivells=nivells+1;
[xn,xm]=size(xx);
n=log(xm)/log(2);
f=zeros(nivells,xm);
x=zeros(nivells,xm);
f(nivells,1:xm)=ff(1:xm);
x(nivells,1:xm)=xx(1:xm);
ntot=nl+nr+1;

for k=nivells:-1:2   
   %bajamos niveles mediante Medias en Celda
   for i=1:2^(n-(nivells-k+1)) 
      f(k-1,i)=(1/2)*(f(k,2*i-1)+f(k,2*i));
      x(k-1,i)=(x(k,2*i)+x(k,2*i-1))/2;  % solo lo uso al final para el plot, ya que trabajo en nodos enteros
   end
end
    
   fr=zeros(nivells,xm);
   fr(1,1:2^(n-(nivells-1)))=f(1,1:2^(n-(nivells-1))); %para asegurarme de que toda la reconstruccion es nueva

for k=2:nivells
   ultim=2^(n-(nivells-k+1));
   %%%%%%%%%%%%ENO ADAPTADO A E2 PARA MEDIAS EN CELDA %%%%%%%%%%%%%%%%
   for i=1:ultim

      for l=1:ntot
          suma=0;
           if(i+l-1<=ultim && i-ntot+l>0)
               for s=1:ntot           
                  suma=(capiamincmccoefisolupunto(ntot-l,l-1,grau,i-ntot+l+s-1-i+1)*fr(k-1,i-ntot+l:i+l-1)'-fr(k-1,i-ntot+l+s-1))^2+suma; 
               end
                   dif(l)=suma;
           else
                  dif(l)=1e100;  %le doy un valor excesivo para tocar menos codigo en el tratamiento de fronteras  
           end
      end
      [minim,pos]=min(dif(1:ntot));
      vnl(k,i)=ntot-pos; %num nodos a la izquierda de f^{k}_{i} para ENO
      vnr(k,i)=nl+nr-vnl(k,i);
      
   end 
   % FIN ENO %%%%%%%%%%%%%%%

   for i=2:ultim-1
      if(vnr(k,i-1)==0 && vnl(k,i+1)==0)

                   Gl=capiaiamccoefsrgenusol(0,fr,nl,nr,grau,i,k);
                   Gc=capiaiamccoefsrgenusol(0.5,fr,nl,nr,grau,i,k);
                   Gr=capiaiamccoefsrgenusol(1,fr,nl,nr,grau,i,k);
                    
                   if(Gl*Gr<0) %Sospechosa
                        if(Gl*Gc<0)  %Singular
                              %debemos extender de la derecha
                              fr(k,2*i)=capiamincsrdretau(fr,nl,nr,grau,i,k);     
                              fr(k,2*i-1)=capiamincsrdretauidea(fr,nl,nr,grau,i,k);
                                 %para consistencia
                                 fra=fr(k,2*i-1);
                                 frb=fr(k,2*i);
                                 fr(k,2*i-1)=fra+fr(k-1,i)-(1/2)*(fra+frb);
                                 fr(k,2*i)=frb+fr(k-1,i)-(1/2)*(fra+frb);                           
                         else
                             %debemos extender de la izquierda
                             fr(k,2*i-1)=capiamincsresquerrau(fr,nl,nr,grau,i,k);
                             fr(k,2*i)=capiamincsresquerrauidea(fr,nl,nr,grau,i,k);   
                                 %para consistencia
                                 fra=fr(k,2*i-1);
                                 frb=fr(k,2*i);
                                 fr(k,2*i-1)=fra+fr(k-1,i)-(1/2)*(fra+frb);
                                 fr(k,2*i)=frb+fr(k-1,i)-(1/2)*(fra+frb);
                         end
                   else % hacemos ENO
                             fr(k,2*i-1)=capiamincmccoefisolu(vnl(k,i),vnr(k,i),grau)*fr(k-1,i-vnl(k,i):i+vnr(k,i))';  
                             fr(k,2*i)=capiamincmccoefisoluidea(vnl(k,i),vnr(k,i),grau)*fr(k-1,i-vnl(k,i):i+vnr(k,i))';
                                  %para consistencia
                                  fra=fr(k,2*i-1);
                                  frb=fr(k,2*i);
                                  fr(k,2*i-1)=fra+fr(k-1,i)-(1/2)*(fra+frb);
                                  fr(k,2*i)=frb+fr(k-1,i)-(1/2)*(fra+frb);
                   end
      else % hacemos ENO
         
         fr(k,2*i-1)=capiamincmccoefisolu(vnl(k,i),vnr(k,i),grau)*fr(k-1,i-vnl(k,i):i+vnr(k,i))';  
         fr(k,2*i)=capiamincmccoefisoluidea(vnl(k,i),vnr(k,i),grau)*fr(k-1,i-vnl(k,i):i+vnr(k,i))';
               %para consistencia
               fra=fr(k,2*i-1);
               frb=fr(k,2*i);
               fr(k,2*i-1)=fra+fr(k-1,i)-(1/2)*(fra+frb);
               fr(k,2*i)=frb+fr(k-1,i)-(1/2)*(fra+frb);
      end   
   end  
   %Fronteras
   i=1; 
   fr(k,2*i-1)=capiamincmccoefisolu(vnl(k,i),vnr(k,i),grau)*fr(k-1,i-vnl(k,i):i+vnr(k,i))';  
   fr(k,2*i)=capiamincmccoefisoluidea(vnl(k,i),vnr(k,i),grau)*fr(k-1,i-vnl(k,i):i+vnr(k,i))';
   i=ultim; 
   fr(k,2*i-1)=capiamincmccoefisolu(vnl(k,i),vnr(k,i),grau)*fr(k-1,i-vnl(k,i):i+vnr(k,i))';  
   fr(k,2*i)=capiamincmccoefisoluidea(vnl(k,i),vnr(k,i),grau)*fr(k-1,i-vnl(k,i):i+vnr(k,i))';    
   %Fin_fronteras  
end   

    disp('Entre fr y fsin')
    error=mse(fsin,fr(nivells,1:xm));
    error
 
    disp('Entre fr y fn')
    error=mse(ff,fr(nivells,1:xm));
    error
    
    disp('Entre fr y fn')
    error=calculosnr(fr(nivells,1:xm),ff);
    error
    
    disp('Entre fin y fr')
    error=calculosnr(fsin,fr(nivells,1:xm));
    error
   
hold on

%axis([-1,1,-0.1,1.5])
plot(xx, ff,'k.','LineWidth',5,'MarkerSize',7)   
plot(xx, fsin,'--k','LineWidth',2) 
plot(xx,fr(nivells,1:xm),'-k','LineWidth',3);

%axis([-1 1 -1.2 1.2 ]);
%set(gca,'Xtick',[-1 -0.5 0 0.5 1], 'FontSize',15)
%set(gca,'Ytick',[-1 -0.5 0 0.5 1], 'FontSize',15)

%axis([0 1 -0.25 1.25 ]);
%set(gca,'Xtick',[0 .5 1], 'FontSize',15)
%set(gca,'Ytick',[0 .5 1], 'FontSize',15)
box
hold off
