/* Last edited on 2008-08-15 11:05:17 by stolfi */
/***************************************************************************
 *   Copyright (C) 2008 by Douglas Castro,,,   *
 *   douglas@douglas-laptop   *
 *                                                                         *
 ***************************************************************************/

#include "arvore.h"

void mostra_pacote(FILE *arq, VReg pac[], int k, double xmin[], double xmax[])
  {
    fprintf(stderr, "------------------------------------------------\n");
    int ex;
    for (ex = 0; ex < 2; ex++) 
      { fprintf(stderr, "  x[%d] = [ %14.9f %14.9f ]\n", ex, xmin[ex], xmax[ex]); }
    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");
      }
    fprintf(stderr, "------------------------------------------------\n");
  }

void altera_valor_em_pacote(VReg pac[], int k, int i0, int i1, double xmin[], double xmax[], double x[], double v)
  { int i = k*i1 + i0; 
    pac[i].fv += v;
    if (pac[i].p != NULL)
      { altera_valor_em_arvore(pac[i].p, 0, xmin, xmax, x, v); }
  }

void altera_valor_em_arvore(Reg *p, int prof, double xmin[], double xmax[], double x[], double v)
  {
    double xlo[2], xhi[2];
    int ex; 
    for (ex = 0; ex < 2; ex++)
      { xlo[ex] = xmin[ex]; 
        xhi[ex] = xmax[ex];
        if (! ((xlo[ex] <= x[ex]) && (x[ex] <= xhi[ex])))
          { fprintf(stderr, "%24.16e %24.16e %24.16e\n", xlo[ex], xhi[ex], x[ex]); }
        assert((xlo[ex] <= x[ex]) && (x[ex] <= xhi[ex]));
      }
    do { 
      p->fval += v;
      /* Decide em que direção descer: */
      ex = prof % 2;
      double xmd = (xlo[ex] + xhi[ex])/2;
      if (x[ex] < xmd)
        { p = p->esq; 
          xhi[ex] = xmd;
        }
      else
        { p = p->dir;
          xlo[ex] = xmd;
        }
      v *= 2;
      prof++;
    } while (p != NULL);
  }



double avalia_arvore(VReg pac[], double xmin[], double xmax[], int niv, int profund, Interp_quad interpola_quad, double x[], int debug)
{
  int k=5; // Todo pacote tem {k*k} elementos (depende da ordem de interpolação). 
  //int ctr = (k*k)/2;
  if (debug)
    { mostra_pacote(stderr, pac, k, xmin, xmax); }
  
  if ((niv <= 0) && pacote_trivial(pac,k)) 
    { return interpola_ponto_c(pac, xmin, xmax, x); }
  /* Pacote precisa ser dividido: */
  
  /* Decide o eixo {ex} perpendicular ao corte: */
  int ex = profund % 2;
  
  /* Decide qual filho queremos (0 = esq/inf, 1 = dir/sup): */
  int qual = (x[ex] > (xmin[ex] + xmax[ex])/2);
  
  /* Constrói o pacote filho: */
  VReg fil[k*k]; /* Pacote filho relevante. */
  double xmin_fil[2], xmax_fil[2]; /* Coordenadas da célula central do filho. */
  divide_pacote(pac, xmin, xmax, ex, interpola_quad, qual, fil, xmin_fil, xmax_fil);
  
  /* Recursão: */
  double res = avalia_arvore(fil, xmin_fil, xmax_fil, niv-1, profund+1, interpola_quad, x, debug);
  
  return res;
  
}

int pacote_trivial(VReg pac[], int k)
{
  int i,j = 1;
  for(i=0;i<k*k;i++)
    {
      if(pac[i].p != NULL)
        {
          j = 0;
        }
    }
  return j;
}

/** nessa versao fiz uma brincadeira atribuindo o valor interpolado {fint} nao ao filho direito,
 * mas ao filho esquerdo para ver no que da
 */
void divide_pacote(VReg pac[],double xmin[], double xmax[], int ex,	Interp_quad interpola_quad,	int qual,	VReg fil[], double xmin_fil[], double xmax_fil[])
{
	int k=5;
        // int ctr = (k*k)/2;

	double xmed = (xmin[ex] + xmax[ex])/2.0;
	xmin_fil[ex] = (qual == 0 ? xmin[ex] : xmed);
	xmax_fil[ex] = (qual == 0 ? xmed : xmax[ex]);
    
	xmin_fil[1-ex] = xmin[1-ex];
	xmax_fil[1-ex] = xmax[1-ex];
	
	int i,j;
	if(ex == 0)
	{ // montando gabarito
		for(i=0; i<k; i++)
		{ 
			if (qual == 0)
			{ /* Filho esquerdo: */
				for(j=0; j<k; j++)
				{ 
					int m = j/2 + 1;
					Reg *pl = pac[i*k+m].p;
					double fa = pac[i*k+m-1].fv;
					double fb = pac[i*k+m].fv;
					double fc = pac[i*k+m+1].fv;
					double fint = interpola_quad(fa, fb, fc);

					Reg *p_fil;
					double fv_fil;
					if(j%2 == 0)
					{ 
						p_fil = (pl == NULL ? NULL : pl->esq);
						fv_fil = (p_fil == NULL ? fint : p_fil->fval);
					}
					else
					{
						p_fil = (pl == NULL ? NULL : pl->dir);
						fv_fil = (p_fil == NULL ? 2.0*fb - fint : p_fil->fval);
					}
					fil[i*k+j].p = p_fil;
					fil[i*k+j].fv = fv_fil;
				}
			}
			else
			{ /* filho direito: */
				for(j=0; j<k; j++)
				{
					int m = (j+1)/2 + 1;
					Reg *pl = pac[i*k+m].p;
					double fa = pac[i*k+m-1].fv;
					double fb = pac[i*k+m].fv;
					double fc = pac[i*k+m+1].fv;
					double fint = interpola_quad(fa, fb, fc);
					
					Reg *p_fil;
					double fv_fil;
					if(j%2==0)
					{
						p_fil = (pl == NULL ? NULL : pl->dir);
						fv_fil = (p_fil == NULL ? 2.0*fb - fint : p_fil->fval);
					}
					else
					{
						p_fil = (pl == NULL ? NULL : pl->esq);
						fv_fil = (p_fil == NULL ? fint : p_fil->fval);
					}
					fil[i*k+j].p = p_fil;
					fil[i*k+j].fv = fv_fil;
				}
			}
		}
	}
	else
	{
      // montando gabarito
		int i,j;
		for(i=0;i<k;i++)
		{
			if(qual==1)
			{ // filho direito
				for(j=0;j<k;j++)
				{
					int z = (i+1)/2;
					Reg *pl = pac[(z+1)*k+j].p;
					double fa = pac[z*k+j].fv;
					double fb = pac[(z+1)*k+j].fv;
					double fc = pac[(z+2)*k+j].fv;
					double fint = interpola_quad(fa, fb, fc);
					Reg *p_fil;
					double fv_fil;
					if(i%2==0)
					{ 
						p_fil = (pl == NULL ? NULL : pl->dir);
						fv_fil = (p_fil == NULL ? 2.0*fb - fint : p_fil->fval);
					}
					else
					{
						p_fil = (pl == NULL ? NULL : pl->esq);
						fv_fil = (p_fil == NULL ? fint : p_fil->fval);
					}
					fil[i*k+j].p = p_fil;
					fil[i*k+j].fv = fv_fil;
				}
			}
			else
			{//filho esquerdo
				for(j=0;j<k;j++)
				{
					int z = i/2;
					Reg *pl = pac[(z+1)*k+j].p;
					double fa = pac[z*k+j].fv;
					double fb = pac[(z+1)*k+j].fv;
					double fc = pac[(z+2)*k+j].fv;
					double fint = interpola_quad(fa, fb, fc);
					Reg *p_fil;
					double fv_fil;
					if(i%2==0)
					{ 
						p_fil = (pl == NULL ? NULL : pl->esq);
						fv_fil = (p_fil == NULL ? fint : p_fil->fval);
					}
					else
					{
						p_fil = (pl == NULL ? NULL : pl->dir);
						fv_fil = (p_fil == NULL ? 2.0*fb - fint : p_fil->fval);
					}
					fil[i*k+j].p = p_fil;
					fil[i*k+j].fv = fv_fil;
				}
			}
		}
	}
        
}

/** Nessa versao de divide pacote, o valor interpolado {fint} esta atribuido ao filho direito,
 * como proposto no principio
 */

// void divide_pacote(VReg pac[],double xmin[], double xmax[], int ex,	Interp_quad interpola_quad,	int qual,	VReg fil[], double xmin_fil[], double xmax_fil[])
// {
//   int k=5;
// 
//   double xmed = (xmin[ex] + xmax[ex])/2.0;
//   xmin_fil[ex] = (qual == 0 ? xmin[ex] : xmed);
//   xmax_fil[ex] = (qual == 0 ? xmed : xmax[ex]);
//     
//   xmin_fil[1-ex] = xmin[1-ex];
//   xmax_fil[1-ex] = xmax[1-ex];
// 	
//   int i,j,ind=1;
//   if(ex == 0)
//     { // montando gabarito
//       for(i=0; i<k; i++)
//         { 
//           if (qual == 0)
//             { /* Filho esquerdo: */
//               for(j=0; j<k; j++)
//                 { 
//                   int m = j/2 + 1;
//                   Reg *pl = pac[i*k+m].p;
//                   double fa = pac[i*k+m-1].fv;
//                   double fb = pac[i*k+m].fv;
//                   double fc = pac[i*k+m+1].fv;
//                   double fint = interpola_quad(fa, fb, fc);
// 
//                   Reg *p_fil;
//                   double fv_fil;
//                   if(j%2 == 0)
//                     { 
//                       p_fil = (pl == NULL ? NULL : pl->esq);
//                       fv_fil = ((pl == NULL || pl->esq == NULL) ? 2.0*fb - fint : pl->esq->fval);
//                     }
//                   else
//                     {
//                       p_fil = (pl == NULL ? NULL : pl->dir);
//                       fv_fil = ((pl == NULL || pl->dir == NULL) ? fint : pl->dir->fval);
//                     }
//                   fil[i*k+j].p = p_fil;
//                   fil[i*k+j].fv = fv_fil;
//                 }
//             }
//           else
//             { /* filho direito: */
//               for(j=0; j<k; j++)
//                 {
//                   int m = (j+1)/2 + 1;
//                   Reg *pl = pac[i*k+m].p;
//                   double fa = pac[i*k+m-1].fv;
//                   double fb = pac[i*k+m].fv;
//                   double fc = pac[i*k+m+1].fv;
//                   double fint = interpola_quad(fa, fb, fc);
// 					
//                   Reg *p_fil;
//                   double fv_fil;
//                   if(j%2==0)
//                     {
//                       p_fil = (pl == NULL ? NULL : pl->dir);
//                       fv_fil = ((pl == NULL || pl->dir == NULL) ? fint : pl->dir->fval);
//                     }
//                   else
//                     {
//                       p_fil = (pl == NULL ? NULL : pl->esq);
//                       fv_fil = ((pl == NULL || pl->esq == NULL) ? 2.0*fb - fint : pl->esq->fval);
//                     }
//                   fil[i*k+j].p = p_fil;
//                   fil[i*k+j].fv = fv_fil;
//                 }
//             }
//         }
//     }
//   else
//     {
//       // montando gabarito
//       int i,j;
//       for(i=0;i<k;i++)
//         {
//           if(qual==1)
//             { // filho direito
//               for(j=0;j<k;j++)
//                 {
//                   int z = (i+1)/2;
//                   Reg *pl = pac[(z+1)*k+j].p;
//                   double fa = pac[z*k+j].fv;
//                   double fb = pac[(z+1)*k+j].fv;
//                   double fc = pac[(z+2)*k+j].fv;
//                   double fint = interpola_quad(fa, fb, fc);
//                   Reg *p_fil;
//                   double fv_fil;
//                   if(i%2==0)
//                     { 
//                       p_fil = (pl == NULL ? NULL : pl->dir);
//                       fv_fil = ((pl == NULL || pl->dir == NULL) ? fint : pl->dir->fval);
//                     }
//                   else
//                     {
//                       p_fil = (pl == NULL ? NULL : pl->esq);
//                       fv_fil = ((pl == NULL || pl->esq == NULL) ? 2.0*fb - fint : pl->esq->fval);
//                     }
//                   fil[i*k+j].p = p_fil;
//                   fil[i*k+j].fv = fv_fil;
//                 }
//             }
//           else
//             {//filho esquerdo
//               for(j=0;j<k;j++)
//                 {
//                   int z = i/2;
//                   Reg *pl = pac[(z+1)*k+j].p;
//                   double fa = pac[z*k+j].fv;
//                   double fb = pac[(z+1)*k+j].fv;
//                   double fc = pac[(z+2)*k+j].fv;
//                   double fint = interpola_quad(fa, fb, fc);
//                   Reg *p_fil;
//                   double fv_fil;
//                   if(i%2==0)
//                     { 
//                       p_fil = (pl == NULL ? NULL : pl->esq);
//                       fv_fil = ((pl == NULL || pl->esq == NULL) ? 2.0*fb - fint : pl->esq->fval);
//                     }
//                   else
//                     {
//                       p_fil = (pl == NULL ? NULL : pl->dir);
//                       fv_fil = ((pl == NULL || pl->dir == NULL) ? fint : pl->dir->fval);
//                     }
//                   fil[i*k+j].p = p_fil;
//                   fil[i*k+j].fv = fv_fil;
//                 }
//             }
//         }
//     }
// }

bool_t poda_arvore(VReg pac[], double xmin[], double xmax[], int profund, Interp_quad interpola_quad, double fpoda, double eps)
{
  int  k=5; // posteriormente k pode se tornar argumento da funcao
  int ctr = (k*k)/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 = profund % 2;
  
  int qual;
  bool_t filho_superfluo[2]; /* TRUE para filhos que podem ser podados. */
  for (qual = 0; qual <= 1; qual++)
    { 
      VReg pac_fil[k*k];
      double xmin_fil[2], xmax_fil[2];
  
      /* Constrói o pacote do filho {qual}: */
			divide_pacote(pac, xmin, xmax, ex, interpola_quad, qual, pac_fil, xmin_fil, xmax_fil);
		
      /* Calcula o valor que {pac_fil[ctr].fv} teria se esse nó fosse NULL: */
      double fa, fb, fc, fpoda_fil;
      if (ex == 0)
        { 
          fa = pac[ctr-1].fv;
          fb = pac[ctr].fv;
          fc = pac[ctr+1].fv;
        }
      else
        {
          fa = pac[ctr-k].fv;
          fb = pac[ctr].fv;
          fc = pac[ctr+k].fv;
        }
      double fint = interpola_quad(fa, fb, fc);
      if (qual == 0)
        { fpoda_fil = 2.0*fb - fint; }
      else 
        { fpoda_fil = fint; }
		
			filho_superfluo[qual] = poda_arvore(pac_fil, xmin_fil, xmax_fil, profund+1, interpola_quad, fpoda_fil, eps);
    }
        
  /* Elimina os filhos, se s�o n�o-nulos mas ambos sup�rfluos: */
  if ((filho_superfluo[0] && filho_superfluo[1]))
    { /* Elimina os dois filhos: */
      free(pac[ctr].p->esq); pac[ctr].p->esq = NULL;
      free(pac[ctr].p->dir); pac[ctr].p->dir = NULL;
    }
          
  /* Paranoia: */
  assert((pac[ctr].p->esq == NULL) == (pac[ctr].p->dir == NULL));
        
  /* Se {pac[ctr].p} tem filhos, n�o � sup�rfluo: */
  if (pac[ctr].p->esq != NULL) { return FALSE; }

  /* Vale a pena manter a raiz? */
  double detalhe = fpoda - pac[ctr].fv;
  /* fprintf(stderr, "%*sdetalhe = %9.6f\n", 2*profund, "", detalhe); */

  if (fabs(detalhe) < eps)
    { /* 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;
    }
}


/**
 * Cria uma arvore completa. A media de {f} na célula é armazenada 
 * no nó correspondente.
 * @param niv [in] quantidade de niveis que tera a arvore
 * @param profund [in] profundidade de uma certa celula
 * @param ind [in] indice da celula
 * @param xmin [in] coordenada do vertice inferior esquerdo
 * @param xmax [in] coordenada do vertice superior direito
 * @param (f) funcao arbitraria
 */
Reg *completa(int niv, int prof, int indice, double xmin[], double xmax[], Funcao *f, Quadratura integra)
{ 
  if (niv == 0) { return NULL; }
  Reg *e = malloc(sizeof(Reg));

  if (f != NULL) 
    { e->fval = integra(xmin[0], xmax[0], xmin[1], xmax[1], f)/(fabs(xmax[0]-xmin[0])*fabs(xmax[1]-xmin[1])); }
  else
    { e->fval = 0.0; }
  e->ind = indice;

  double xmin_fil[2], xmax_fil[2];
  int ex = prof%2;
  if(ex==0)
    {
      double xmed = (xmin[ex] + xmax[ex])/2.0;
      xmin_fil[0] = xmin[0];
      xmax_fil[0] = xmed;
		
      xmin_fil[1] = xmin[1];
      xmax_fil[1] = xmax[1];
	
      e->esq = completa(niv - 1, prof + 1, 2*indice + 0, xmin_fil, xmax_fil, f,integra);
		
      xmin_fil[0] =  xmed;
      xmax_fil[0] = xmax[0];
		
      xmin_fil[1] = xmin[1];
      xmax_fil[1] = xmax[1];
      e->dir = completa(niv - 1, prof + 1, 2*indice + 1, xmin_fil, xmax_fil, f,integra);
    }
  else
    {
      double xmed = (xmin[1] + xmax[1])/2.0;
      xmin_fil[0] = xmin[0];
      xmax_fil[0] = xmax[0];
      xmin_fil[1] = xmin[1];
      xmax_fil[1] = xmed;
	
      e->esq = completa(niv - 1, prof + 1, 2*indice + 0, xmin_fil, xmax_fil, f,integra);
		
      xmin_fil[0] = xmin[0];
      xmax_fil[0] = xmax[0];
      xmin_fil[1] = xmed;
      xmax_fil[1] = xmax[1];
      e->dir = completa(niv - 1, prof + 1, 2*indice + 1, xmin_fil, xmax_fil, f,integra);
    }
	
  return e;

}

double interpola_ponto_4(VReg pac[], double xmin[], double xmax[], double x[])
{
  /* !!! Encontrar a f�rmula adequada !!! */
  int k=5; // Todo pacote tem {k*k} elementos (depende da ordem de interpolação). 
  int ctr = (k*k)/2;
  double hx = xmax[0] - xmin[0];
  double hy = xmax[1] - xmin[1];
  double a = xmin[0] - hx;
  double b = xmin[1] - hy;
  //aplica [xmin[qual]-h,xmax[qual]+h] em [-1,2]
  double z =(x[0] - a)/hx - 1.0;
  double p1x =  z*z*(z-1.0+1.0/4.0);    //pol no interv [-1,0]  p1 = (z + 1.0)*(z + 1.0)*(z + 1.0/4.0);
  double p2x =  15.0*z*z*z*z - 30.0*z*z*z + (27.0/2.0)*z*z + (3.0/2.0)*z + 1.0/4.0;
  double p3x =  -(z - 1.0)*(z - 1.0)*(z - 5.0/4.0+1.0);
  z = (x[1] - b)/hy - 1.0;
  double p1y =  z*z*(z-1.0+1.0/4.0);    //pol no interv [-1,0]  p1 = (z + 1.0)*(z + 1.0)*(z + 1.0/4.0);
  double p2y =  15.0*z*z*z*z - 30.0*z*z*z + (27.0/2.0)*z*z + (3.0/2.0)*z + 1.0/4.0;
  double p3y =  -(z - 1.0)*(z - 1.0)*(z - 5.0/4.0+1.0);  //pol no interv [1,2]  -(z - 2.0)*(z - 2.0)*(z - 5.0/4.0);
  printf("h=%f z=%f p1=%f p2=%f p3=%f\n",hx,z,p1x,p2x,p3x);
	
  return p2x*p2y*pac[ctr].fv + p3x*p1y*pac[ctr-1-k].fv + p2x*p1y*pac[ctr-k].fv + p1x*p1y*pac[ctr-k+1].fv + p1x*p2y*pac[ctr+1].fv + p1x*p3y*pac[ctr+k+1].fv + p2x*p3y*pac[ctr+k].fv + p3x*p3y*pac[ctr+k-1].fv + p3x*p2y*pac[ctr-1].fv;;
}

double interpola_ponto_c(VReg pac[], double xmin[], double xmax[], double x[])
{
  /* !!! Encontrar a f�rmula adequada !!! */
  int k=5; // Todo pacote tem {k*k} elementos (depende da ordem de interpolação). 
  int ctr = (k*k)/2;
  return pac[ctr].fv;
}