/***************************************************************************
 *   Copyright (C) 2009 by Douglas Castro   *
 *   douglas@ime.unicamp.br   *
 *                                                                         *
 *   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 "inicializa.h" 
#include "assert.h" 

double condicaoinicial(int d, double x, double y, double z) 
{
  double xc, yc, dx, dy, r, dist;
  switch(d){
    case 1:
      if(x<-0.5) { return 0.0;}
      else if(x>0.5) { return 0.0;}
      else if(fabs(x)==0.5) { return 0.5;}
      else if(fabs(x)<0.5) { return 1.0;}
      break;
    case 2:
//       double xc = 0.500, yc = 0.500*M_SQRT1_2;
      xc = 0.0, yc = 0.0; // so pq o dominio vai ser [-1,1]x[-1,1]
      r = 0.2;
      /* Gaussiana centrada em {xc,yc} com raio medio {r}: */
      dx = x - xc;
      dy = y - yc;
      dist = dx*dx + dy*dy;
      if( dist< r*r)
      {
//     return exp(-dx*dx/r-dy*dy/r);
        return 1.0;
      }
      else
      { return 0.0; }
      break;
    case 3:
      return 1.0;
      break;
    default:
      printf("cheque a implementacao de condicao inicial em inicializa.c\n");
      assert(FALSE);
  }
}

void inicia(int d, VNo pac[], double xmin[], double xmax[], int profund, Preditor pred, int niv, int lim, Funcao *f, Quadratura integra, double eps)
{
  int ord = 5; // posteriormente ord pode se tornar argumento da funcao,eh o tamanho necessrio para interpolacao
  int ctr = pow(ord,d)/2;
  
  assert((pac[ctr].p->fil[0] == NULL) == (pac[ctr].p->fil[1] == NULL));
  
  if(profund == niv-1) 
  { /* Não há o que fazer: */
    return;
  }
  if(pac[ctr].p->leaf)
  { return; }
  
  int qual;
  if(pac[ctr].p->fil[0] == NULL) 
  {
    pac[ctr].p->fil[0] = malloc(sizeof(No)); // cria o filho
    pac[ctr].p->fil[1] = malloc(sizeof(No)); // cria o filho
    for(qual=0;qual<2;qual++)
    {
      free(pac[ctr].p->fil[qual]->fil[0]); pac[ctr].p->fil[qual]->fil[0] = NULL; // diz que o neto e null
      free(pac[ctr].p->fil[qual]->fil[1]); pac[ctr].p->fil[qual]->fil[1] = NULL; // diz que o neto e null
      pac[ctr].p->fil[qual]->deletable = FALSE;
      pac[ctr].p->fil[qual]->leaf = profund + 1 == niv ? TRUE : FALSE;
      pac[ctr].p->fil[qual]->virtual_leaf = /*profund + 1 == niv ? TRUE :*/ FALSE;
    }
  }
  for(qual=0;qual<2;qual++)
  {
    //calcula a posicao da celula filho
    int i;
    int ex = profund%d;
    double xmin_fil[d], xmax_fil[d];

    limites_celula( d, ex, qual, xmin, xmax, xmin_fil, xmax_fil);
    //calcular o valor para armazenar na celula filho
    double fint = integra( d, xmin_fil, xmax_fil, f);
    double cvol = 1.0;
    for(i=0;i<d;i++)
      cvol *= fabs(xmin_fil[i]-xmax_fil[i]);
    pac[ctr].p->fil[qual]->fval = fint/cvol;
    
    int corr = pow(ord,ex);
    if(profund>=4)
    {
      double x0 = pac[ctr-corr].fv; // considera a condicao de contorno, estipulada em divide_pacote
      double x1 = pac[ctr].fv; 
      double x2 = pac[ctr+corr].fv; 
      double prv = pred( x0, x1, x2 ); // previsao do filho esquerdo
      double fprv = qual==0 ? prv : 2.0*x1 - prv;
      double det = fabs(fint/cvol - fprv);
      double l = profund;
      double tol = pow(2.0,d*(l-(niv+1)))*eps;
      if(det<tol)
      {
        pac[ctr].p->fil[qual]->deletable = TRUE;
        if(pac[ctr].p->deletable)//nao precisa testar o filho pq acabamos de torna-lo deletavel
        {
          int i;
          bool_t pode = TRUE; // o filho pode ser folha virtual e o pai folha? 
          auto bool_t diz(VNo pac[], int ord, int d);
          bool_t diz(VNo pac[], int ord, int d)
          {
            bool_t temp = FALSE, tmp = TRUE;
            int i;
            for(i=0;i<d;i++)
            {
              int corr = pow(ord,i);
              if(pac[ctr - corr].p == NULL || pac[ctr - corr].p->fil[0] == NULL)
              { temp = TRUE; }
              else if(pac[ctr + corr].p == NULL || pac[ctr + corr].p->fil[1] == NULL)
              { temp = TRUE; } 
              else if(pac[ctr - corr].p->leaf || pac[ctr - corr].p->fil[0]->leaf)
              { temp = TRUE; }
              else if(pac[ctr + corr].p->leaf || pac[ctr + corr].p->fil[1]->leaf)
              { temp = TRUE; }
              else
              { tmp = FALSE; }
            }
            if(temp == TRUE)
            { return temp; }
            else{ return tmp; }
          }
          pode = diz(pac, ord, d);
          if(pode == TRUE)
          {
            pac[ctr].p->leaf = TRUE;
            pac[ctr].p->virtual_leaf = FALSE;
            pac[ctr].p->fil[qual]->leaf = FALSE;
            pac[ctr].p->fil[qual]->virtual_leaf = TRUE;
          }
          else
          {
            pac[ctr].p->leaf = (profund+1 == niv ? TRUE : FALSE); // se chegar no limite do nivel eh folha
            pac[ctr].p->virtual_leaf = (/*profund+1 == niv ? TRUE :*/ FALSE);
            pac[ctr].p->fil[qual]->leaf = (profund+1 == niv ? TRUE : FALSE);
            pac[ctr].p->fil[qual]->virtual_leaf = (profund+1 == niv ? TRUE : FALSE);
          }
          free(pac[ctr].p->fil[qual]->fil[0]); pac[ctr].p->fil[qual]->fil[0] = NULL;
          free(pac[ctr].p->fil[qual]->fil[1]); pac[ctr].p->fil[qual]->fil[1] = NULL;
        }
      }
    }
  }
  // nao avanca mais de um nivel por vez
  if(profund == lim)
  { return; }
  for (qual = 0; qual <= 1; qual++)
  { 
    //calcula a posicao da celula filho
    int ex = profund%d;
    double xmin_fil[d], xmax_fil[d];
    
    limites_celula( d, ex, qual, xmin, xmax, xmin_fil, xmax_fil);
    
    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);
    
    inicia(d, pac_fil, xmin_fil, xmax_fil, profund+1, pred, niv, lim, f, integra, eps);
  }
}

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;
  for(i=0;i<d;i++)
  {
    xmin_fil[i] = xmin[i];
    xmax_fil[i] = xmax[i];
  }
  xmax_fil[ex] = qual == 0 ? xmed : xmax_fil[ex];
  xmin_fil[ex] = qual == 1 ? xmed : xmin_fil[ex];
}

/**
 * Cria malha inicial...
 * @param d [in] dimensao do dominio
 * @param pac[] [in] pacote de vponteiros e medias celulares
 * @param xmin[] [in] limite inferior da celula filho
 * @param xmax[] [in] limite superior da celula filho
 * @param prof [in] profundidade do no na malha
 * @param pred [in] operador de previsao
 * @param niv [in] numero de niveis na malha
 * @param lim [in] limite na profundidade ate onde a funcao pode avancar
 * @param f [in] funcao a ser discretizada
 * @param integra [in] regra de integracao, quadratura de Gauss 5 pontos
 * @param eps [in] parametro de tolerancia
 * @return
 */
void inicia2(int d, VNo pac[], double xmin[], double xmax[], int prof, Preditor pred, int niv, int lim, Funcao *f, Quadratura integra, double eps)
{
  int ord = 5; // posteriormente ord pode se tornar argumento da funcao,eh o tamanho necessario para interpolacao
  int ctr = pow(ord,d)/2;
  
  if(prof==lim){return;}
  if(prof==niv){return;}
  assert((pac[ctr].p->fil[0] == NULL) == (pac[ctr].p->fil[1] == NULL));
  
  if((pac[ctr].p->fil[0] == NULL))
  {
    pac[ctr].p->fil[0] = malloc(sizeof(No)); // cria o filho
    pac[ctr].p->fil[1] = malloc(sizeof(No)); // cria o filho
    int qual;
    for(qual=0;qual<2;qual++)
    {
      pac[ctr].p->fil[qual]->fil[0] = NULL; // diz que o neto e null
      pac[ctr].p->fil[qual]->fil[1] = NULL; // diz que o neto e null
//       pac[ctr].p->fil[qual]->deletable = FALSE;
//       pac[ctr].p->fil[qual]->leaf = FALSE;
//       pac[ctr].p->fil[qual]->virtual_leaf = FALSE;
    }
  }
  
  int ex = prof % d;
  int qual;
  for(qual=0;qual<2;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);
    /* calcula os limites da celula filho */
    double xmin_fil[d]; double xmax_fil[d];
    limites_celula( d, ex, qual, xmin, xmax, xmin_fil, xmax_fil);
    //calcular o valor para armazenar na celula filho
    double fint = integra( d, xmin_fil, xmax_fil, f);
    double cvol = 1.0;
    int i;
    for(i=0;i<d;i++)
      cvol *= fabs(xmin_fil[i]-xmax_fil[i]);
    
    pac[ctr].p->fil[qual]->fval = fint/cvol;
    
    int corr = pow(ord,ex);
    double x0 = pac[ctr-corr].fv; // considera a condicao de contorno, estipulada em divide_pacote
    double x1 = pac[ctr].fv; 
    double x2 = pac[ctr+corr].fv; 
    double prv = pred( x0, x1, x2 ); // previsao do filho esquerdo
    double fprv = qual==0 ? prv : 2.0*x1 - prv;
    
    double det = fabs(fint/cvol - fprv);
    double l = prof;
    double tol = pow(2.0,d*(l-(niv+1)))*eps;
    if((det<tol)&&prof>3)
    {
      pac[ctr].p->fil[qual]->deletable = TRUE;
      pac[ctr].p->fil[qual]->leaf = pac[ctr].p->deletable ? TRUE : FALSE;
      pac[ctr].p->fil[qual]->virtual_leaf = FALSE;
    }
    else
    {
      pac[ctr].p->fil[qual]->deletable = FALSE;
      pac[ctr].p->fil[qual]->leaf = FALSE;
      pac[ctr].p->fil[qual]->virtual_leaf = FALSE;
    }
    inicia2( d, pac_fil, xmin_fil, xmax_fil, prof+1, pred, niv, lim, f, integra, eps);
  }
  
}