function principal_avalia_hier_adap(parametros,plota)
%Encontra uma aproximacao de uma funcao dada com base radial ou tensorial 
%hierarquica adaptativa

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;
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
assert(parametros.n_pontos == 0)
%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/E4-AHthetameio-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_ini = discretiza_funcao(parametros,p_ini,CD,LD);
assert(size(f_ini,2) == 1);
nf = 1;
if plota.funcao
    figure(9);
    for j = 1:nf
        plota_funcao(p_ini,f_ini(:,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));

NC_ant = [NaN NaN];
NC_v = NC_ini;
KP_v = KP_ini;
nh = 0; %numero de malhas e bases criadas 
nb_h =0; %numero total de elementos nessas bases
p_h = p_ini; %sao os pontos de amonstragem da aproximacao acumulada
np_h = np_ini; %numero de pontos de amonstragem da aproximacao acumulada
A_h = []; %valores de todos os elementos nos pontos de amostragem de p
alpha_h = []; %coeficientes desses elementos na aproximacao acumulada 
s_h = zeros(np_ini,1); %valores da aproximacao acumulada nos pontos de p %aqui
e_h = f_ini; %residuo da aprox acumulada nos pontos de p
inib_h = []; %inib_h(iv) eh o comecco do nivel iv em A_h, alpha_h  %aqui
beta = 1; %reducao de tolerancia por nivel 
theta = 0.5; %tolerancia relativa para a escolha das celulas criticas
for iv = parametros.min_nivel:passov:parametros.max_nivel
    fprintf(2,'=== nivel %d ======================================\n',iv);
    %neste ponto e_h e s_h sao o residuo e aproximacao acumulada ate o
    %nivel anterior. malha_ant eh a malha formada pelas celulas relevantes
    %usandas no nivel anterior. 
    %calcula a tolerancia a usar neste nivel
    tol = parametros.tol*beta^(parametros.max_nivel - iv);
    fprintf(2,'tolerancia = %8.6f \n',tol);
    
    %determina uma malha inicial malha_iv refinando a malha do nivel
    %anterior (malha_ant)
    if iv == parametros.min_nivel
        malha_iv = malha_ini; %malha inicial
        KP_v = KP_ini;
    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));
        malha_iv = refina_malha(malha_ant,KF,p_h);
        assert(max(mod(KP_ant,KF)) == 0)
        KP_v = KP_ant ./ KF;
    end
    NC_v = malha_iv.NC;
    fprintf(2,'numero de celulas na malha inicial do nivel = %d \n',...
        nnz(malha_iv.ixc));
    
    fprintf(2,'numero de celulas na malha completa = [%d %d] \n',...
        NC_v(1),NC_v(2));
    raio_abs_v = (LD./NC_v)*raio_rel;%raio absoluto do elemento em cada eixo
    fprintf(2,'raio absoluto = [%8.5f %8.5f]\n',...
        raio_abs_v(1),raio_abs_v(2));
    
    plota.janela_base = 20;
    pasta_nivel = cria_pasta(sprintf('%s/iv%02d',...
        pasta_global,iv));
    arq = fopen(sprintf('%s/parametros.tex',pasta_nivel),'wt');
    fprintf(arq,'\\newcommand{\\PARMraioX}{%s}\n',...
        converte_raio_tex(raio_abs_v(1)));
    fprintf(arq,'\\newcommand{\\PARMraioY}{%s}\n',...
        converte_raio_tex(raio_abs_v(2)));
    fprintf(arq,'\\newcommand{\\PARMcentrosX}{%d}\n',NC_v(1));
    fprintf(arq,'\\newcommand{\\PARMcentrosY}{%d}\n',NC_v(2));
    fprintf(arq,'\\newcommand{\\PARMptporcelX}{%d}\n',KP_v(1)); 
    fprintf(arq,'\\newcommand{\\PARMptporcelY}{%d}\n',KP_v(2));
    fprintf(arq,'\\newcommand{\\PARMptporcelula}{%d}\n',KP_v(1)*KP_v(2));

    fprintf(2,'--- construcao da pre-base do nivel ----------------------\n');
    estagio = sprintf('nivel %02d (preb)',iv);
    
    %plota_malhas_hier2(malha_iv,...
    %parametros.min_nivel,passov,iv-passov,pasta_nivel,plota,estagio,'b');

    %determina o vetor K_v das posicoes das celulas criticas deste nivel
    K_v = determina_celulas_criticas(malha_iv,p_h,e_h,theta*tol); 
    nk_v = size(K_v,1);
    fprintf(2,'numero de celulas criticas antes da aproximacao %d\n',nk_v);
    fprintf(arq,'\\newcommand{\\PARMncelCrit}{%d}\n',nk_v);
    if nk_v == 0 
        fprintf(2,'criterio satisfeito antes do nivel %d\n',iv);
        break
    end
    
    [E_v,Q_v] = determina_elementos_relevantes(K_v,NC_v,LD,raio_abs_v,...
        parametros.tipobase,parametros.prebase_gorda);
     ne_v = size(E_v,1); %numero de elementos escolhidos neste nivel
    fprintf(2,'numero de elementos escolhidos %d\n',ne_v);
    fprintf(arq,'\\newcommand{\\PARMnelEsc}{%d}\n',ne_v);
    nr_v = size(Q_v,1); %numero de celulas relevantes neste nivel
    fprintf(2,'numero de celulas relevantes %d\n',nr_v);
    fprintf(arq,'\\newcommand{\\PARMncelRel}{%d}\n',nr_v);
    
    %constroi uma malha com as celulas relevantes e determina o vetor PR_v
    %que seleciona os pontos de amostragem relevantes
    fprintf(2,'contruindo a malha da pre-base\n');
    [PR_v,malha_v] = constroi_malha_da_prebase(Q_v,CD,LD,NC_v,p_h);
    
    %p_v sao os pontos de amostragem relevantes
    p_v = p_h(PR_v,:);
    np_v = size(p_v,1);
    
    fprintf(2,'%d pontos de amostragem relevantes\n',np_v);
    fprintf(arq,'\\newcommand{\\PARMnpPre}{%d}\n',np_v);
    
    %constroi a base com os elementos relevantes para este nivel
    base_v = constroi_base(malha_v,E_v,raio_abs_v,parametros.tipobase);
%    confere_pontos_da_base(malha_v,base_v,NP,PR_v,...
 %       parametros.tensorial,parametros.dobrado,p_v,p_ini);
    nb_v = size(base_v.elementos,2);
    fprintf(2,'tamanho da pre-base %d\n',nb_v);
    fprintf(arq,'\\newcommand{\\PARMdimApre}{%d}\n',nb_v);
    assert(nb_v == ne_v)
    %mostra_base(base_v)
    
    [A_v,M_v] = obtem_base_discreta(parametros,p_v,malha_v,base_v,...
        parametros.tensorial,KP_v,...
        NP,np_comp,PR_v);
    if plota.base
        plota_base(parametros,malha_v,NP,p_v,A_v,base_v,pasta_nivel,...
            plota,estagio);
    end

    figure(10 + plota.janela_base); clf;
    plota_centros(base_v,p_v,A_v,malha_v,NP,parametros,pasta_nivel,...
    plota,estagio)

% 
% plota_centros_hier2(base_v,p_v,A_v,malha_v,NP,...
%             parametros,parametros.min_nivel,passov,iv,inib_h,nb_h,...
%             pasta_nivel,plota,estagio)
        

    assert((parametros.metodo == 1) || (parametros.metodo == 3)...
         || (parametros.metodo == 4)) 
    
    %o problema eh aproximar o residuo e_h pela base A_v nos pontos p_v
    e_hv = e_h(PR_v); %residuo nos pontos p_v
    grava_tex = true;
    [b_v,alpha_v,es_v,ee_v,map_v,condM_v] = aproximacao(parametros,NP,...
    p_v,np_comp,CD,LD,base_v,A_v,e_hv,M_v,pasta_nivel,plota,...
    estagio,grava_tex);

    erro_v_L2 = norma_erro_L2(ee_v,np_comp);
    erro_v_Linf = norma_erro_Linf(ee_v);
    fprintf(2,'erro residual com a pre-base do nivel = %9.6f (L2) %9.6f (Linf)\n', ...
        erro_v_L2, erro_v_Linf);
    
    if parametros.adap_por_nivel
        fprintf(2,'--- aproximacao adaptativa 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;
        [ix_a,base_a,A_a,alpha_a] = ...
            escolhe_sub_base(parametros,base_v,A_v,alpha_v,ee_v);
        M_a = M_v(ix_a,ix_a);
        nb_a = size(A_a,2);
        fprintf(2,'tam base adaptativa %d\n',nb_a);
        if nb_a == 0
            ee_a = e_hv;
            erro_a_L2 = norma_erro_L2(ee_a,np_comp);
            erro_a_Linf = norma_erro_Linf(ee_a);
            fprintf(2,'erro com base reduzida e alphas originais %9.6f (L2) %9.6f (Linf)\n',...
                erro_a_L2, erro_a_Linf);
            arq = fopen(sprintf('%s/paramadapt.tex',pasta_nivel),'wt');
            fprintf(arq,'\\newcommand{\\PARMerroaHilb}{%s}\n',...
                converte_erro_tex(erro_a_L2));
            fprintf(arq,'\\newcommand{\\PARMerroaLinf}{%s}\n',...
                converte_erro_tex(erro_a_Linf));
            fprintf(arq,'\\newcommand{\\PARMtol}{%s}\n',...
                converte_erro_tex(tol));
            fprintf(arq,'\\newcommand{\\PARMnba}{%d}\n',nb_a);
            fclose(arq);
            break
        end
        figure(100); clf;
        plota_centros(base_a,p_v,A_a,malha_v,NP,parametros,pasta_nivel,...
            plota,estagio)
        % plota_centros_hier2(base_v,p_v,A_v,malha_v,NP,...
         %   parametros,parametros.min_nivel,passov,iv,inib_h,nb_h,...
          %  pasta_nivel,plota,estagio)

        %calcula residuo com coeficientes originais e base reduzida
        es_a = A_a*alpha_a;
        ee_a = e_hv - es_a;
        erro_a_L2 = norma_erro_L2(ee_a,np_comp);
        erro_a_Linf = norma_erro_Linf(ee_a);
        fprintf(2,'erro com base reduzida e alphas originais %9.6f (L2) %9.6f (Linf)\n',...
            erro_a_L2, erro_a_Linf);
        arq = fopen(sprintf('%s/paramadapt.tex',pasta_nivel),'wt');
        fprintf(arq,'\\newcommand{\\PARMerroaHilb}{%s}\n',...
            converte_erro_tex(erro_a_L2));
        fprintf(arq,'\\newcommand{\\PARMerroaLinf}{%s}\n',...
            converte_erro_tex(erro_a_Linf));
        fprintf(arq,'\\newcommand{\\PARMtol}{%s}\n',...
            converte_erro_tex(tol));
        fprintf(arq,'\\newcommand{\\PARMnba}{%d}\n',nb_a);
        
        %substitui a base pela base reduzida
        base_v = base_a;
        A_v = A_a;
        M_v = M_a;
        nb_v = nb_a;
        if (parametros.hier_adap_recalc)
            estagio = '(recalc)';
            pasta_nivel_recalc = cria_pasta(sprintf('%s/iv%02dr',...
                pasta_global,iv));
            grava_tex = true;
            [b_v,alpha_v,es_v,ee_v,map_v,condM_v] = aproximacao(parametros,NP,...
                p_v,np_comp,CD,LD,base_v,A_v,e_hv,M_v,pasta_nivel_recalc,...
                plota,estagio,grava_tex);            
            erro_v_L2 = norma_erro_L2(ee_v,np_comp);
            erro_v_Linf = norma_erro_Linf(ee_v);
            fprintf(2,'erro com base reduzida e alphas recalculados %9.6f (L2) %9.6f (Linf)\n',...
                erro_a_L2, erro_a_Linf);
            fprintf(arq,'\\newcommand{\\PARMerroarHilb}{%s}\n',...
                converte_erro_tex(erro_a_L2));
            fprintf(arq,'\\newcommand{\\PARMerroarLinf}{%s}\n',...
                converte_erro_tex(erro_a_Linf));
        else 
             alpha_v = alpha_a;
             es_v = es_a;
             ee_v = ee_a;
             erro_v_L2 = erro_a_L2;
             erro_v_Linf = erro_a_Linf;
        end
    end
    if erro_v_Linf <= tol
        fprintf(2,'criterio satisfeito no nivel %d\n',iv);
        break
    end
       
    fclose(arq);
    
    if iv == parametros.min_nivel
        ib_v = 1;
    else
        ib_v = inib_h(iv-passov) + tamb_h(iv-passov);
    end
    fb_v = ib_v + nb_v - 1;
    
    s_h(PR_v) = s_h(PR_v) + es_v;
    e_h(PR_v) = ee_v;
    
    A_h(:,ib_v:fb_v) = zeros(np_ini,nb_v);
    A_h(PR_v,ib_v:fb_v) = A_v;  
    alpha_h(ib_v:fb_v) = alpha_v;
   
    inib_h(iv) = ib_v;
    tamb_h(iv) = nb_v;
    nb_h = nb_h + nb_v
    fprintf(2,'tamanho da base hierarquica acumulada = %d\n', nb_h);
    
    malha_ant = malha_v;
    p_ant = p_v;
    PR_ant = PR_v;
    KP_ant =KP_v;
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