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

#include "amr.h"

void mata_mato(Reg *r)
{
  if(r->esq->esq==NULL)
  {
    free(r->esq); r->esq=NULL;
    free(r->dir); r->dir=NULL;
    return;
  }
  mata_mato(r->esq);
  mata_mato(r->dir);
}

double trunca_detalhe(double detalhe);

double norma_inf(double vetexa[], double vetapr[],int niv)
{
  int comp = (int) pow(2,niv-1)/*(sizeof vetexa)/sizeof(double)*/;
  double norm = 0.0;
  for(comp-1;comp>=0;comp--)
  {
    double aux = fabs(vetexa[comp]-vetapr[comp]);
    if(aux > norm)
    {
      norm = aux;
    }
  }
  printf("contra prova norma = %f\n",norm);
  return norm;
}

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

void compara_folhas(Reg *r, double vet[], int *cont,  double *norm)
{
  assert((r->esq == NULL) == (r->dir == NULL));
  if(r->esq == 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->esq, vet, cont, norm);
  compara_folhas(r->dir, vet, cont, norm);
}


void copia_folhas(Reg *r, double vet[], int *cont)
{
  assert((r->esq == NULL) == (r->dir == NULL));
  if(r->esq == NULL)
  {
    int c = *cont;
//    printf("cccc = %d\n", c);
    double db = r->fval;
    vet[c] = db;
    *cont = c + 1;
    return;
  }
  copia_folhas(r->esq, vet, cont);
  copia_folhas(r->dir, vet, cont);
}

void amr_analise_haar(Reg *pac, int profund, double eps, int lim, bool_t sim)
{
  
  assert((pac->esq == NULL) == (pac->dir == NULL));
  
  if (pac->esq == NULL) 
  { /* Não há o que fazer: */
    return;
  }
  
  int qual;
  
  for (qual = 0; qual <= 1; qual++)
  { 
    if(qual==0)
    {
      amr_analise_haar(pac->esq, profund+1, eps, lim, sim);
    }
    else
    {
      amr_analise_haar(pac->dir, profund+1, eps, lim, sim);
    }
  }
  
  // lembrando: |C^k| = |C^{k-1}|/2
  //            |f^k| = (|f^{k+1}_e| + |f^{k+1}_d|)/2
  
  double restricao = ( pac->esq->fval + pac->dir->fval )/2.0;
  
  double detalhe = ( pac->esq->fval - pac->dir->fval )/2.0;
  //so para ver se preserva as bordas, e funciona
  double diff = pac->esq->fval - pac->dir->fval;
  
  pac->esq->fval = detalhe;
  pac->fval = restricao;
  
  /* 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
    if (fabs(detalhe) < eps && (pac->esq->esq==NULL) && fabs(diff)<eps)
    { /* Elimina os dois filhos: */
      free(pac->esq); pac->esq = NULL;
      free(pac->dir); pac->dir = NULL;
    }
  }
  
  return;
}

/**
 * Aqui o processo e um pouco mais delicado. tem que passar um parametro bool_t sim
 * para que a funcao saiba qdo criar novos nos na arvore, que foram retirados devido 
 * truncamento
 */
void amr_sintese_haar(Reg *pac, int profund, int niv)
{
  
  assert((pac->esq == NULL) == (pac->dir == NULL));
  if (profund == niv-1) 
  { /* Não há o que fazer: */
    return;
  }
  if (pac->esq == NULL) 
  { 
    Reg *l = malloc(sizeof(Reg));
    Reg *r = malloc(sizeof(Reg));
    pac->esq = l;
    pac->dir = r;
  }
  
  
  double detalhe = pac->esq->fval;
  pac->esq->fval = pac->fval + detalhe;
  pac->dir->fval = pac->fval - detalhe;
  
  int qual;
  
  for (qual = 0; qual <= 1; qual++)
  { 
    if(qual==0)
    {
      amr_sintese_haar(pac->esq, profund+1, niv);
    }
    else
    {
      amr_sintese_haar(pac->dir, profund+1, niv);
    }
  }
  
}

void amr_analise_harten(VReg pac[],int niv, double xmin[], double xmax[], int profund, Preditor pred, double eps, int lim, bool_t sim)    
{
  int  k=5; 
  int ctr = (k*k)/2;
  
  assert((pac[ctr].p->esq == NULL) == (pac[ctr].p->dir == NULL));
  
  if (pac[ctr].p->esq == NULL) 
  { /* Não há o que fazer: */
    return;
  }
  if( profund == lim )
  { /* tentativa de fazer analise um nivel por vez*/
    return;
  }
  /* Decide o eixo {ex} perpendicular ao corte: */
  int ex = profund % 2;
  
  int qual;
  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, pred, qual, pac_fil, xmin_fil, xmax_fil);

    amr_analise_harten(pac_fil, niv, xmin_fil, xmax_fil, profund+1, pred, eps, lim, sim);
  }
  
  // lembrando: |C^k| = |C^{k-1}|/2
  //            |f^k| = (|f^{k+1}_e| + |f^{k+1}_d|)/2
  
  double restricao = ( pac[ctr].p->esq->fval + pac[ctr].p->dir->fval )/2.0;
  double fa, fb, fc;
  
  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 fprev = pred(fa, fb, fc);/// por enquanto usamos a interplacao pra prever o filho esquerdo; 
  
  double detalhe = pac[ctr].p->esq->fval - fprev;
  double diff = pac[ctr].p->esq->fval - pac[ctr].p->dir->fval;
  
  //serao alteradas apenas os valores dos filhos das celulas de profundidade lim - 1
  if(profund == lim - 1 && pac[ctr].p->esq != NULL && pac[ctr].p->dir != NULL)
  { 
    pac[ctr].p->esq->fval = detalhe;
//    pac[ctr].p->fval = restricao;
    pac[ctr].p->dir->fval = restricao;
  }
  
  /* truncamento dos detalhes, se necessario */
  assert((pac[ctr].p->esq->esq==NULL) == (pac[ctr].p->esq->dir==NULL));
  assert((pac[ctr].p->dir->esq==NULL) == (pac[ctr].p->dir->dir==NULL));
  if(sim)
  {// se sim == TRUE, entra aqui 
    // o parametro de truncamento para medias celulares nao e constante
    // depende da profundidade na malha
    double pf = (double) profund;
    double q = 0.5; // paramentro no intervalo aberto (0,1) metade do usado na malha isotropica para fins de comparacao
    double kk = niv - pf;
    double eps_k = (1.0 - q)*eps*pow(q, kk);
    if((fabs(detalhe) < eps_k) && (pac[ctr].p->esq->esq==NULL) && (pac[ctr].p->dir->esq==NULL)/* && (fabs(diff)<eps)*/)
    { /* Elimina os dois filhos: */
      free(pac[ctr].p->esq); pac[ctr].p->esq = NULL;
      free(pac[ctr].p->dir); pac[ctr].p->dir = NULL;
    }
  }
  
  return;
}

void mostra_pac(FILE *arq, VReg 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(VReg pac[], double xmin[], double xmax[], int profund, Preditor pred, int niv, int lim)
{
  int  k=5; // posteriormente k pode se tornar argumento da funcao
  int ctr = (k*k)/2;
  
  assert((pac[ctr].p->esq == NULL) == (pac[ctr].p->dir == NULL));
  
  if (profund == niv-1) 
  { /* Não há o que fazer: */
    return;
  }
  
  
  // cria uma arvore, sera uasada se o filho for nulo
  Reg *e = malloc(sizeof(Reg));
  Reg *r = malloc(sizeof(Reg));
  if (pac[ctr].p->esq == NULL) 
  { 
    pac[ctr].p->esq = e;
    pac[ctr].p->esq->fval = 0.0;
    pac[ctr].p->esq->esq = NULL;
    pac[ctr].p->esq->dir = NULL;
  }
  if( pac[ctr].p->dir == NULL) 
  {
    pac[ctr].p->dir = r;
    pac[ctr].p->dir->fval = pac[ctr].fv;
    pac[ctr].p->dir->esq = NULL;
    pac[ctr].p->dir->dir = NULL;
  }
  
  /* Decide o eixo {ex} perpendicular ao corte: */
  int ex = profund % 2;
  
  double fa, fb, fc;
  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;
  }
  // a sintese sera feita a partir de onde a analise parou
  // lim e 
  if(profund ==  lim )///?????????
  {
    double detalhe = pac[ctr].p->esq->fval;
  
    double fprev = pred(fa, fb, fc);
    double esq = fprev + detalhe;
    pac[ctr].p->esq->fval = esq;
    double dir = 2.0*pac[ctr].p->dir->fval - esq;
    pac[ctr].p->dir->fval = dir;
    return;
  }
  int qual;
  //double fprev_fil;
  
  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, pred, qual, pac_fil, xmin_fil, xmax_fil);
    
    amr_sintese_harten(pac_fil, xmin_fil, xmax_fil, profund+1, pred, niv, lim);
  }
  
}