/* Last edited on 2009-11-09 17:36:50 by stolfilocal */
/***************************************************************************
 *   Copyright (C) 2009 by Douglas Castro   *
 *   douglas@ime.unicamp.br   *
 *   Last edited on 2009-07-30 16:17:49 by stolfilocal
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#include "arvore.h"

// este esta correto, a versao abaixo eh pra tentar um codigo mais compacto.
void divide_pacote_geral(int d, VNo pac[], int ex, Preditor pred, int qual, VNo fil[], bool_t op)
{
  int k = op ? 9 : 5;

  int i,j,l;
  for(l=0; l<(d==3 ? k:1); l++) // reservado para o eixo 2
  {
    for(j=0;j<(d>=2 ? k:1); j++) // reservado para o eixo 1
    {
      for(i=0;i<k;i++) // reservado para o eixo 0
      {
        int m,n,p;
        //filho direito corte ortog ex=0
        m = ( ex == 0 ? (op ? 2 : 1) + (i+qual)/2 : i);
        //filho direito corte ortog ex=1
        n = ( ex == 1 ? (op ? 2 : 1) + (j+qual)/2 : j);
        //filho direito corte ortog ex=2
        p = ( ex == 2 ? (op ? 2 : 1) + (l+qual)/2 : l);
//fprintf(stderr,"m = %d, n = %d, p = %d\n", m,n,p);
				
        No *pl;
				
        pl = pac[l*k*k + n*k + m].p;
        int sz = op ? 5 : 3;
        int tp = sz/2, t = 0, s;
        double fi[sz];
        for(s=-tp; s<=tp; s++)
        {
          fi[t] = pac[l*k*k + n*k + m + s].fv;
          t = t+1;
        }
        // montando os pacotes
        double fb = fi[sz/2]; 
        double fint = pred(op, fi);
        No *p_fil;
        double fv_fil;
        if(ex==0)
        {
          if(qual==0)
          {
            p_fil = (pl == NULL ? NULL : pl->fil[i%2]);
            fv_fil = (p_fil == NULL ? (i%2 == 0 ? 2.0*fb - fint : fint) : p_fil->fval);
//fprintf(stderr,"ex==0 q = 0, pacf[%d] = pac[%d]->[%d]\n", l*k*k+j*k+i, l*k*k + n*k + m, i%2);
          }
          else
          {
            p_fil = (pl == NULL ? NULL : pl->fil[(i+1)%2]);
            fv_fil = (p_fil == NULL ? (i%2 == 0 ? fint : 2.0*fb - fint) : p_fil->fval);
//fprintf(stderr,"ex==0 q = 1, pacf[%d] = pac[%d]->[%d]\n", l*k*k+j*k+i, l*k*k + n*k + m, (i+1)%2);
          }
        }
        else if(ex==1)
        {
          if(qual==0)
          {
            p_fil = (pl == NULL ? NULL : pl->fil[j%2]);
            fv_fil = (p_fil == NULL ? (j%2 == 0 ? 2.0*fb - fint : fint) : p_fil->fval);
//fprintf(stderr,"ex==1 q = 0, pacf[%d] = pac[%d]->[%d]\n", l*k*k+j*k+i, l*k*k + n*k + m, j%2);
          }
          else
          {
            p_fil = (pl == NULL ? NULL : pl->fil[(j+1)%2]);
            fv_fil = (p_fil == NULL ? (j%2 == 0 ? fint : 2.0*fb - fint) : p_fil->fval);
//fprintf(stderr,"ex==1 q = 1, pacf[%d] = pac[%d]->[%d]\n", l*k*k+j*k+i, l*k*k + n*k + m, (j+1)%2);
          }
        }
        else
        {
          if(qual==0)
          {
            p_fil = (pl == NULL ? NULL : pl->fil[l%2]);
            fv_fil = (p_fil == NULL ? (l%2 == 0 ? 2.0*fb - fint : fint)  : p_fil->fval);
//fprintf(stderr,"ex==2 q = 0, pacf[%d] = pac[%d]->[%d]\n", l*k*k+j*k+i, l*k*k + n*k + m, l%2);
          }
          else
          {
            p_fil = (pl == NULL ? NULL : pl->fil[(l+1)%2]);
            fv_fil = (p_fil == NULL ? (l%2 == 0 ? fint : 2.0*fb - fint) : p_fil->fval);
//fprintf(stderr,"ex==2 q = 1, pacf[%d] = pac[%d]->[%d]\n", l*k*k+j*k+i, l*k*k + n*k + m, (l+1)%2);
          }
        }
        fil[l*k*k+j*k+i].p = p_fil;
        fil[l*k*k+j*k+i].fv = fv_fil;
      }
    }
  }
}


No *completa(int d, int niv, int prof, int indice, double xmin[], double xmax[], Funcao *f, Quadratura integra)
{ 
  if (niv == 0) { /*fprintf(stderr,"prof = %d\n",prof);*/ return NULL; }
  No *e = malloc(sizeof(No));

  if (f != NULL) 
  { 
    double fint =  integra( d, xmin, xmax, f);
    double cvol = 1.0;
    int i;
    for(i=0;i<d;i++)
      cvol *= fabs(xmax[i]-xmin[i]);
    e->fval = fint/cvol;
    e->del = FALSE;
  }
  else
  { e->fval = 0.0; }
  
  /* Para ODEs: */
  e->fv[0] = e->fval;
  e->fv[1] = e->fval;
  e->fv[2] = e->fval;

  int ex = prof%d;
  int qual;
  for(qual=0;qual<2;qual++)
  {
    double xmin_fil[d], xmax_fil[d];
    limites_celula( d, ex, qual, xmin, xmax, xmin_fil, xmax_fil);
    e->fil[qual] = completa(d, niv - 1, prof + 1, 2*indice + qual, xmin_fil, xmax_fil, f, integra);
  }

  return e;
}

No *comp_sol_ana(int d, int niv, int prof, double xmin[], double xmax[], Funcao *f, QuadraturaSA integraSA, double t)
{ 
  if (niv == 0) { return NULL; }
  No *e = malloc(sizeof(No));

  if (f != NULL) 
  { 
    double fint =  integraSA( d, xmin, xmax, f, t);
    double cvol = 1.0;
    int i;
    for(i=0;i<d;i++)
      cvol *= fabs(xmax[i]-xmin[i]);
    e->fval = fint/cvol;
    e->del = FALSE;
  }
  else
  { e->fval = 0.0; }
  
  /* Para ODEs: */
  e->fv[0] = e->fval;
  e->fv[1] = e->fval;
  e->fv[2] = e->fval;

  int ex = prof%d;
  int qual;
  for(qual=0;qual<2;qual++)
  {
    double xmin_fil[d], xmax_fil[d];
    limites_celula( d, ex, qual, xmin, xmax, xmin_fil, xmax_fil);
    e->fil[qual] = comp_sol_ana(d, niv - 1, prof +1,  xmin_fil, xmax_fil, f, integraSA, t);
  }

  return e;
}

No *comp_sol_ana_burg(int d, int niv, int prof, double xmin[], double xmax[], Funcao *f, QuadraturaSA integraSA)
{ 
  if (niv == 0) { return NULL; }
  No *e = malloc(sizeof(No));

  if (f != NULL) 
  {
    double fint =  integraSA( d, xmin, xmax, f, 0.0);
    double cvol = 1.0;
    int i;
    for(i=0;i<d;i++)
      cvol *= fabs(xmax[i]-xmin[i]);
    e->fval = fint/cvol;
    e->del = FALSE;
  }
  else
  { e->fval = 0.0; }
  
  /* Para ODEs: */
  e->fv[0] = e->fval;
  e->fv[1] = e->fval;
  e->fv[2] = e->fval;

  int ex = prof%d;
  int qual;
  for(qual=0;qual<2;qual++)
  {
    double xmin_fil[d], xmax_fil[d];
    limites_celula( d, ex, qual, xmin, xmax, xmin_fil, xmax_fil);
    e->fil[qual] = comp_sol_ana_burg(d, niv - 1, prof +1,  xmin_fil, xmax_fil, f, integraSA);
  }

  return e;
}
void conta_folhas(No *r, int *cont)
{
  //if (r == NULL) { return;}
  assert((r->fil[0] == NULL) == (r->fil[1] == NULL));
  
  if(r->fil[0] == NULL)
  {
    int c = *cont;
    *cont = c + 1;
    return;
  }
  else
  {
    conta_folhas(r->fil[0], cont);
    conta_folhas(r->fil[1], cont);
  }
}

void conta_nos(No *r, int *cont)
{
  if (r == NULL) { return;}
  assert((r->fil[0] == NULL) == (r->fil[1] == NULL));
  
  int c = *cont;
  *cont = c + 1;
  conta_nos(r->fil[0], cont);
  conta_nos(r->fil[1], cont);
}

void copia_folhas(No *r, double vet[], int *cont)
{
  assert((r->fil[0] == NULL) == (r->fil[1] == NULL));
  if(r == NULL){return;}
  if(r->fil[0] == NULL)
  {
    int c = *cont;
    double db = r->fval;
    vet[c] = db;
    *cont = c + 1;
    return;
  }
  int qual;
  for(qual=0;qual<2;qual++)
    copia_folhas(r->fil[qual], vet, cont);
}

/* necessaria */
void atualiza_arvore(No *r)
{
  /* Paranoia */
  assert((r->fil[0] == NULL) == (r->fil[1] == NULL));
  if(r->fil[0]==NULL) { return; }

  int qual;
  for (qual = 0; qual <= 1; qual++)
  { 
    atualiza_arvore(r->fil[qual]);
  }
  
  assert((r->fil[0] != NULL) == (r->fil[1] != NULL));
  if(r->fil[0] != NULL)
  {
    /* fprintf(stderr,"atualiza_arvore - ent\n"); */
    double ce = r->fil[0]->fval;
    double cd = r->fil[1]->fval;
    r->fval = 0.5*(ce + cd);
    r->fv[0] = 0.5*(r->fil[0]->fv[0] + r->fil[1]->fv[0]);
    r->fv[1] = r->fval;
    if(r->fv[1] != r->fval)
    {
      fprintf(stderr,"fval = %f\tfv1 = %f\n",r->fval,r->fv[1]);
      fprintf(stderr,"fval = %f\tfv1 = %f\n",r->fval,r->fv[1]);
    } 
    assert(r->fv[1] == r->fval);
  }
}


void atualiza_nos_internos(No *r)
{
  /* Paranoia */
  assert((r->fil[0] == NULL) == (r->fil[1] == NULL));
  if(r->fil[0]==NULL) { return; }

  int qual;
  for (qual = 0; qual <= 1; qual++)
  { 
    atualiza_nos_internos(r->fil[qual]);
  }
  
  assert((r->fil[0] != NULL) == (r->fil[1] != NULL));
  if(r->fil[0] != NULL)
  {
    /* fprintf(stderr,"atualiza_arvore - ent\n"); */
    double ce = r->fil[0]->fval;
    double cd = r->fil[1]->fval;
    r->fval = 0.5*(ce + cd);
    r->fv[0] = r->fval;
    r->fv[1] = r->fval;
    r->fv[2] = r->fval;
  }
}

void atualiza_folhas( No *r)
{
  if(r == NULL)
  { return; }
 
  int qual;
  for(qual = 0; qual <= 1; qual++)
  { 
    atualiza_folhas( r->fil[qual] );
  }
  assert((r->fil[0]==NULL)==(r->fil[1]==NULL));
  
  if(r->fil[0]==NULL)
  {
    double tn = r->fval;
    double tnmaismeio = r->fv[0];
    r->fv[1] = tnmaismeio;
    r->fval = tnmaismeio;
    assert(r->fv[1] == r->fval );
    r->fv[0] = tn;
  }
}

void atualiza_folhas_v2( No *r, int t)
{
  if(r == NULL){ return; }
 
  int qual;
  for(qual = 0; qual <= 1; qual++)
  { 
    atualiza_folhas_v2( r->fil[qual] , t);
  }
  assert((r->fil[0]==NULL)==(r->fil[1]==NULL));
  
  if(r->fil[0]==NULL)
  {
    double tn = r->fval;
    double tnmais = r->fv[t];
    r->fval = tnmais;
    r->fv[t] = tn;
    if(t==0){r->fv[1] = tnmais; r->fv[2] = tnmais;}
    else if(t==1){r->fv[2] = tnmais; }
    else if(t==2){r->fv[0] = tnmais; r->fv[1] = tnmais; r->fv[2] = tnmais;}
    else{fprintf(stderr,"breo\n");}
  }
}

void atualiza_folhas_rk3p1( No *r)
{
  if(r == NULL)
  { return; }
 
  int qual;
  for(qual = 0; qual <= 1; qual++)
  { 
    atualiza_folhas_rk3p1( r->fil[qual] );
  }
  assert((r->fil[0]==NULL)==(r->fil[1]==NULL));
  
  if(r->fil[0]==NULL)
  {
    double tn = r->fval;
    double tn13 = r->fv[0];
    r->fval = tn13;
    r->fv[0] = tn;
    r->fv[1] = tn13;
    r->fv[2] = tn13;
  }
}

void atualiza_folhas_rk3p2( No *r)
{
  if(r == NULL)
  { return; }
 
  int qual;
  for(qual = 0; qual <= 1; qual++)
  { 
    atualiza_folhas_rk3p2( r->fil[qual] );
  }
  assert((r->fil[0]==NULL)==(r->fil[1]==NULL));
  
  if(r->fil[0]==NULL)
  {
    double tn13 = r->fval;
    double tn23 = r->fv[1];

    r->fval = tn23;
    r->fv[1] = tn13;
    r->fv[2] = tn23;
  }
}

void atualiza_folhas_rk3p3( No *r)
{
  if(r == NULL)
  { return; }
 
  int qual;
  for(qual = 0; qual <= 1; qual++)
  { 
    atualiza_folhas_rk3p3( r->fil[qual] );
  }
  assert((r->fil[0]==NULL)==(r->fil[1]==NULL));
  
  if(r->fil[0]==NULL)
  {
    r->fval = r->fv[2];
    r->fv[0] = r->fv[2];
    r->fv[1] = r->fv[2];
  }
}

void atualiza_arvore_rk3(No *r)
{
  /* Paranoia */
  assert((r->fil[0] == NULL) == (r->fil[1] == NULL));
  if(r->fil[0]==NULL) { return; }

  int qual;
  for (qual = 0; qual <= 1; qual++)
  { 
    atualiza_arvore_rk3(r->fil[qual]);
  }
  
  assert((r->fil[0] != NULL) == (r->fil[1] != NULL));
  if(r->fil[0] != NULL)
  {
    /* fprintf(stderr,"atualiza_arvore - ent\n"); */
    double ce = r->fil[0]->fval;
    double cd = r->fil[1]->fval;
    r->fval = 0.5*(ce + cd);
    r->fv[0] = r->fval;
    r->fv[1] = r->fval;
    r->fv[2] = r->fval;

  }
}

bool_t poda_arvore(int d, bool_t op, int niv, VNo pac[], int prof, Preditor pred, double fpoda, double eps)
{
  int  ord=5; // posteriormente k pode se tornar argumento da funcao
	if(op){ord=9;}
  int ctr = pow(ord,d)/2;
  
  if (pac[ctr].p == NULL) 
  { /* Não há o que podar: */
    return FALSE;
  }
  
  /* fprintf(stderr, "%*sf(int) = %9.6f f(arv)= %9.6f\n", 2*profund, "", pac[ctr].fv, pac[ctr].p->fval); */
  
  /* Precisa podar as sub-árvores, depois ver se a raiz é supérflua. */

  /* Decide o eixo {ex} perpendicular ao corte: */
  int ex = prof % d;
  
  int qual;
  bool_t filho_superfluo[2]; /* TRUE para filhos que podem ser podados. */
  for (qual = 0; qual <= 1; qual++)
  { 
    int tam = pow(ord,d);
    VNo pac_fil[tam];
  
    /* Constrói o pacote do filho {qual}: */
    divide_pacote_geral(d, pac, ex, pred, qual, pac_fil,op);
    
    /* Calcula o valor que {pac_fil[ctr].fv} teria se esse nó fosse NULL: */
    int sz = op ? 5 : 3;
		int tp = sz/2, t = 0, s;
    int corr = ipow(ord,ex);
    double fi[sz];
    for(s=-tp;s<=tp;s++)
    {
      fi[t] = pac[ctr + s*corr].fv;
      t = t+1;
    }
    double fb = fi[sz/2]; 
    double fint = pred(op, fi);
    double fpoda_fil = qual==1 ? fint : 2.0*fb - fint; 

    filho_superfluo[qual] = poda_arvore(d, op, niv, pac_fil, prof+1, pred, fpoda_fil, eps);
  }
  /* Elimina os filhos, se sao nao-nulos mas ambos superfluos: */
  if ((filho_superfluo[0] && filho_superfluo[1]) && prof>3)
  { 
    /*fprintf(stderr, "%*s\n", 2*profund, "+");*/
    /* 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;
  }
  
  /* Paranoia: */
  assert((pac[ctr].p->fil[0] == NULL) == (pac[ctr].p->fil[1] == NULL));
        
  /* Se {pac[ctr].p} tem filhos, nao e superfluo: */
  if (pac[ctr].p->fil[0] != NULL) { return FALSE; }

  /* Vale a pena manter a raiz? */
  double detalhe = fpoda - pac[ctr].fv;
  /* fprintf(stderr, "%*sdetalhe = %9.6f\n", 2*profund, "", detalhe); */
  double l = prof;
  double tol = eps*pow(2.0,(l+1-niv)); //double tol = pow(2.0,d*(l-(niv+1)))*eps;
  if ((fabs(detalhe) < tol) && (prof>=4))
  { /* O nó central é supérfluo: */
    /* fprintf(stderr, "%*ssuperfluo!\n", 2*profund, ""); */
    return TRUE;
  }
  else
  { /* O nó central precisa ser mantido: */
    /* fprintf(stderr, "%*simportante!\n", 2*profund, ""); */
    return FALSE;
  }
}

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

void copia_arvore(No *r, No *s)
{
  if(r==NULL) { return;}
  /* Paranoia */
  assert((s->fil[0]==NULL)==(s->fil[1]==NULL));
  assert((r->fil[0]==NULL)==(r->fil[1]==NULL));
  if((s->fil[0]==NULL) && (r->fil[0]!=NULL))
  {
    s->fval = r->fval;
    s->fv[0] = r->fv[0];
    s->fv[1] = r->fv[1];
    s->fil[0] = malloc(sizeof(No));
    s->fil[1] = malloc(sizeof(No));
    s->fil[0]->fil[0] = NULL;
    s->fil[0]->fil[1] = NULL;
    s->fil[1]->fil[0] = NULL;
    s->fil[1]->fil[1] = NULL;
  }
  else
  {
    s->fval = r->fval;
    s->fv[0] = r->fv[0];
    s->fv[1] = r->fv[1];
  }
    
  int qual;
  for(qual=0;qual<2;qual++)
    copia_arvore(r->fil[qual],s->fil[qual]);
  
}

void encontra_vizinhos_inferiores_esq_refinada(int d, No *r, int prof, int dir, double xmin[], double xmax[], double xmm[], double *somafluxo)
{
  if(r == NULL) { return; }
  double delta3 = (xmax[dir] - xmin[dir])/3.0;
  if((xmax[dir] + delta3) <= xmm[dir])
  { /* isto eh, o centro do pacote nao contribui com o fluxo na celula de intteresse, com coord inferior xmm */
    return;
  }
  
  assert((r->fil[0] == NULL)==(r->fil[1] == NULL));

  if((r->fil[0] == NULL)&&(xmm[dir] <= (xmax[dir]+delta3)))// qto a desigualdade, estou pedindo demais?
  {
    /* */
    double tmp = *somafluxo;
    tmp += r->flxsup[dir];
    *somafluxo = tmp;
    return;
  }

  int qual;
  for(qual=0;qual<2;qual++)
  {
    int ex = prof % d;
    double xminf[d], xmaxf[d];
    limites_celula( d, ex, qual, xmin, xmax, xminf, xmaxf);
    encontra_vizinhos_inferiores_esq_refinada( d, r->fil[qual], prof+1, dir, xminf, xmaxf, xmm, somafluxo);
  }
}

void encontra_vizinhos_superiores_dir_refinada(int d, No *r, int prof, int dir, double xmin[], double xmax[], double xMM[], double *somafluxo)
{
  if(r == NULL) { return; }
  double delta3 = (xmax[dir] - xmin[dir])/3.0;
  double lv = xmin[dir] - delta3;
  if( lv >= xMM[dir])
  { /* isto eh, o centro do pacote nao contribui com o fluxo na celula de interesse, com coord superior xMM */
    return;
  }
  
  if((r->fil[0] == NULL)/*&&(lv < xMM[dir])*/) // qto a desigualdade
  {
    /* */
    double tmp = *somafluxo;
    tmp += r->flxinf[dir];
    *somafluxo = tmp;
    return;
  }

  int qual;
  for(qual=0;qual<2;qual++)
  {
    int ex = prof % d;
    double xminf[d], xmaxf[d];
    limites_celula( d, ex, qual, xmin, xmax, xminf, xmaxf);
    encontra_vizinhos_superiores_dir_refinada( d, r->fil[qual], prof+1, dir, xminf, xmaxf, xMM, somafluxo);
  }
}

void encontra_media_inferior_refinada(int d, No *r, int prof, int dir, double xmin[], double xmax[], double xmm[], double *media)
{
  if(r == NULL) { return; }
  double delta2 = (xmax[dir] - xmin[dir])/2.0;
  if((xmax[dir] + delta2) <= xmm[dir])
  { /* isto eh, o centro do pacote nao contribui com o fluxo na celula de intteresse, com coord inferior xmm */
    return;
  }
  
  if((r->fil[0] == NULL)/*&&(xmm[dir] <= (xmax[dir]+delta2))*/)// qto a desigualdade, estou pedindo demais?
  {
    /* */
    *media = r->fval;
    return;
  }

  int qual;
  for(qual=0;qual<2;qual++)
  {
    int ex = prof % d;
    double xminf[d], xmaxf[d];
    limites_celula( d, ex, qual, xmin, xmax, xminf, xmaxf);
    encontra_media_inferior_refinada( d, r->fil[qual], prof+1, dir, xminf, xmaxf, xmm, media);
  }
}

void encontra_media_superior_refinada(int d, No *r, int prof, int dir, double xmin[], double xmax[], double xMM[], double *media)
{
  if(r == NULL) { return; }
  double delta2 = (xmax[dir] - xmin[dir])/2.0;
  double lv = xmin[dir] - delta2;
  if( lv >= xMM[dir])
  { /* isto eh, o centro do pacote nao contribui com o fluxo na celula de interesse, com coord superior xMM */
    return;
  }
  
  if((r->fil[0] == NULL)/*&&(lv < xMM[dir])*/) // qto a desigualdade
  {
    /* */
    *media = r->fval;
    return;
  }

  int qual;
  for(qual=0;qual<2;qual++)
  {
    int ex = prof % d;
    double xminf[d], xmaxf[d];
    limites_celula( d, ex, qual, xmin, xmax, xminf, xmaxf);
    encontra_media_superior_refinada( d, r->fil[qual], prof+1, dir, xminf, xmaxf, xMM, media);
  }
}

void divide_pacote_menor(int d, VFl pac[], int ex, int qual, VFl fil[])
{
  assert((d >= 1) && (d <= 3));
	
  auto void escolhe_filhos(int d, int ex, VFl pac[], int qual, VFl fil[]);
  void escolhe_filhos(int d,  int ex, VFl pac[], int qual, VFl fil[])
  {
		//controle de natalidade!
    int i, j, m, t, s;

    /* nesse tipo de pacote consideramos apenas os vizinhos inferiores mais proximos */
    switch(d){
      case 1:
        for(i=0; i<=d; i++)
        { 
          j = qual == 0 ? i : qual;
          m = qual == 0 ? 1-i : i;
          fil[i].p = pac[j].p == NULL ? NULL : pac[j].p->fil[m];
          fil[i].fvs[0] = fil[i].p == NULL ? pac[j].fvs[0] : fil[i].p->flxsup[0];
          fil[i].fvi[0] = fil[i].p == NULL ? pac[j].fvi[0] : fil[i].p->flxinf[0];
        }
        break;
      case 2:
        /* aqui as estradas representam os seguintes nos
           | 0 | 1 |
           | x | 2 |
         isto eh, so considero um vizinho inferior em cada direcao */ 
        for(i=0;i<=d;i++)
        {
          m = ex == 0 ? i/2 + qual : (i+qual)/2;
          j = qual == 0 ? i : m;
          t = ex == 0 ? 1 + (ex-qual) : 1 - (ex - qual);
          s = ex == qual ? -1 : +1;
          fil[i].p = pac[j].p == NULL ? NULL : pac[j].p->fil[t+s*(i+1-ex)/2];
          if(ex==0)
          {
            fil[i].fvs[0] = fil[i].p == NULL ? pac[j].fvs[0] : fil[i].p->flxsup[0];
            fil[i].fvs[1] = fil[i].p == NULL ? pac[j].fvs[1]/*/2.0*/ : fil[i].p->flxsup[1];
            fil[i].fvi[0] = fil[i].p == NULL ? pac[j].fvi[0] : fil[i].p->flxinf[0];
            fil[i].fvi[1] = fil[i].p == NULL ? pac[j].fvi[1]/*/2.0*/ : fil[i].p->flxinf[1];
          }
          else
          {
            fil[i].fvs[0] = fil[i].p == NULL ? pac[j].fvs[0]/*/2.0*/ : fil[i].p->flxsup[0];
            fil[i].fvs[1] = fil[i].p == NULL ? pac[j].fvs[1] : fil[i].p->flxsup[1];
            fil[i].fvi[0] = fil[i].p == NULL ? pac[j].fvi[0]/*/2.0*/ : fil[i].p->flxinf[0];
            fil[i].fvi[1] = fil[i].p == NULL ? pac[j].fvi[1] : fil[i].p->flxinf[1];
          }
        }
        break;
      case 3:
        /* aqui as estradas representam os seguintes nos
        | 0 | 1 |     | x | 3 |   (celula 3 abaixo da cel 1)
        | x | 2 |     | x | x |   ( x indica quem nao eh usado, nem passado )
        isto eh, so considero um vizinho inferior em cada direcao */ 
        for(i=0;i<=d;i++)
        {
          if(i==3)
          {
            j = /*ex == 2 && qual == 1 ? 1 :*/ i;
            t = qual;
            if( (ex == 2) && (qual == 0)){ t = 1;}
            if( (ex == 2) && (qual == 1)){ t = 0; j = 1;}
            fil[i].p = pac[j].p == NULL ? NULL : pac[j].p->fil[t];
          }
          else
          {
            if(ex==2)
            {
              j = i;
              fil[i].p = pac[j].p == NULL ? NULL : pac[j].p->fil[qual];
            }
            else
            {
              m = ex == 0 ? i/2 + qual : (i+qual)/2;
              j = qual == 0 ? i : m;
              t = ex == 0 ? 1 + (ex-qual) : 1 - (ex - qual);
              s = ex == qual ? -1 : +1;
              fil[i].p = pac[j].p == NULL ? NULL : pac[j].p->fil[t+s*(i+1-ex)/2];
            }
          }
          
          if(ex==0)
          {
            fil[i].fvs[0] = fil[i].p == NULL ? pac[j].fvs[0] : fil[i].p->flxsup[0];
            fil[i].fvs[1] = fil[i].p == NULL ? pac[j].fvs[1]/2.0 : fil[i].p->flxsup[1];
            fil[i].fvs[1] = fil[i].p == NULL ? pac[j].fvs[2]/2.0 : fil[i].p->flxsup[2];
          }
          else if(ex==1)
          {
            fil[i].fvs[0] = fil[i].p == NULL ? pac[j].fvs[0]/2.0 : fil[i].p->flxsup[0];
            fil[i].fvs[1] = fil[i].p == NULL ? pac[j].fvs[1] : fil[i].p->flxsup[1];
            fil[i].fvs[2] = fil[i].p == NULL ? pac[j].fvs[2]/2.0 : fil[i].p->flxsup[2];
          }
          else
          {
            fil[i].fvs[0] = fil[i].p == NULL ? pac[j].fvs[0]/2.0 : fil[i].p->flxsup[0];
            fil[i].fvs[1] = fil[i].p == NULL ? pac[j].fvs[1]/2.0 : fil[i].p->flxsup[1];
            fil[i].fvs[2] = fil[i].p == NULL ? pac[j].fvs[2] : fil[i].p->flxsup[2];
          }
        }
        break;
      default:
        fprintf(stderr," divide_pacote_menor -arvore.c- nao implementado para dimensao superior a 3");
        break;
    }
  }//fim funcao

  escolhe_filhos(d, ex, pac, qual, fil);
}

void limites_celula(int d, int ex, int qual, double xmin[], double xmax[], double xmin_fil[], double xmax_fil[])
{
  double xmed = (xmin[ex] + xmax[ex])/2.0;

  int i;
  // vertices da celula filho
  for(i=0;i<d;i++)
  {
    if(qual == 0)
    {
      xmin_fil[i] = xmin[i];
      xmax_fil[i] = i == ex ? xmed : xmax[i];
    }
    else
    {
      xmin_fil[i] = i == ex ? xmed : xmin[i];
      xmax_fil[i] = xmax[i];
    }
  }
}


bool_t indica_poda_arvore(int d, bool_t op, int niv, VNo pac[], int prof, Preditor pred, double fpoda, double eps)
{
  int  ord=5; // posteriormente k pode se tornar argumento da funcao
	if(op){ord=9;}
  int ctr = pow(ord,d)/2;
  
  if (pac[ctr].p == NULL) 
  { /* Não há o que podar: */
    return FALSE;
  }
  
  /* fprintf(stderr, "%*sf(int) = %9.6f f(arv)= %9.6f\n", 2*profund, "", pac[ctr].fv, pac[ctr].p->fval); */
  
  /* Precisa podar as sub-árvores, depois ver se a raiz é supérflua. */

  /* Decide o eixo {ex} perpendicular ao corte: */
  int ex = prof % d;
  
  int qual;
  //bool_t filho_superfluo[2]; /* TRUE para filhos que podem ser podados. */
  for (qual = 0; qual <= 1; qual++)
  { 
    int tam = pow(ord,d);
    VNo pac_fil[tam];
  
    /* Constrói o pacote do filho {qual}: */
    divide_pacote_geral(d, pac, ex, pred, qual, pac_fil,op);
    
    /* Calcula o valor que {pac_fil[ctr].fv} teria se esse nó fosse NULL: */
    int sz = op ? 5 : 3;
		int tp = sz/2, t = 0, s;
    int corr = ipow(ord,ex);
    double fi[sz];
    for(s=-tp;s<=tp;s++)
    {
      fi[t] = pac[ctr + s*corr].fv;
      t = t+1;
    }
    double fb = fi[sz/2]; 
    double fint = pred(op, fi);
    double fpoda_fil = qual==1 ? fint : 2.0*fb - fint; 

    pac_fil[ctr].p->fil[qual]->del = indica_poda_arvore(d, op, niv, pac_fil, prof+1, pred, fpoda_fil, eps);
  }
  /* Elimina os filhos, se sao nao-nulos mas ambos superfluos: */
//  if ((filho_superfluo[0] && filho_superfluo[1]) && prof>3)
//  { 
    /*fprintf(stderr, "%*s\n", 2*profund, "+");*/
    /* 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;
//  }
  
  /* Paranoia: */
  assert((pac[ctr].p->fil[0] == NULL) == (pac[ctr].p->fil[1] == NULL));
        
  /* Se {pac[ctr].p} tem filhos, nao e superfluo: */
  if (pac[ctr].p->fil[0] != NULL) { return FALSE; }

  /* Vale a pena manter a raiz? */
  double detalhe = fpoda - pac[ctr].fv;
  /* fprintf(stderr, "%*sdetalhe = %9.6f\n", 2*profund, "", detalhe); */
  double l = prof;
  double tol = eps*pow(2.0,(l+1-niv)); //double tol = pow(2.0,d*(l-(niv+1)))*eps;
  if ((fabs(detalhe) < tol) && (prof>=4))
  { /* O nó central é supérfluo: */
    /* fprintf(stderr, "%*ssuperfluo!\n", 2*profund, ""); */
    return TRUE;
  }
  else
  { /* O nó central precisa ser mantido: */
    /* fprintf(stderr, "%*simportante!\n", 2*profund, ""); */
    return FALSE;
  }
}

void threshold(int d, bool_t op, int niv, VNo pac[], int prof, Preditor pred, double fpoda, double eps)
{
  int  ord=5; // posteriormente k pode se tornar argumento da funcao
	if(op){ord=9;}
  int ctr = pow(ord,d)/2;
  
  if (pac[ctr].p == NULL) 
  { /* Não há o que podar: */
    return;
  }
  
  /* fprintf(stderr, "%*sf(int) = %9.6f f(arv)= %9.6f\n", 2*profund, "", pac[ctr].fv, pac[ctr].p->fval); */
  
  /* Precisa podar as sub-árvores, depois ver se a raiz é supérflua. */

  /* Decide o eixo {ex} perpendicular ao corte: */
  int ex = prof % d;
  
  int qual;
  //bool_t filho_superfluo[2]; /* TRUE para filhos que podem ser podados. */
  for (qual = 0; qual <= 1; qual++)
  { 
    int tam = pow(ord,d);
    VNo pac_fil[tam];
  
    /* Constrói o pacote do filho {qual}: */
    divide_pacote_geral(d, pac, ex, pred, qual, pac_fil,op);
    
    /* Calcula o valor que {pac_fil[ctr].fv} teria se esse nó fosse NULL: */
    int sz = op ? 5 : 3;
		int tp = sz/2, t = 0, s;
    int corr = ipow(ord,ex);
    double fi[sz];
    for(s=-tp;s<=tp;s++)
    {
      fi[t] = pac[ctr + s*corr].fv;
      t = t+1;
    }
    double fb = fi[sz/2]; 
    double fint = pred(op, fi);
    double fpoda_fil = qual==1 ? fint : 2.0*fb - fint; 

    threshold(d, op, niv, pac_fil, prof+1, pred, fpoda_fil, eps);
  }
  /* Elimina os filhos, se sao nao-nulos mas ambos superfluos: */
  auto bool_t podemesmo(int d, int ex, int ord );
  bool_t podemesmo(int d, int ex, int ord)
  {
    int corr = ipow(ord,ex);
    bool_t veredito = FALSE, l1, l2, l3;
    switch(d){
      case 1:
        veredito = pac[ctr-1].p->fil[0]->del && pac[ctr-1].p->fil[1]->del && pac[ctr+1].p->fil[0]->del && pac[ctr+1].p->fil[1]->del;
        break;
      case 2:
        l2 = pac[ctr-corr-1].p->fil[0]->del && pac[ctr-corr-1].p->fil[1]->del && pac[ctr-corr-1].p->fil[0]->del && pac[ctr-corr-1].p->fil[1]->del && pac[ctr-corr+1].p->fil[0]->del && pac[ctr-corr+1].p->fil[1]->del;
        l1 = pac[ctr-1].p->fil[0]->del && pac[ctr-1].p->fil[1]->del && pac[ctr+1].p->fil[0]->del && pac[ctr+1].p->fil[1]->del;
        l3 = pac[ctr+corr-1].p->fil[0]->del && pac[ctr+corr-1].p->fil[1]->del && pac[ctr+corr-1].p->fil[0]->del && pac[ctr+corr-1].p->fil[1]->del && pac[ctr+corr+1].p->fil[0]->del && pac[ctr+corr+1].p->fil[1]->del;
        veredito = l1 && l2 && l3;
        break;
      //case 3:
        
      //  break;
      default:
        fprintf(stderr,"caso nao contemplado - arvore.c - podemesmo \n");
        break;
    }
    return veredito;
  }
  bool_t pode = podemesmo( d, ex, ord );
  if ( pode && (prof>3))
  { 
    /* 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;
  }


}


/*
void calcula_divergente(int d, int prof, VNo pac[], double xmin[], double xmax[], Fluxo f, Preditor pred, bool_t op)
{
  int ord = op ? 9 : 5;

  int tam = ipow(ord,d);
  int ctr = tam/2;
  assert((pac[ctr].p->fil[0] == NULL) == (pac[ctr].p->fil[1] == NULL));
  if(pac[ctr].p->fil[0] == NULL)
  { 
    int i;
    for(i=0;i<d;i++)
    {
      int corr = ipow(ord,i);
      double h = xmax[i] - xmin[i];
      // primeira derivada de ordem 2 na direcao i 
      pac[ctr].p->divf[i] = 1.0;//(f( pac[ctr+corr].fv ) - f( pac[ctr-corr].fv ))/(2.0*h);
    }
    return; 
  }
 
  int qual;
  for(qual = 0; qual <= 1; qual++)
  { 
    int ex =  prof % d;
    double xminf[d], xmaxf[d];
    VNo pacf[tam];
    limites_celula( d, ex, qual, xmin, xmax, xminf, xmaxf);
    divide_pacote_geral( d, pac, ex, pred, qual, pacf,op);
    calcula_divergente( d, prof+1, pacf, xminf, xmaxf, f, pred, op);
  }
  
}
*/

/** nessa versao fiz uma brincadeira atribuindo o valor interpolado {fint} nao ao filho direito,
 * mas ao filho esquerdo para ver no que da, eh igual, desde que a formula da interpolacao seja alterada.
 */
///void divide_pacote_geral(int d, int tp, VNo pac[], int ex, Preditor pred, int qual, VNo fil[])
//{
//  assert((d >= 1) && (d <= 3));
//  assert(tp == 5);
//
//  int kk = 5;
//	
//  auto void escolhe_filhos(int d, int kk, int ex, VNo pac[], Preditor pred,	int qual,	VNo fil[]);
//  void escolhe_filhos(int d, int kk, int ex, VNo pac[], Preditor pred,	int qual,	VNo fil[])
//  {
//		//controle de natalidade!
//    int i,j,k, ll, kz, ky, kx;
//    if(d==1){ kz = 1; ky = 1; kx = kk; ll = 0;}
//    if(d==2){ kz =  1; ky = kk; kx = kk; ll = 0;}
//   if(d==3){ kz = kk; ky = kk; kx = kk; ll = kk*kk;}
		/*  SO para lembrar que, da pra tirar dessa funcao e colocar apenas o caso 3
    *  mas e preciso lembrar que as 3 consideracoes acima sao necessarias
    */
//    switch(d){
//      case 1:
//        j=0; //for(j=0;j<ky;j++){
//        for(i=0; i<kx; i++)
//        { 
//          int m = (i + qual)/2 + 1;
//          int n = (j + qual)/2;
//          int ce = ( ex == 0 ? j : n  );
//          int cc = ( ex == 0 ? j : n+1);
//          int cd = ( ex == 0 ? j : n+2);
//          int de = ( ex == 0 ? m-1 : i);
//          int dc = ( ex == 0 ? m : i  );
//          int dd = ( ex == 0 ? m+1 : i);
 //         No *pl = pac[cc*kk + dc].p;
//          double fa = pac[ce*kk + de].fv;
//          double fb = pac[cc*kk + dc].fv;
//          double fc = pac[cd*kk + dd].fv;
//          double fint = pred(fa, fb, fc);
//          No *p_fil;
//          double fv_fil;
//          double vfint[] ={ fint, 2.0*fb - fint}; 
//          int q = ( ex == 0 ? i : j);
//          if(q%2 == 0)
//          { 
//            int z = qual;
//            p_fil = (pl == NULL ? NULL : pl->fil[z]);
//            fv_fil = (p_fil == NULL ? vfint[z] : p_fil->fval);
//          }
//          else
//          {
//            int z = (qual == 0 ? 1 : 0); 
//            p_fil = (pl == NULL ? NULL : pl->fil[z]);
//            fv_fil = (p_fil == NULL ? vfint[z] : p_fil->fval);
//          }
//          fil[j*kk + i].p = p_fil;
//          fil[j*kk + i].fv = fv_fil;
//        }//}
//        break;
//      case 2:
//        for(j=0; j<ky; j++)
//        { 
//          for(i=0; i<kx; i++)
//          { 
//            int m = (i + qual)/2 + 1;
//            int n = (j + qual)/2;
//            int ce = ( ex == 0 ? j : n  );
//            int cc = ( ex == 0 ? j : n+1);
//            int cd = ( ex == 0 ? j : n+2);
//            int de = ( ex == 0 ? m-1 : i);
//            int dc = ( ex == 0 ? m : i  );
//            int dd = ( ex == 0 ? m+1 : i);
//            No *pl = pac[cc*kk + dc].p;
//            double fa = pac[ce*kk + de].fv;
//            double fb = pac[cc*kk + dc].fv;
//            double fc = pac[cd*kk + dd].fv;
//            double fint = pred(fa, fb, fc);
//            No *p_fil;
//            double fv_fil;
//            double vfint[] ={ fint, 2.0*fb - fint}; 
//            int q = ( ex == 0 ? i : j);
//            if(q%2 == 0)
//            { 
//              int z = qual;
//              p_fil = (pl == NULL ? NULL : pl->fil[z]);
//              fv_fil = (p_fil == NULL ? vfint[z] : p_fil->fval);
//            }
//            else
//            {
//              int z = (qual == 0 ? 1 : 0); 
//              p_fil = (pl == NULL ? NULL : pl->fil[z]);
//              fv_fil = (p_fil == NULL ? vfint[z] : p_fil->fval);
//            }
//            fil[j*kk + i].p = p_fil;
//            fil[j*kk + i].fv = fv_fil;
//          }
//        }
//        break;
//      case 3:
//        for(k=0;k<kz;k++)
//        {
//          for(j=0;j<ky;j++)
//          {
//            for(i=0;i<kx;i++)
//            {
//              int m = (i + qual)/2 + 1;
//              int n = (j + qual)/2;
//              int o = (k + qual)/2 + 1;
//              int ce = ( ex == 0 ? j : n  );
//              int cc = ( ex == 0 ? j : n+1);
//              int cd = ( ex == 0 ? j : n+2);
//              int de = ( ex == 0 ? m-1 : i);
//              int dc = ( ex == 0 ? m : i  );
//              int dd = ( ex == 0 ? m+1 : i);
//              //int ee = ( ex == 2 ? o-1 : k);
//              int ec = ( ex == 2 ? o :   k);
//              //int ed = ( ex == 2 ? o+1 : k);
//              No *pl = pac[ec*ll + cc*kk + dc].p;
//              double fa = pac[ec*ll + ce*kk + de].fv;
//              double fb = pac[ec*ll + cc*kk + dc].fv;
//              double fc = pac[ec*ll + cd*kk + dd].fv;
//              double fint = pred(fa, fb, fc);
//              No *p_fil;
//              double fv_fil;
//              double vfint[] = { fint, 2.0*fb - fint}; 
							
//              int q = ( ex == 0 ? i : j);
//              q = ( ex == 2 ? k : q);
//              if(q%2 == 0)
//              { 
//                int z = qual;
//                p_fil = (pl == NULL ? NULL : pl->fil[z]);
//                fv_fil = (p_fil == NULL ? vfint[z] : p_fil->fval);
//              }
//              else
//              {
//                int z = (qual == 0 ? 1 : 0); 
// 								z = ();//parece que nao precisa de correcao aqui
//                p_fil = (pl == NULL ? NULL : pl->fil[z]);
//                fv_fil = (p_fil == NULL ? vfint[z] : p_fil->fval);
//              }
//              fil[k*ll + j*kk + i].p = p_fil;
//              fil[k*ll + j*kk + i].fv = fv_fil;
//            }
//          }
//        }
//        break;
//      default:
//        fprintf(stderr," divide_pacote -arvore.c- nao implementado para dimensao superior a 3");
//        break;
//    }
//  }//fim funcao
	
//  escolhe_filhos(d,  kk, ex, pac, pred,	qual,	fil);
	
//}