function principal_avalia_hier3_adap(parametros,plota)
%Encontra uma aproximacao de uma funcao dada com base tensorial
%hierarquica adaptativa.
%Essa versao usa interpolacao de medias em cada nivel

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

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;

%define a malha completa no nivel maximo
fprintf(2,'definindo a malha regular completa do nivel maximo\n');
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
if 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_comp = KP_fin .* NC_fin;
    p_comp = constroi_grade_de_amostras(CD,LD,NP_comp);
else
    NP_comp = [NaN NaN];
    p_comp = obtem_pontos_irregulares(parametros.n_pontos,CD,LD);
end
np_comp = size(p_comp,1);
fprintf(2,'numero de pontos de amostragem  = %d\n',np_comp);

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

%tamanho da malha inicial completa
NC_ini = tamanho_malha_hierarquica(parametros.min_nivel);

assert(max(mod(NC_fin,NC_ini)) == 0)

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-AH3-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,:));
if parametros.n_pontos == 0
    dir2 = sprintf('NX%05d-NY%05d-r%08.5f',NP_comp(1),NP_comp(2),raio_rel);
else
    dir2 = sprintf('NP%08d-r%08.5f',parametros.n_pontos,raio_rel);
end
pasta_global = cria_pasta(sprintf('%s/%s/%s',...
    pasta_raiz,dir1,dir2));

%constroi a matriz f_comp da funcao a aproximar
%f_comp(i,j) valor da funcao numero j ( = 1) no ponto de amostragem i
f_comp = discretiza_funcao(parametros,p_comp,CD,LD);
assert(size(f_comp,2) == 1);
nf = 1;
if plota.funcao && (parametros.n_pontos == 0)
    figure(9);
    for j = 1:nf
        plota_funcao(parametros,p_comp,f_comp(:,j),CD,LD,NP_comp,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);
if parametros.n_pontos == 0
    fprintf(arq_global,'\\newcommand{\\PARMamostrX}{%d}\n',NP_comp(1));
    fprintf(arq_global,'\\newcommand{\\PARMamostrY}{%d}\n',NP_comp(2));
end
fprintf(arq_global,'\\newcommand{\\PARMamostras}{%d}\n',np_comp);
fprintf(arq_global,'\\newcommand{\\PARMraio}{%s}\n',...
    converte_raio_tex(raio_rel));

tol = parametros.tol;
NC_ant = [NaN NaN];

%base_h eh a base hierarquica (multinivel) atual
%nb_h eh numero total de elementos na base hierarquica
%alpha_h eh o vetor de coeficientes para a base base_h
%inib_h(k) eh o comecco do nivel k em alpha_h, base_h
%tamb_h(k) eh o numero de elementos da base no nivel k
%malha_h(k) eh a malha com as celulas relevantes do nivel k
nb_h = 0;
alpha_h = zeros(nb_h,nf);
for iv = parametros.min_nivel:passov:parametros.max_nivel
    %disp('estamos no pause')
    %pause
    fprintf(2,'=== nivel %d ======================================\n',iv);
    estagio = sprintf('nivel %02d (preb)',iv);
    plota.janela_base = 20;
    pasta_nivel = cria_pasta(sprintf('%s/iv%02d',pasta_global,iv));
    
    %constroi a malha do nivel iv (malha_v)
    if iv == parametros.min_nivel
        %constroi a malha inicial
        NC_v = NC_ini;
        malha_v = constroi_malha_para_pontos(CD,LD,NC_v,p_comp);
    else
        %refina a malha anterior
        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));
        NC_v = NC_ant .* KF;
        malha_v = refina_malha(malha_ant,KF,p_comp);
    end
    
    nc_v = size(malha_v.celulas,2);
    malha_h(iv) = malha_v;
    fprintf(2,'numero de celulas na malha completa %d x %d = %d\n',...
        NC_v(1),NC_v(2),prod(NC_v));
    fprintf(2,'numero de celulas ocupadada = %d\n',nc_v);
    %A_h  sao osvalores de todos os elementos nos pontos de p
    %M_h eh a matriz do sistema MQ
    raio_abs_v = (LD./NC_v)*raio_rel;
    
    fprintf(2,'raio absoluto do nivel %d = [%8.5f %8.5f]\n',...
        iv,raio_abs_v(1),raio_abs_v(2));
    
    %obtem os centros das celulas e as medias celulares
    [E_v,p_v,f_v,PR_v] = obtem_centros_e_medias_celulares(malha_v,CD,LD,...
        p_comp,f_comp);
    
    NP_v = NC_v;
    np_v = prod(NP_v);
    assert(np_v == size(p_v,1))
    assert(size(f_v,1) == np_v)
    assert(size(f_v,2) == nf)
    
    if plota.funcao
        figure(6);
        for j = 1:nf
            plota_funcao(parametros,p_v,f_v(:,j),CD,LD,NP_v,NaN,NaN);
            title(sprintf('medias celulares'),'FontSize',20);
            hold off
            if plota.salva_eps
                saveas(gcf,sprintf('%s/mc%04d.eps',pasta_nivel,j));
            end
        end
       % disp('segundo pause')
        %pause
    end
    
    %recalcula a matriz da base A_h para o conjunto de pontos p_v
    A_h = sparse(np_v,0);
    fprintf(2,'---------------------------------\n');
    fprintf(2,'recalculando a aproximacao anterior\n');
    for jv = parametros.min_nivel:passov:iv-passov
        fprintf(2,'-- expandindo nivel %d para pontos do nivel %d\n',jv,iv);
        ini_j = inib_h(jv);
        fim_j = ini_j + tamb_h(jv) - 1;
        base_j.elementos = base_h.elementos(ini_j:fim_j);
        base_j.uniforme = true;
        malha_j = malha_h(jv);
        NC_j = malha_j.NC;
        KP_j = NP_v ./ NC_j;
        fprintf(2,'obtendo base do nivel %d no pontos correntes\n', jv);
        fprintf(2,'raio_abs = %9.6f %9.6f \n', base_j.elementos(1).raio_abs);
        fprintf(2,'pontos por celula = %d %d\n', KP_j);
        [A0_j,A_j] = obtem_matriz_A(parametros,p_v,malha_j,base_j,...
            parametros.tensorial,KP_j,NP_v,PR_v);
        
        temp = plota.salva_eps;
        plota.salva_eps = false;
        estagio_j = sprintf('%s:%d',estagio,j); 
%         plota_base(parametros,malha_j,NP_v,p_v,A_j,base_j,pasta_nivel,...
%             plota,estagio_j);
        plota.salva_eps = temp;
        
        A_h(:, ini_j:fim_j) = A_j;
        assert(issparse(A_h))
    end
    fprintf(2,'fim do recalculo da aproximacao anterior\n');
    fprintf(2,'---------------------------------\n');
    
    %calcula a aproximaçao corrente s_v nos pontos p_v
    if iv == parametros.min_nivel
        s_v = zeros(np_v,nf);
    else
        s_v = A_h(:,1:nb_h)*alpha_h(1:nb_h, :);
    end
    assert(size(s_v,1) == np_v)
    assert(size(s_v,2) == nf)
    
    if plota.aproximacao
        figure(7);
        for j = 1:nf
            plota_funcao(parametros,p_v,s_v(:,j),CD,LD,NP_v,NaN,NaN);
            title(sprintf('aproximacao anterior'),'FontSize',20);
            hold off
            if plota.salva_eps
                saveas(gcf,sprintf('%s/sv%04d.eps',pasta_nivel,j));
            end
        end
%         disp('terceiro pause')
%         pause
    end
    
    %calcula o residuo e_v das medias celulars f_v relativo a aprox s_v
    e_v = f_v - s_v;
    
    if plota.erro
        figure(8);
        for j = 1:nf
            plota_funcao(parametros,p_v,e_v(:,j),CD,LD,NP_v,NaN,NaN);
            title(sprintf('residuo corrente'),'FontSize',20);
            hold off
            if plota.salva_eps
                saveas(gcf,sprintf('%s/ev%04d.eps',pasta_nivel,j));
            end
        end
%         disp('quarto pause')
%         pause
    end
    
    %constroi a prebase do nivel iv
    base_v = constroi_base(malha_v,E_v,raio_abs_v,parametros.tipobase);
    nb_v = size(base_v.elementos,2);
    fprintf(2,'tamanho da pre-base deste nivel %d\n',nb_v);
    KP_v = [1, 1];
    [A0_comp,A_v] = obtem_matriz_A(parametros,p_v,malha_v,base_v,...
        parametros.tensorial,KP_v,NP_v,PR_v);
    
    if plota.base
        plota_base(parametros,malha_v,NP_v,p_v,A_v,base_v,pasta_nivel,...
            plota,estagio);
    end
    
    %interpola o residuo e_v usando a base do nivel iv
    assert(parametros.metodo == 2)
    grava_tex = true;
    M_v = A_v;
    assert(nb_v == np_v)
    [b_v,alpha_v,es_v,ee_v,map,condM] = ...
        aproximacao(parametros,NP_v,p_v,np_v,CD,LD,base_v,A_v,e_v,M_v,...
        pasta_nivel,plota,estagio,grava_tex);
    
    %     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)
    
    erro_v_L2 = norma_erro_L2(ee_v,np_v);
    erro_v_Linf = norma_erro_Linf(ee_v);
    fprintf(2,'erro residual com a base  do nivel = %9.6f (L2) %9.6f (Linf)\n', ...
        erro_v_L2, erro_v_Linf);
    
    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{\\PARMdimAnivel}{%d}\n',nb_v);
    
    %disp('quinto pause')
    %pause
        
    if parametros.adap_por_nivel
        fprintf(2,'--- aproximacao adaptativa no 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,alpha_a] = ...
            escolhe_sub_base(parametros,base_v,A_v,M_v,...
            alpha_v,ee_v);
        nb_a = size(A_a,2);
        fprintf(2,'tam base adaptativa do nivel = %d\n',nb_a);
        
        %         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
        es_a = A_a*alpha_a;
        ee_a = f_v - es_a;
        erro_a_L2 = norma_erro_L2(ee_a,np_v);
        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_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);
        fclose(arq_a);
        
        %substitui a prebase pela base reduzida
        base_v = base_a;
        nb_v = nb_a;
        A_v = A_a;
        M_v = M_a;
        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
    fclose(arq);
    assert(size(base_v.elementos,2) == nb_v)
    
    %plota a aproximacao acumulada incluindo este nivel
    s_v = s_v + es_v;
    e_v = ee_v;
    if plota.aproximacao
        assert(nf == 1)
        plota.janela_base = 40;
        estagio = '(global)';
        nomef = 'funcao';
        pasta = pasta_global;
        figure(5+plota.janela_base);
        plota_funcao(parametros,p_v,s_v,CD,LD,NP_v,NaN,NaN);
        title(sprintf('aproximacao da %s %s',...
            nomef,estagio),'FontSize',20);
        hold off
        if plota.salva_eps
            saveas(gcf,sprintf('%s/s_v.eps',pasta));
        end
        
        figure(6+plota.janela_base);
        if plota.erro_curva_nivel
            plot_contour_funcao(parametros,p_v,e_v,CD,LD,NP_v);
        else
            emax = max([plota.emax, max(abs(e_v))]);
            plota_funcao(parametros,p_v,e_v,CD,LD,NP_v,-emax,emax);
        end
        title(sprintf('erro da %s %s',...
            nomef,estagio),'FontSize',20);
        hold off
        if plota.salva_eps
            saveas(gcf,sprintf('%s/e_v.eps',pasta));
        end
    end
    
    %acrescenta elementos do nivel iv na base multinivel
    if nb_h == 0
        base_h = base_v;
        alpha_h = alpha_v;
        inib_h(iv) = 1;
    else 
        ini_v = nb_h + 1;
        fim_v = ini_v + nb_v - 1;
        base_h.elementos(1,ini_v:fim_v) = base_v.elementos(1,:);
        %alpha_h(1,ini_v:fim_v) = alpha_v(1,:);
         alpha_h(ini_v:fim_v,1) = alpha_v(:,1);
        base_h.uniforme = false;
        inib_h(iv) = ini_v;
    end
    tamb_h(iv) = nb_v;
    nb_h = nb_h + nb_v;
    
    malha_ant = malha_v;
    NC_ant = NC_v;
end
fprintf(arq_global,'\\newcommand{\\PARMnbh}{%d}\n',nb_h);
fclose(arq_global);
profile viewer
%deveria tentar achar a aproximacao e/ou plotar com a base inteira