/***************************************************************************
 *   Copyright (C) 2008 by Douglas Castro   *
 *   douglas@localhost.localdomain   *
 *                                                                         *
 ***************************************************************************/

#include "amr.h"

void mata_mato(No *r)
{
  if(r->fil[0]->fil[0]==NULL)
  {
    free(r->fil[0]); r->fil[0]=NULL;
    free(r->fil[1]); r->fil[1]=NULL;
    return;
  }
  mata_mato(r->fil[0]);
  mata_mato(r->fil[1]);
}

void troca( int *a, int *b)
{
  int temp = *a;
  *a = *b;
  *b = temp;
}

void compara_folhas(No *r, double vet[], int *cont,  double *norm)
{
  assert((r->fil[0] == NULL) == (r->fil[1] == NULL));
  if(r->fil[0] == NULL)
  {
    int c = *cont;
    double db = r->fval;
    double nn = *norm;
    double diff = fabs(vet[c]-db);
    if(diff>nn){ *norm = diff;}
    *cont = c + 1;
    return;
  }
  compara_folhas(r->fil[0], vet, cont, norm);
  compara_folhas(r->fil[1], vet, cont, norm);
}

void amr_analise_haar(No *r, int niv, int prof, double eps, bool_t sim)
{
  
  assert((r->fil[0] == NULL) == (r->fil[1] == NULL));
  
  if (r->fil[0] == NULL) 
  { /* Não há o que fazer: */
    //fprintf(stderr,"niv = %d  prof = %d \n",niv,prof);
    return;
  }
  
  int qual;
  for (qual = 0; qual <= 1; qual++)
  { 
    amr_analise_haar(r->fil[qual], niv, prof+1, eps, sim);
  }
  
  // lembrando: |C^k| = |C^{k-1}|/2
  //            |f^k| = (|f^{k+1}_e| + |f^{k+1}_d|)/2
  
  double restricao = ( r->fil[0]->fval + r->fil[1]->fval )/2.0;
  
  double detalhe = r->fil[1]->fval - restricao;//r->fval;
  
  r->fil[0]->fval = detalhe;
  r->fil[0]->fv[1] = detalhe;
  r->fil[1]->fval = restricao;
  r->fil[1]->fv[1] = restricao;
  r->fval = restricao;
  r->fv[1] = r->fval;

  
  /* Elimina os filhos, se sao nao-nulos mas ambos superfluos e tem netos nulos:  */
  /* truncamento dos detalhes, se necessario */
  if(sim)
  {// se sim == TRUE, entra aqui
    int pot = prof+1 - niv;
    double eps_k = eps*pow(2.0, pot);
    if((fabs(detalhe) < eps_k) && ((r->fil[0]->fil[0]==NULL) && (r->fil[1]->fil[0]==NULL)))
    { /* Elimina os dois filhos: */
      free(r->fil[0]); r->fil[0] = NULL;
      free(r->fil[1]); r->fil[1] = NULL;
    }
  }
  //return;
}

void amr_sintese_haar(No *r, int prof, int niv)
{
  assert((r->fil[0] == NULL) == (r->fil[1] == NULL));
  if (prof == (niv-1)) 
  { /* Não há o que fazer: */
    //fprintf(stderr,"prof = %d niv = %d\n",prof, niv);
    return;
  }
  //fprintf(stderr,"%*s prof = %d\n",2*prof,"",prof);
  double detalhe = r->fil[0] == NULL ? 0.0 : r->fil[0]->fval;
  double previsao = r->fil[0] == NULL ? r->fval : r->fil[1]->fval;

  if (r->fil[0] == NULL) 
  { 
    No *e = malloc(sizeof(No));
    No *d = malloc(sizeof(No));
    r->fil[0] = e;
    r->fil[0]->fil[0] = NULL;
    r->fil[0]->fil[1] = NULL;
    r->fil[1] = d;
    r->fil[1]->fil[0] = NULL;
    r->fil[1]->fil[1] = NULL;
  }
  //fprintf(stderr,"problema aqui oh %d \n", profund);
  r->fil[0]->fval = previsao - detalhe;
  r->fil[0]->fv[1] = previsao - detalhe;
  r->fil[1]->fval = previsao + detalhe;
  r->fil[1]->fv[1] = previsao + detalhe;
  
  int qual;
  for (qual = 0; qual <= 1; qual++)
  { 
    amr_sintese_haar(r->fil[qual], prof+1, niv);
  }
}

void amr_analise_harten(int d, bool_t op, VNo pac[], int niv, double xmin[], double xmax[], int prof, Preditor pred, double eps, int lim, bool_t sim)    
{
  int k = 5;
  if(op){k=9;} 
  int ctr = ipow(k,d)/2;
  No *r = pac[ctr].p;
  assert((r->fil[0] == NULL) == (r->fil[1] == NULL));
  
  if(pac[ctr].p->fil[0] == NULL) 
  { /* Não há o que fazer: */
    //fprintf(stderr,"niv = %d  prof = %d lim = %d \n",niv,prof,lim);
    return;
  }
  if(lim < prof) 
  { /* Não há o que fazer: */
    //fprintf(stderr,"niv = %d  prof = %d lim = %d \n",niv,prof,lim);
    return;
  }
  /* Decide o eixo {ex} perpendicular ao corte: */
  int ex = prof % d;
  
  int qual;
  for (qual = 0; qual <= 1; qual++)
  { 
    int tam = ipow(k,d);
    VNo pacf[tam];
    double xminf[d], xmaxf[d];
    /* Constrói o pacote do filho {qual}: */
    divide_pacote_geral(d, pac, ex, pred, qual, pacf, op);
    limites_celula( d, ex, qual, xmin, xmax, xminf, xmaxf);
    amr_analise_harten(d, op, pacf, niv, xminf, xmaxf, prof+1, pred, eps, lim, sim);
  }
  // lembrando: |C^k| = |C^{k-1}|/2
  //            |f^k| = (|f^{k+1}_e| + |f^{k+1}_d|)/2
  if( prof == lim )
  { /* tentativa de fazer analise um nivel por vez*/
    //fprintf(stderr,"%*s niv = %d  prof = %d lim = %d \n",2*prof,"",niv,prof,lim);
    
    double restricao = ( pac[ctr].p->fil[0]->fval + pac[ctr].p->fil[1]->fval )/2.0;
  
    int sz = op ? 5 : 3;
    int tp = sz/2, t = 0, s;
    double fi[sz];
    for(s=-tp;s<=tp;s++)
    {
      fi[t] = pac[ctr + s].fv;
      t = t+1;
    }
    
    double fprev = pred(op, fi);
  
    double detalhe = pac[ctr].p->fil[1]->fval - fprev;
    //double diff = pac[ctr].p->fil[0]->fval - pac[ctr].p->fil[1]->fval;
    
    /* como indicado em amr (funcao), eventualmente precisa ser comentado */
    //serao alteradas apenas os valores dos filhos das celulas de profundidade lim
    assert((pac[ctr].p->fil[0] == NULL) == (pac[ctr].p->fil[1] == NULL));
    //if(pac[ctr].p->fil[0] != NULL)
    //{ 
      pac[ctr].p->fil[1]->fval = detalhe;
      pac[ctr].p->fil[1]->fv[1] = detalhe;
      pac[ctr].p->fil[0]->fval = restricao;
      pac[ctr].p->fil[0]->fv[1] = restricao;
    //}
    /* ate aqui */
    /* truncamento dos detalhes, se necessario */
    assert((pac[ctr].p->fil[0]->fil[0] == NULL) == (pac[ctr].p->fil[0]->fil[1] == NULL));
    assert((pac[ctr].p->fil[1]->fil[0] == NULL) == (pac[ctr].p->fil[1]->fil[1] == NULL));
    if(sim)
    { // se sim == TRUE, entra aqui 
      // o parametro de truncamento para medias celulares nao e constante
      // depende da profundidade na malha
      double pot = ( prof+1 - niv ); 
      double eps_l = eps*pow(2.0, pot);
      if((fabs(detalhe) < eps_l) && ((pac[ctr].p->fil[0]->fil[0]==NULL) && (pac[ctr].p->fil[1]->fil[0]==NULL)))
      { /* Elimina os dois filhos: */
        free(pac[ctr].p->fil[0]); pac[ctr].p->fil[0] = NULL;
        free(pac[ctr].p->fil[1]); pac[ctr].p->fil[1] = NULL;
      }
    }
    return;
  }
}

void mostra_pac(FILE *arq, VNo pac[], int k)
{ int i0, i1;
  for (i1 = k-1; i1 >= 0; i1--)
  { for (i0 = 0; i0 < k; i0++)
    { int i = k*i1 + i0;
      fprintf(stderr, " %14.9f%s", pac[i].fv, (pac[i].p == NULL ? " " : "*"));
    }
    fprintf(stderr, "\n");
  }
}

void amr_sintese_harten(int d,  bool_t op, VNo pac[], double xmin[], double xmax[], int prof, Preditor pred, int niv, int lim)
{
  int  k = 5; // posteriormente k pode se tornar argumento da funcao
  if(op){k=9;}
  int tm = pow(k,d);
  int ctr = (tm)/2;
  
  if(prof == (niv-1)) 
  { /* Não há o que fazer: */
    //fprintf(stderr,"niv = %d  prof = %d lim = %d \n",niv,prof,lim);
    pac[ctr].p->fil[0] = NULL;
    pac[ctr].p->fil[1] = NULL;
    return;
  }
  //fprintf(stderr,"niv = %d  prof = %d lim = %d \n",niv,prof,lim);
  
  // cria uma arvore, sera usada se o filho for nulo
  assert((pac[ctr].p->fil[0] == NULL) == (pac[ctr].p->fil[1] == NULL));
  if(pac[ctr].p->fil[0] == NULL)
  { //return;
    pac[ctr].p->fil[0] = malloc(sizeof(No));
    pac[ctr].p->fil[0]->fval = 0.0;
    pac[ctr].p->fil[0]->fv[1] = 0.0;
    pac[ctr].p->fil[0]->fil[0] = NULL;
    pac[ctr].p->fil[0]->fil[1] = NULL;
    pac[ctr].p->fil[1] = malloc(sizeof(No));
    pac[ctr].p->fil[1]->fval = pac[ctr].p->fval;
    pac[ctr].p->fil[1]->fv[1] = pac[ctr].p->fval;
    pac[ctr].p->fil[1]->fil[0] = NULL;
    pac[ctr].p->fil[1]->fil[1] = NULL;
  }
  /* Decide o eixo {ex} perpendicular ao corte: */
  int ex = prof % d;

  // a sintese sera feita a partir de onde a analise parou
    // lim e 
  if(prof == lim )///?????????
  {
    int sz = op ? 5 : 3;
    int tp = sz/2, t = 0, s;
    double fi[sz];
    for(s=-tp;s<=tp;s++)
    {
      fi[t] = pac[ctr + s].fv;
      t = t+1;
    }
    
    double fprev = pred(op, fi);
    
    double detalhe = pac[ctr].p->fil[1]->fval;
  
    double dir = fprev + detalhe;
    pac[ctr].p->fil[1]->fval = dir;
    pac[ctr].p->fil[1]->fv[1] = dir;
    double esq = 2.0*pac[ctr].p->fil[1]->fval - dir;
    pac[ctr].p->fil[0]->fval = esq;
    pac[ctr].p->fil[0]->fv[1] = esq;
    return;
  }
  
  int qual;
  for (qual = 0; qual <= 1; qual++)
  { 
    VNo pacf[k*k];
    double xminf[d], xmaxf[d];
    /* Constrói o pacote do filho {qual}: */
    divide_pacote_geral(d, pac, ex, pred, qual, pacf, op);
    limites_celula( d, ex, qual, xmin, xmax, xminf, xmaxf);
    amr_sintese_harten(d, op, pacf, xminf, xmaxf, prof+1, pred, niv, lim);
  }
}

void analise_multiresolucao(int d, bool_t op, VNo pac[], int niv, double xmin[], double xmax[], int prof, Preditor pred, double eps, bool_t sim, bool_t haar)
{
  if(haar)
  {
    int ord = 5;
    if(op){ord=9;}
    int ctr = ipow(ord,d)/2;
    No *r = pac[ctr].p;
    amr_analise_haar( r, niv, prof, eps, sim);
  }
  else
  {
    int lim;
    for(lim = niv-1;lim >= 0;lim--) /* pq niv-2 ? */
    {
      fprintf(stderr, "lim = %d\n", lim);
      amr_analise_harten( d, op, pac, niv, xmin, xmax, prof, pred, eps, lim, sim);
    }
  }
}

void sintese_multiresolucao(int d, bool_t op, VNo pac[], int niv, double xmin[], double xmax[], int prof, Preditor pred, double eps, bool_t haar)
{
  
  if(haar)
  {
    int ord = 5;
    int ctr = ipow(ord,d)/2;
    No *r = pac[ctr].p;
    amr_sintese_haar( r, prof, niv);
  }
  else
  {
    int lim;
    for(lim = 0; lim<=niv-2; lim++)
    {
      amr_sintese_harten( d, op, pac, xmin, xmax, prof, pred, niv, lim);
    }
  }
}


void amr(int d, bool_t op, VNo pac[], int niv, double xmin[], double xmax[], int prof, Preditor pred, double eps, bool_t sim, bool_t haar, bool_t recs)
{
  int k = 5;
  if(op){k=9;}
  int kc = ipow(k,d)/2;
  /* numero de celulas folha */
  int ncf = 0;
  conta_folhas( pac[kc].p, &ncf);
  fprintf(stderr,"numero de folhas malha regular = %d\n", ncf);
  /* Guarda estado inicial para verificar a solucao analitica */
  double mecf[ncf];
  int pos = 0;
  copia_folhas(pac[kc].p, mecf, &pos);
  /* teste, somente d==1 */
  auto void escreveteste(FILE *arq, double vet[], double xmin[], double dx, int ncf);
  void escreveteste(FILE *arq, double vet[], double xmin[], double dx ,int ncf)
  {
    int i;
    for(i=0; i<ncf;i++)
    {
      fprintf(arq, "%10.6f  %10.6f\n", xmin[0]+i*dx, vet[i]);
    }
  }
  //double dx = fabs(xmax[0]-xmin[0])/ncf;
  //FILE *arq = fopen("med.txt", "w");
  //escreveteste(arq, mecf, xmin, dx, ncf);
  //fclose(arq);
  //fprintf(stderr,"escreveu med.txt\n");

  analise_multiresolucao( d, op, pac, niv, xmin, xmax, 0, pred, eps, sim, haar);
  fprintf(stderr,"fez analise\n");
 
  //FILE *disco = fopen("solucao.txt", "w");
 // escreve_em_arquivo( d, 0, disco, pac, xmin, xmax, pred, 1, FALSE);
  //fclose(disco);

  /* Inicia gnuplot */
  //gnuplot_ctrl *h3;
  //h3 = gnuplot_init();
  //exportando( d, h3);
  //gnuplot_close(h3) ;
  
  /* conta folhas */
  int qtd = 0;
  conta_folhas( pac[kc].p, &qtd);
  fprintf(stderr,"quantidade de folhas depois da poda = %d\n", qtd);
  double porc = 100.0 - 100.0*qtd/ncf;
  fprintf(stderr,"compressao = %f\n", porc);
  
//   MakeTree(OUT_PREFIX2, pac[kc].p, niv,0);
  //DoMesh(OUT_PREFIX2, pac[kc].p);

  if(recs)
  {
    sintese_multiresolucao( d, op, pac, niv, xmin, xmax, 0, pred, eps, haar);
    fprintf(stderr,"fez reconstrucao\n");
    //qtd = 0;
    //conta_folhas( pac[kc].p, &qtd);
    //fprintf(stderr,"quantidade de folhas depois da rec = %d\n", qtd);
  
    /* exporta tipo eps a reconstrucao das medias */
    //FILE *disco = fopen("solucao.txt", "w");
    //escreve_em_arquivo( d, 0, disco, pac, xmin, xmax, pred, 1, FALSE);
    //fclose(disco);
    //fprintf(stderr,"escreveu solucao.txt\n");

    /* Inicia gnuplot */
    //gnuplot_ctrl *h2;
    //h2 = gnuplot_init();
    //exportando( d, h2);
    //gnuplot_close(h2) ;
    /* fim exportacao */
    /*
    fprintf(stderr,"copiar folhas reconstruidas\n");
    double mrcf[ncf];
    pos = 0;
    copia_folhas(pac[kc].p, mrcf, &pos);
    fprintf(stderr,"copiou folhas reconstruidas\n");
  
    FILE *arq2 = fopen("rec.txt", "w");
    escreveteste(arq2, mrcf, xmin, dx, ncf);
    fclose(arq2);
    fprintf(stderr,"escreveu rec.txt\n");
  */
    /* calculo da norma L1 do erro */
    double cvol = 1.0;
    int i;
    for(i=0;i<d;i++)
    {
      int qce = qtd_cel_eixo( d, i, niv-1);
    //fprintf(stderr,"cel por eixo = %d\n",qce);
      cvol *= fabs(xmax[i]-xmin[i])/qce;
    }
    fprintf(stderr,"cvol = %f\n",cvol);
    //double erro = erro_norma_l1(mecf, mrcf, ncf, cvol);
    //fprintf(stderr,"erro norma l1 = %16.13f\n", erro);
    double soma = 0.0;
    int pos = 0;
    erro_norma_l1_v2(mecf, pac[kc].p, &pos, &soma);
    double norma = soma*cvol;
    fprintf(stderr,"erro norma l1 = %.20f\n", norma);
  }
  
}

void reconstroi(int d, bool_t op, VNo pac[], int niv, double xmin[], double xmax[], int prof, Preditor pred, bool_t haar)
{
  int  k = op ? 9 : 5; 
  int tm = ipow(k,d);
  int ctr = tm/2;
  
  if(prof == (niv-1)) 
  { /* Não há o que fazer: */
    //fprintf(stderr,"niv = %d  prof = %d lim = %d \n",niv,prof,lim);
    pac[ctr].p->fil[0] = NULL;
    pac[ctr].p->fil[1] = NULL;
    return;
  }
  
  /* Decide o eixo {ex} perpendicular ao corte: */
  int ex = prof % d;
  // cria uma arvore, sera usada se o filho for nulo
  assert((pac[ctr].p->fil[0] == NULL) == (pac[ctr].p->fil[1] == NULL));
  if(pac[ctr].p->fil[0] == NULL)
  { //return;
    if(haar)
    {
      /* conferir */
      double fprev = pac[ctr].fv;
      
      pac[ctr].p->fil[1] = malloc(sizeof(No));
      pac[ctr].p->fil[1]->fval = fprev;
      pac[ctr].p->fil[1]->fv[0] = fprev;
      pac[ctr].p->fil[1]->fv[1] = fprev;
      pac[ctr].p->fil[1]->fil[0] = NULL;
      pac[ctr].p->fil[1]->fil[1] = NULL;
      pac[ctr].p->fil[0] = malloc(sizeof(No));
      pac[ctr].p->fil[0]->fval = 2.0*pac[ctr].p->fval - fprev;
      pac[ctr].p->fil[0]->fv[0] = 2.0*pac[ctr].p->fval - fprev;
      pac[ctr].p->fil[0]->fv[1] = 2.0*pac[ctr].p->fval - fprev;
      pac[ctr].p->fil[0]->fil[0] = NULL;
      pac[ctr].p->fil[0]->fil[1] = NULL;
    }
    else
    {
      int sz = op ? 5 : 3;
      int tp = sz/2, t = 0, s;
      double fi[sz];
      for(s=-tp;s<=tp;s++)
      {
        fi[t] = pac[ctr + s].fv;
        t = t+1;
      }
    
      double fprev = pred(op, fi);
      
      pac[ctr].p->fil[1] = malloc(sizeof(No));
      pac[ctr].p->fil[1]->fval = fprev;
      pac[ctr].p->fil[1]->fv[0] = fprev;
      pac[ctr].p->fil[1]->fv[1] = fprev;
      pac[ctr].p->fil[1]->fil[0] = NULL;
      pac[ctr].p->fil[1]->fil[1] = NULL;
      pac[ctr].p->fil[0] = malloc(sizeof(No));
      pac[ctr].p->fil[0]->fval = 2.0*pac[ctr].p->fval - fprev;
      pac[ctr].p->fil[0]->fv[0] = 2.0*pac[ctr].p->fval - fprev;
      pac[ctr].p->fil[0]->fv[1] = 2.0*pac[ctr].p->fval - fprev;
      pac[ctr].p->fil[0]->fil[0] = NULL;
      pac[ctr].p->fil[0]->fil[1] = NULL;
    }
  } 
  
  int qual;
  for (qual = 0; qual <= 1; qual++)
  { 
    int tam = pow(k,d);
    VNo pacf[tam];
    double xminf[d], xmaxf[d];
    /* Constrói o pacote do filho {qual}: */
    divide_pacote_geral(d, pac, ex, pred, qual, pacf, op);
    limites_celula( d, ex, qual, xmin, xmax, xminf, xmaxf);
    reconstroi( d, op, pacf, niv, xminf, xmaxf, prof+1, pred, haar);
  }
}

/* ----------------------------- poda --------------------------------*/

void poda_haar(No *r, int niv, int prof, double eps, bool_t sim)
{
  
  assert((r->fil[0] == NULL) == (r->fil[1] == NULL));
  
  if (r->fil[0] == NULL) 
  { /* Não há o que fazer: */
    //fprintf(stderr,"niv = %d  prof = %d \n",niv,prof);
    return;
  }
  
  int qual;
  for (qual = 0; qual <= 1; qual++)
  { 
    poda_haar(r->fil[qual], niv, prof+1, eps, sim);
  }
  
  // lembrando: |C^k| = |C^{k-1}|/2
  //            |f^k| = (|f^{k+1}_e| + |f^{k+1}_d|)/2

  double detalhe = r->fil[1]->fval - r->fval; //( pac->fil[0]->fval - pac->fil[1]->fval )/2.0;
  
  /* Elimina os filhos, se sao nao-nulos mas ambos superfluos e tem netos nulos:  */
  /* truncamento dos detalhes, se necessario */
  if(sim)
  {// se sim == TRUE, entra aqui
    int pot = prof+1 - niv;
    double eps_k = eps*pow(2.0, pot);
    if((fabs(detalhe) < eps_k) && ((r->fil[0]->fil[0]==NULL) && (r->fil[1]->fil[0]==NULL)))
    { /* Elimina os dois filhos: */
      free(r->fil[0]); r->fil[0] = NULL;
      free(r->fil[1]); r->fil[1] = NULL;
    }
  }
}

void refina_haar(No *r, int prof, int niv)
{
  assert((r->fil[0] == NULL) == (r->fil[1] == NULL));
  if (prof == niv-1) 
  { /* Não há o que fazer: */
    //fprintf(stderr,"prof = %d niv = %d\n",prof, niv);
    return;
  }

  double previsao = r->fil[0] == NULL ? r->fval : r->fil[1]->fval;
  No *e;
  No *d; 
  if (r->fil[0] == NULL) 
  { 
    e = malloc(sizeof(No));
    d = malloc(sizeof(No));
    r->fil[0] = e;
    r->fil[0]->fil[0] = NULL;
    r->fil[0]->fil[1] = NULL;
    r->fil[1] = d;
    r->fil[1]->fil[0] = NULL;
    r->fil[1]->fil[1] = NULL;
    r->fil[0]->fval = previsao;
    r->fil[0]->fv[1] = previsao;
    r->fil[1]->fval = previsao;
    r->fil[1]->fv[1] = previsao;
     //fprintf(stderr,"problema aqui oh %d \n", profund);
    return;
  }
  
  int qual;
  for (qual = 0; qual <= 1; qual++)
  { 
    refina_haar(r->fil[qual], prof+1, niv);
  }
}

void poda_harten(int d, bool_t op, VNo pac[], int niv, double xmin[], double xmax[], int prof, Preditor pred, double eps, int lim, bool_t sim)    
{
  int k = op ? 9 : 5; 
  int ctr = ipow(k,d)/2;
  No *r = pac[ctr].p;
  assert((r->fil[0] == NULL) == (r->fil[1] == NULL));
  
  if(r->fil[0] == NULL) 
  { /* Não há o que fazer: */
    return;
  }
  if(lim < prof) 
  { /* Não há o que fazer: */
    return;
  }
  /* Decide o eixo {ex} perpendicular ao corte: */
  int ex = prof % d;
  
  int qual;
  for (qual = 0; qual <= 1; qual++)
  { 
    int tam = ipow(k,d);
    VNo pacf[tam];
    double xminf[d], xmaxf[d];
    /* Constrói o pacote do filho {qual}: */
    divide_pacote_geral(d, pac, ex, pred, qual, pacf, op);
    limites_celula( d, ex, qual, xmin, xmax, xminf, xmaxf);
    poda_harten(d, op, pacf, niv, xminf, xmaxf, prof+1, pred, eps, lim, sim);
  }
  // lembrando: |C^k| = |C^{k-1}|/2
  //            |f^k| = (|f^{k+1}_e| + |f^{k+1}_d|)/2
  if( prof == lim )
  { /* tentativa de fazer analise um nivel por vez*/
    //fprintf(stderr,"%*s niv = %d  prof = %d lim = %d \n",2*prof,"",niv,prof,lim);

    int sz = op ? 5 : 3;
    int tp = sz/2, t = 0, s;
    int corr = ipow(k,ex);
    double fi[sz];
    for(s=-tp;s<=tp;s++)
    {
      fi[t] = pac[ctr + s*corr].p == NULL ? pac[ctr + s*corr].fv : pac[ctr + s*corr].p->fval;
      t = t+1;
    }
    
    double fprev = pred(op, fi);
  
    double detalhe = pac[ctr].p->fil[1]->fval - fprev;
    //double diff = pac[ctr].p->fil[0]->fval - pac[ctr].p->fil[1]->fval;
    
    /* como indicado em amr (funcao), eventualmente precisa ser comentado */
    //serao alteradas apenas os valores dos filhos das celulas de profundidade lim
    assert((pac[ctr].p->fil[0] == NULL) == (pac[ctr].p->fil[1] == NULL));

    /* truncamento dos detalhes, se necessario */
    assert((pac[ctr].p->fil[0]->fil[0] == NULL) == (pac[ctr].p->fil[0]->fil[1] == NULL));
    assert((pac[ctr].p->fil[1]->fil[0] == NULL) == (pac[ctr].p->fil[1]->fil[1] == NULL));
    if(sim)
    { // se sim == TRUE, entra aqui 
      // o parametro de truncamento para medias celulares nao e constante
      // depende da profundidade na malha
      double pot = ( prof+1 - niv ); 
      double eps_l = eps*pow(2.0, pot);
      if((fabs(detalhe) < eps_l) && ((pac[ctr].p->fil[0]->fil[0]==NULL) && (pac[ctr].p->fil[1]->fil[0]==NULL)))
      { /* Elimina os dois filhos: */
        pac[ctr].p->fval = 0.5*(pac[ctr].p->fil[0]->fval + pac[ctr].p->fil[1]->fval);
        free(pac[ctr].p->fil[0]); pac[ctr].p->fil[0] = NULL;
        free(pac[ctr].p->fil[1]); pac[ctr].p->fil[1] = NULL;
      }
    }
    return;
  }
}

void poda_hartenv2(int d, bool_t op, VNo pac[], int niv, double xmin[], double xmax[], int prof, Preditor pred, double eps, int lim, bool_t sim)    
{
  int k = op ? 9 : 5; 
  int ctr = ipow(k,d)/2;
  No *r = pac[ctr].p;
  assert((r->fil[0] == NULL) == (r->fil[1] == NULL));
  
  if(r->fil[0] == NULL) 
  { /* Não há o que fazer: */
    return;
  }
  if(lim < prof) 
  { /* Não há o que fazer: */
    return;
  }
  /* Decide o eixo {ex} perpendicular ao corte: */
  int ex = prof % d;
  
  int qual;
  for (qual = 0; qual <= 1; qual++)
  { 
    int tam = ipow(k,d);
    VNo pacf[tam];
    double xminf[d], xmaxf[d];
    /* Constrói o pacote do filho {qual}: */
    divide_pacote_geral(d, pac, ex, pred, qual, pacf, op);
    limites_celula( d, ex, qual, xmin, xmax, xminf, xmaxf);
    poda_hartenv2(d, op, pacf, niv, xminf, xmaxf, prof+1, pred, eps, lim, sim);
  }
  // lembrando: |C^k| = |C^{k-1}|/2
  //            |f^k| = (|f^{k+1}_e| + |f^{k+1}_d|)/2
  if( prof == lim )
  { /* tentativa de fazer analise um nivel por vez*/
    //fprintf(stderr,"%*s niv = %d  prof = %d lim = %d \n",2*prof,"",niv,prof,lim);

    int sz = op ? 5 : 3;
    int tp = sz/2, t = 0, s;
    int corr = ipow(k,ex);
    double fi[sz];
    for(s=-tp;s<=tp;s++)
    {
      fi[t] = pac[ctr + s*corr].p == NULL ? pac[ctr + s*corr].fv : pac[ctr + s*corr].p->fval;
      t = t+1;
    }
    
    double fprev = pred(op, fi);
  
    double detalhe = pac[ctr].p->fil[1]->fval - fprev;
    //double diff = pac[ctr].p->fil[0]->fval - pac[ctr].p->fil[1]->fval;
    
    /* como indicado em amr (funcao), eventualmente precisa ser comentado */
    //serao alteradas apenas os valores dos filhos das celulas de profundidade lim
    assert((pac[ctr].p->fil[0] == NULL) == (pac[ctr].p->fil[1] == NULL));

    /* truncamento dos detalhes, se necessario */
    assert((pac[ctr].p->fil[0]->fil[0] == NULL) == (pac[ctr].p->fil[0]->fil[1] == NULL));
    assert((pac[ctr].p->fil[1]->fil[0] == NULL) == (pac[ctr].p->fil[1]->fil[1] == NULL));
    if(sim)
    { // se sim == TRUE, entra aqui 
      // o parametro de truncamento para medias celulares nao e constante
      // depende da profundidade na malha
      //fprintf(stderr,"opa\n");
      double pot = ( prof+1 - niv ); 
      double eps_l = eps*pow(2.0, pot);
      bool_t descN;
      if((pac[ctr].p->fil[0]->fil[0]==NULL) && (pac[ctr].p->fil[1]->fil[0]==NULL))
      { descN = TRUE; }
      else
      { descN = FALSE;}
      bool_t pode = TRUE;
      
      auto bool_t podemesmo(int d, int k, int niv, int prof);
      bool_t podemesmo(int d, int k, int niv, int prof)
      {
        int i;
        int lm = k == 5 ? 1 : 2;
        int corr = ipow(k,d-1);
        int fil;
        switch(d)
        {
          case 1:
            fil = i==1 ? 0 : 1;
            for(i=-lm;i<=lm;i++)
            {
              if(prof < (niv-3))
              {
                if(pac[ctr+i].p!=NULL)
                { 
                  if(pac[ctr+i].p->fil[fil]!=NULL)
                  {
                    if(pac[ctr+i].p->fil[fil]->fil[0]!=NULL)
                    { pode = FALSE; }
                  }
                }
              }
            }
            break;
          case 2:
            fil = i==1 ? 0 : 1;
            for(s=-lm;s<=lm;s++)
            {
              for(i=-lm;i<=lm;i++)
              {
                if(prof < (niv-3))
                {
                  if(pac[ctr+s*corr+i].p!=NULL)
                  { 
                    if(pac[ctr+s*corr+i].p->fil[fil]!=NULL)
                    {
                      if(pac[ctr+s*corr+i].p->fil[fil]->fil[0]!=NULL)
                      { pode = FALSE;}
                    }
                  }
                }
              }
            }
            break;
          case 3:
            fil = i==1 ? 0 : 1;
            for(t=-lm;t<=lm;t++)
            {
              for(s=-lm;s<=lm;s++)
              {
                for(i=-lm;i<=lm;i++)
                {
                  if(prof < (niv-3))
                  {
                    if(pac[ctr+s*corr+t*corr*corr+i].p!=NULL)
                    { 
                      if(pac[ctr+s*corr+t*corr*corr+i].p->fil[fil]!=NULL)
                      {
                        if(pac[ctr+s*corr+t*corr*corr+i].p->fil[fil]->fil[0]!=NULL)
                        { pode = FALSE; }
                      }
                    }
                  }
                }
              }
            }
            break;
          default: 
            fprintf(stderr,"caso nao contemplado - poda_hartenv2\n");
            break;
        }
        return pode;
      }
      bool_t podeM = podemesmo(d, k, niv, prof);
      if((fabs(detalhe) < eps_l) && (podeM && descN))
      { 
        pac[ctr].p->fval = 0.5*(pac[ctr].p->fil[0]->fval + pac[ctr].p->fil[1]->fval);
        free(pac[ctr].p->fil[0]); pac[ctr].p->fil[0] = NULL;
        free(pac[ctr].p->fil[1]); pac[ctr].p->fil[1] = NULL;
      }
    }
    return;
  }
}


void refina_harten(int d, bool_t op, VNo pac[], double xmin[], double xmax[], int prof, Preditor pred, int niv)
{
  int  k = op ? 9 : 5;
  int tm = ipow(k,d);
  int ctr = tm/2;
  
  if(prof == (niv-1)) 
  { /* Não há o que fazer: */
    //fprintf(stderr,"niv = %d  prof = %d lim = %d \n",niv,prof,lim);
    pac[ctr].p->fil[0] = NULL;
    pac[ctr].p->fil[1] = NULL;
    return;
  }
  
  /* Decide o eixo {ex} perpendicular ao corte: */
  int ex = prof % d;
  // cria uma arvore, sera usada se o filho for nulo
  assert((pac[ctr].p->fil[0] == NULL) == (pac[ctr].p->fil[1] == NULL));
  if(pac[ctr].p->fil[0] == NULL)
  { //return;
    
    int sz = op ? 5 : 3;
    int tp = sz/2, t = 0, s;
    double fi[sz];
    int corr = ipow(k,ex);
//fprintf(stderr, "corr = %d",corr);
    for(s=-tp;s<=tp;s++)
    {
      fi[t] = pac[ctr + s*corr].p == NULL ? pac[ctr + s*corr].fv : pac[ctr + s*corr].p->fval;
      t = t+1;
    }
    
    double fprev = pred(op, fi);
  
    pac[ctr].p->fil[0] = malloc(sizeof(No));
    pac[ctr].p->fil[0]->fval = 2.0*pac[ctr].p->fval - fprev;
    pac[ctr].p->fil[0]->fv[0] = 2.0*pac[ctr].p->fval - fprev;
    pac[ctr].p->fil[0]->fv[1] = 2.0*pac[ctr].p->fval - fprev;
    pac[ctr].p->fil[0]->fil[0] = NULL;
    pac[ctr].p->fil[0]->fil[1] = NULL;
    pac[ctr].p->fil[1] = malloc(sizeof(No));
    pac[ctr].p->fil[1]->fval = fprev;
    pac[ctr].p->fil[1]->fv[0] = fprev;
    pac[ctr].p->fil[1]->fv[1] = fprev;
    pac[ctr].p->fil[1]->fil[0] = NULL;
    pac[ctr].p->fil[1]->fil[1] = NULL;
    
    return;
  }
  
  int qual;
  for (qual = 0; qual <= 1; qual++)
  { 
    int tam = pow(k,d);
    VNo pacf[tam];
    double xminf[d], xmaxf[d];
    /* Constrói o pacote do filho {qual}: */
    divide_pacote_geral(d, pac, ex, pred, qual, pacf, op);
    limites_celula( d, ex, qual, xmin, xmax, xminf, xmaxf);
    refina_harten(d, op, pacf, xminf, xmaxf, prof+1, pred, niv);
  }
}

void refina(int d, bool_t op, VNo pac[], int niv, double xmin[], double xmax[], int prof, Preditor pred, bool_t haar)
{
  if(haar)
  {
    int ord = 5;
    int ctr = ipow(ord,d)/2;
    No *r = pac[ctr].p;
    refina_haar( r, prof, niv);
  }
  else
  {
    refina_harten( d, op, pac, xmin, xmax, prof, pred, niv);
  }
}


void poda(int d, bool_t op, VNo pac[], int niv, double xmin[], double xmax[], int prof, Preditor pred, double eps, bool_t sim, bool_t haar)
{
  if(haar)
  {
    int ord = 5;
    int ctr = ipow(ord,d)/2;
    No *r = pac[ctr].p;
    poda_haar( r, niv, prof, eps, sim);
  }
  else
  {
    int lim;
    for(lim = niv-2;lim >= 0;lim--)
    {
      //fprintf(stderr, "lim = %d\n", lim);
//      poda_harten( d, op, pac, niv, xmin, xmax, prof, pred, eps, lim, sim);
      poda_hartenv2( d, op, pac, niv, xmin, xmax, prof, pred, eps, lim, sim);
    }
  }
  //fprintf(stderr,"podou\n");
}