function principal_avalia_hier2_adap(parametros,plota)
%Encontra uma aproximacao de uma funcao dada com base radial ou tensorial
%hierarquica adaptativa
%Essa versao elimina elementos grandes trocando por elementos menores. E
%resolve o sistema com a base multinivel adaptativa (todos os niveis)

assert(parametros.funcao > 0);
%determina o passo dos niveis
passov = passo_nivel(parametros);

%define a malha completa no nivel maximo
fprintf(2,'definindo a malha regular completa do nivel maximo\n');
assert(parametros.dim_dominio == 2);
%ajusta o tamanho do dominio conforme o tipo de arvore
area = prod(parametros.tam_dominio);
if parametros.tipoarvore == 1
    LD = [sqrt(sqrt(2)), sqrt(sqrt(0.5))]*sqrt(area);
else
    LD = [1, 1]*sqrt(area);
end
fprintf(2,'tamanho do dominio ajustado para %.3f x %.3f\n',LD(1),LD(2));
parametros.tam_dominio = LD;
CD = 0.5*LD;

assert(parametros.n_pontos == 0)
NC_fin = tamanho_malha_hierarquica(parametros.max_nivel);
nc_fin = prod(NC_fin);
fprintf(2,'numero de celulas no nivel maximo %d x %d = %d\n',...
    NC_fin(1),NC_fin(2),nc_fin);

profile on
%determina pontos por celula da malha maxima
KP_fin = calcula_pontos_por_celula(parametros.pts_por_celula,LD,NC_fin);
fprintf(2,'numero de pontos por celula no nivel maximo %d x %d = %d\n',...
    KP_fin(1),KP_fin(2),KP_fin(1)*KP_fin(2));
NP = KP_fin .* NC_fin;
p_comp = constroi_grade_de_amostras(CD,LD,NP);

%tamanho da grade de centros varia com o nivel
assert(isnan(parametros.tam_grade_centros(1)))
assert(isnan(parametros.tam_grade_centros(2)))

%define os pontos de amostragem e a malha do nivel inicial
NC_ini = tamanho_malha_hierarquica(parametros.min_nivel);
nc_ini = prod(NC_ini);
p_ini = p_comp;
np_ini = size(p_ini,1);
np_comp = np_ini; %denominador do produto escalar
malha_ini = constroi_malha_completa(CD,LD,NC_ini,p_ini);
fprintf(2,'numero de celulas no nivel inicial %d x %d = %d\n',...
    NC_ini(1),NC_ini(2),nc_ini);

assert(max(mod(NC_fin,NC_ini)) == 0)
KP_ini = KP_fin .* NC_fin ./ NC_ini;
fprintf(2,'numero de pontos por celula no nivel inicial %d x %d = %d\n',...
    KP_ini(1),KP_ini(2),KP_ini(1)*KP_ini(2));

assert(parametros.raios == parametros.raio)
raio_rel = parametros.raio;

%%define as pastas para os plotes
sigla_metodo = ['MQ';'SS'; 'SO'; 'WQ'];
pasta_raiz = cria_pasta(sprintf(...
    '~/Desktop/dout2012/Texto/Comum/E4apagar-AH2-F%02d',parametros.funcao));
dir1 = sprintf('arv%d-tb%02d-ten%d-dob%d-gs%d-%s',parametros.tipoarvore,...
    parametros.tipobase,parametros.tensorial,1,...
    0,sigla_metodo(parametros.metodo,:));
dir2 = sprintf('NX%05d-NY%05d-r%08.5f',NP(1),NP(2),raio_rel);
pasta_global = cria_pasta(sprintf('%s/%s/%s',...
    pasta_raiz,dir1,dir2));

%constroi a matriz f da funcao a aproximar
%f(i,j) valor da funcao numero j ( = 1) no ponto de amostragem i
f = discretiza_funcao(parametros,p_ini,CD,LD);
assert(size(f,2) == 1);
nf = 1;
if plota.funcao
    figure(9); 
    for j = 1:nf
        plota_funcao(p_ini,f(:,j),true,CD,LD,NP,NaN,NaN);
        %define_zlim2(alpha,f,s);
        title(sprintf('funcao'),'FontSize',20);
        hold off
        if plota.salva_eps
            saveas(gcf,sprintf('%s/f%04d.eps',pasta_raiz,j));
        end
        %         pause(1)
    end
end

arq_global = fopen(sprintf('%s/parametros.tex',pasta_global),'wt');
fprintf(arq_global,'\\newcommand{\\PARMtamdomX}{%d}\n',parametros.tam_dominio(1));
fprintf(arq_global,'\\newcommand{\\PARMtamdomY}{%d}\n',parametros.tam_dominio(2));
fprintf(arq_global,'\\newcommand{\\PARMtensorial}{%d}\n',parametros.tensorial);
%fprintf(arq_global,'\\newcommand{\\PARMdobrado}{%d}\n',parametros.dobrado);
fprintf(arq_global,'\\newcommand{\\PARMamostrX}{%d}\n',NP(1));
fprintf(arq_global,'\\newcommand{\\PARMamostrY}{%d}\n',NP(2));
fprintf(arq_global,'\\newcommand{\\PARMamostras}{%d}\n',np_ini);
fprintf(arq_global,'\\newcommand{\\PARMraio}{%s}\n',...
    converte_raio_tex(raio_rel));

tol = parametros.tol;
NC_ant = [NaN NaN];
NC_v = NC_ini;
KP_v = KP_ini;
p_h = p_ini; %sao os pontos de amonstragem
np_h = np_ini; %numero de pontos de amonstragem
e_h = f;
for iv = parametros.min_nivel:passov:parametros.max_nivel
    fprintf(2,'=== nivel %d ======================================\n',iv);
    %malha_h(iv) eh a malha com as celulas relevantes do nivel iv
    %base_h eh a base hierarquica (multinivel) atual
    %nb_h eh numero total de elementos na base hierarquica
    %A_h  sao osvalores de todos os elementos nos pontos de p
    %M_h eh a matriz do sistema MQ
    %inib_h(iv) eh o comecco do nivel iv em A_h, alpha_h, base_h
    %tamb_h(iv) eh o numero de elementos da base no nivel iv
    
    plota.janela_base = 20;
    pasta_nivel = cria_pasta(sprintf('%s/iv%02d',...
        pasta_global,iv));
    
    %ajusta as malhas usando ate o nivel iv se necessario
    if iv == parametros.min_nivel
        malha_h(iv) = malha_ini; %malha inicial
        NC_iv = malha_h(iv).NC;
        raio_abs_iv = (LD./NC_iv)*raio_rel;
        %base_h = constroi_base(malha_ini,[],raio_abs_iv,...
        %   parametros.tipobase);
        KP_iv = KP_ini;
        %[A_h,M_h] = obtem_base_discreta(parametros,p_h,malha_ini,base_h,...
        %   parametros.tensorial,parametros.dobrado,KP_iv,NP,np_comp,1:np_h);
        inib_h(iv) = 1;
        %   tamb_h(iv) = size(A_h,2);
        tamb_h(iv) = 0;
        A_h = sparse(np_comp,0);
        base_h.elementos = struct(zeros(1,0));
        base_h.uniforme = false;
    else
        KF = numero_de_filhos(parametros.tipoarvore,parametros.dim_dominio,...
            iv-passov);
        fprintf(2,'numero de filhos por celula da malha anterior = [%d %d] \n',...
            KF(1),KF(2));
        KP_iv = KP_ant ./ KF;
        NC_iv = NC_ant .* KF;
        raio_abs_iv = (LD./NC_iv)*raio_rel;
        [malha_h,inib_h,tamb_h] = acrescenta_nivel_vazio_hier2(malha_h,...
            base_h,inib_h,tamb_h,iv,CD,LD,NC_iv);
    end
    fprintf(2,'raio absoluto do nivel %d = [%8.5f %8.5f]\n',...
        iv,raio_abs_iv(1),raio_abs_iv(2));
    %acrescenta elementos no nivel iv ateh que todas as celulas criticas
    %deste nivel tenham elementos centrados nelas.
    estagio = sprintf('nivel %02d (preb)',iv);
    iter = 0;
    while true
        fprintf(2,'+++ inicio da iteracao %d +++++++\n',iter);
        nb_ant = size(base_h.elementos,2);
        [malha_h,base_h,A_h,inib_h,tamb_h]  = ...
            acrescenta_elementos_criticos_hier2(parametros,malha_h,p_h,...
            e_h,base_h,A_h,inib_h,tamb_h,iv,NC_iv,KP_iv,raio_rel,...
            plota,pasta_nivel,estagio);
        [malha_h,base_h,A_h,inib_h,tamb_h] = ...
            elimina_redundantes_hier2(parametros,malha_h,base_h,...
            iv,NC_iv,CD,LD,NP,A_h,inib_h,tamb_h,...
            p_h,plota,pasta_nivel,estagio);

        nb_h = size(base_h.elementos,2);
        assert(issparse(A_h))
        for jv = parametros.min_nivel:passov:iv
            NC_jv = malha_h(jv).NC;
            fprintf(2,'celulas na malha ajustada do nivel %2d = %6d',...
                jv,nnz(malha_h(jv).ixc));
            
            fprintf(2,' na malha completa = [%6d %6d] = %6d\n',...
                NC_jv(1),NC_jv(2),prod(NC_jv));
        end
        M_h = (A_h'*A_h)/np_comp;
        assert(issparse(M_h))
        M_h = sparse(M_h);
        if nb_h == nb_ant
            break
        end
        
        figure(10 + plota.janela_base); clf;
        plota_centros_hier2(base_h,p_h,A_h,malha_h,NP,...
            parametros,parametros.min_nivel,passov,iv,inib_h,tamb_h,...
            pasta_nivel,plota,estagio)
        
        assert((parametros.metodo == 1) || (parametros.metodo == 3)...
             || (parametros.metodo == 4)) 
        
        %o problema eh aproximar a funcao f  pela base A_h nos pontos p_h
        grava_tex = true;
        [b_h,alpha_h,s_h,e_h,map_h,condM_h] = aproximacao(parametros,NP,...
            p_h,np_comp,CD,LD,base_h,A_h,f,M_h,pasta_nivel,plota,...
            estagio,grava_tex);
        %alpha_h coeficientes desses elementos na aproximacao
        %s_h sao os valores da aproximacao nos pontos de p
        %e_h eh oresiduo da aprox s_h nos pontos de p
        
        erro_h_L2 = norma_erro_L2(e_h,np_comp);
        erro_h_Linf = norma_erro_Linf(e_h);
        fprintf(2,'erro residual com a base  do nivel = %9.6f (L2) %9.6f (Linf)\n', ...
            erro_h_L2, erro_h_Linf);
        fprintf(2,'+++ fim da iteracao %d +++++++\n',iter);
        iter = iter + 1;
    end  
    
    arq = fopen(sprintf('%s/parametros.tex',pasta_nivel),'wt');
    fprintf(arq,'\\newcommand{\\PARMraioX}{%s}\n',...
        converte_raio_tex(raio_abs_iv(1)));
    fprintf(arq,'\\newcommand{\\PARMraioY}{%s}\n',...
        converte_raio_tex(raio_abs_iv(2)));
    fprintf(arq,'\\newcommand{\\PARMptporcelX}{%d}\n',KP_iv(1));
    fprintf(arq,'\\newcommand{\\PARMptporcelY}{%d}\n',KP_iv(2));
    fprintf(arq,'\\newcommand{\\PARMptporcelula}{%d}\n',KP_iv(1)*KP_iv(2));
    fprintf(arq,'\\newcommand{\\PARMcentrosX}{%d}\n',NC_iv(1));
    fprintf(arq,'\\newcommand{\\PARMcentrosY}{%d}\n',NC_iv(2));
    
    fprintf(2,'tamanho da base atual %d\n',nb_h);
    fprintf(arq,'\\newcommand{\\PARMdimA}{%d}\n',nb_h);
    
    if plota.base
        plota_base_multinivel(parametros,malha_h,NP,p_h,A_h,base_h,...
            pasta_nivel,plota,estagio);
    end
    if parametros.adap_por_nivel
        fprintf(2,'--- aproximacao adaptativa ate nivel %d -------------\n',iv);
        estagio = sprintf('nivel %02d (adap)',iv);
        plota.janela_base = 20;  %reusa as janelas
        pasta_nivel = cria_pasta(sprintf('%s-a',...
            pasta_nivel));
        plota.base_dual = false;
        [base_a,A_a,M_a,inib_a,tamb_a,alpha_a] = ...
            escolhe_sub_base_hier2(parametros,base_h,A_h,M_h,...
            parametros.min_nivel,passov,iv,inib_h,...
            tamb_h,alpha_h,e_h);
        nb_a = size(A_a,2);
        fprintf(2,'tam base adaptativa %d\n',nb_a);
        verifica_base_hier2(base_a,A_a,parametros.min_nivel,passov,...
            iv,inib_a,tamb_a,p_h)
        
        figure(100); clf;
        plota_centros_hier2(base_a,p_h,A_a,malha_h,NP,...
            parametros,parametros.min_nivel,passov,iv,inib_a,tamb_a,...
            pasta_nivel,plota,estagio)
        
        %calcula residuo com coeficientes originais e base reduzida
        s_a = A_a*alpha_a;
        e_a = f - s_a;
        erro_a_L2 = norma_erro_L2(e_a,np_comp);
        erro_a_Linf = norma_erro_Linf(e_a);
        fprintf(2,'erro com base reduzida e alphas originais %9.6f (L2) %9.6f (Linf)\n',...
            erro_a_L2, erro_a_Linf);
        arq_a = fopen(sprintf('%s/paramadapt.tex',pasta_nivel),'wt');
        fprintf(arq_a,'\\newcommand{\\PARMerroaHilb}{%s}\n',...
            converte_erro_tex(erro_a_L2));
        fprintf(arq_a,'\\newcommand{\\PARMerroaLinf}{%s}\n',...
            converte_erro_tex(erro_a_Linf));
        fprintf(arq_a,'\\newcommand{\\PARMtol}{%s}\n',...
            converte_erro_tex(tol));
        fprintf(arq_a,'\\newcommand{\\PARMnba}{%d}\n',nb_a);
        
        assert(nb_a > 0, 'a funcao eh menor que a tolerancia')
        
        if (parametros.hier_adap_recalc)
            estagio = '(recalc)';
            pasta_nivel_recalc = cria_pasta(sprintf('%s/iv%02dr',...
                pasta_global,iv));
            grava_tex = true;
            [b_a,alpha_a,s_a,e_a,map_a,condM_a] = aproximacao(parametros,NP,...
                p_h,np_comp,CD,LD,base_a,A_a,f,M_a,pasta_nivel_recalc,...
                plota,estagio,grava_tex);
            erro_v_L2 = norma_erro_L2(e_h,np_comp);
            erro_v_Linf = norma_erro_Linf(e_h);
            fprintf(2,'erro com base reduzida e alphas recalculados %9.6f (L2) %9.6f (Linf)\n',...
                erro_a_L2, erro_a_Linf);
            fprintf(arq_a,'\\newcommand{\\PARMerroarHilb}{%s}\n',...
                converte_erro_tex(erro_a_L2));
            fprintf(arq_a,'\\newcommand{\\PARMerroarLinf}{%s}\n',...
                converte_erro_tex(erro_a_Linf));
        end
        fclose(arq_a);
        base_h = base_a;
        nb_h = nb_a;
        A_h = A_a;
        M_h = M_a;
        inib_h = inib_a;
        tamb_h = tamb_a;
        alpha_h = alpha_a;
        s_h = s_a;
        e_h = e_a;
        erro_h_L2 = erro_a_L2;
        erro_h_Linf = erro_a_Linf;
        
    end
    fclose(arq);
    
    if erro_h_Linf <= tol
        fprintf(2,'criterio satisfeito no nivel %d\n',iv);
        break
    end
    
    KP_ant = KP_iv;
    NC_ant = NC_iv;
end
fprintf(arq_global,'\\newcommand{\\PARMnbh}{%d}\n',nb_h);
erro_h_L2 = norma_erro_L2(e_h,np_comp);
erro_h_Linf = norma_erro_Linf(e_h);
fprintf(arq_global,'\\newcommand{\\PARMerrohHilb}{%s}\n',...
    converte_erro_tex(erro_h_L2));
fprintf(arq_global,'\\newcommand{\\PARMerrohLinf}{%s}\n',...
    converte_erro_tex(erro_h_Linf));
fclose(arq_global);

plota.janela_base = 40;
estagio = '(global)';
nomef = 'funcao';
pasta = pasta_global;
figure(5+plota.janela_base); 
plota_funcao(p_comp,s_h(:,1),true,CD,LD,NP,NaN,NaN);
title(sprintf('aproximacao da %s %s',...
    nomef,estagio),'FontSize',20);
hold off
if plota.salva_eps
    saveas(gcf,sprintf('%s/s_h.eps',pasta));
end

figure(6+plota.janela_base);
if plota.erro_curva_nivel
    plot_contour_funcao(parametros,p_comp,e_h(:,1),CD,LD,NP);
else
    emax = max([plota.emax, max(abs(e_h(:,1)))]);
    plota_funcao(p_comp,e_h(:,1),true,CD,LD,NP,-emax,emax);
end
title(sprintf('erro da %s %s',...
    nomef,estagio),'FontSize',20);
hold off
if plota.salva_eps
    saveas(gcf,sprintf('%s/e_h.eps',pasta));
end

profile viewer
%deveria tentar achar a aproximacao e/ou plotar com a base inteira