/* vide {base_aprox.h} */
#define base_aprox_C_COPYRIGHT "Copyright © 2008 Danillo Pereira e J. Stolfi, UNICAMP"
/* Last edited on 2008-07-22 12:15:23 by stolfi */ 

#define _GNU_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <assert.h>
#include <math.h>

#include <frgb.h>
#include <argparser.h>
#include <affirm.h>

#include "geometria_basica.h"
#include "matriz_esparsa.h"
#include "objeto.h"
#include "camera_e_filme.h"
#include "base_aprox.h"

BaseAprox CriaBaseAprox(TipoDeDistancia tpDist, TipoDeElemento tpElem, Sitio sb[], double rb[], int nsb)
  { 
    BaseAprox bas = BaseAprox_new(nsb);
    int i;
    for (i = 0; i < nsb; i++)
      { Elemento *eP = &(bas.e[i]); 
        eP->tpDist = tpDist;
        eP->tpElem = tpElem;
        eP->centro = sb[i];
        eP->raio = rb[i];
      }
    return bas;
  }

MatrizEsparsa ConstroiMatrizDaBase(BaseAprox *bas, Sitio sx[], int nsx, int m)
  {
    if (m < 1) { m = 1; }
    if (m > nsx) { m = nsx; }
    MatrizEsparsa M = AlocaMatrizEsparsa(nsx, bas->ne, m*bas->ne);

    double *v = NOTNULL(malloc(bas->ne * sizeof(double)));

    int i;
    int nenM = 0; //número de elementos usados em {M}.
    for (i = 0; i < nsx; i++)
      { //calcula a linha {i} da matriz {M}.
        CalculaValoresDaBase(bas, &(sx[i]), v);
        nenM = AcrescentaLinhaAMatrizEsparsa(&M, nenM, i, v);
      }
    AjustaTamanhoDeMatrizEsparsa(&M, nenM);
    return M;
  }
  
void CalculaValoresDaBase(BaseAprox *bas, Sitio *st, double v[])
  {
    int i;
    bool_t normaliza = FALSE; //diz se é necessário normalizar os valores para soma 1.
    for (i = 0; i < bas->ne; i++)
      { Elemento *eP = &(bas->e[i]); 
        v[i] = CalculaValorDeElemento(eP, st);
        if ((eP->tpElem == TipoDeElemento_Shepard) || (eP->tpElem == TipoDeElemento_RadialShepard))
          { normaliza = TRUE; }
      }
    if (normaliza)
      { NormalizaValoresParaSomaUm(v, bas->ne); }
  }

#define VALOR_MINUSCULO 1e-200
  /* Este valor é usado, em vez de zero, para inicializar somatórias
    de pesos, a fim de evitar divisão por zero. Assim, se todos os
    pesos forem nulos antes de normalizar, eles serão todos nulos (e
    não {NAN]s) depois de normalizados. Por outro lado, se algum peso
    tiver valor significativo (maior que {1e16*VALOR_MINUSCULO}), este
    truque não vai fazer nenhuma diferença no resultado. */
  
void NormalizaValoresParaSomaUm(double v[], int n)
  {
    double sum_v = VALOR_MINUSCULO;
    int ninf = 0; //número de pesos infinitos
    int i;
    for (i = 0; i < n; i++)
      { if (isnan(v[i]))
          { sum_v = NAN; }
        else if (isfinite(v[i]))
          { sum_v += v[i]; }
        else
          { demand(v[i] == +INF, "valor infinito inválido");
            sum_v = +INF; ninf++; 
          }
      }
    demand(ninf <= 1, "dois ou mais elementos com valor infinito no mesmo sítio");
    if (isnan(sum_v))
      { for (i = 0; i < n; i++) { v[i] = NAN; } }
    if (ninf > 0)
      { for (i = 0; i < n; i++) { v[i] = (v[i] == +INF ? 1.0 : 0.0); } }
    else
      { for (i = 0; i < n; i++) { v[i] /= sum_v; } }
  }
  
double CalculaValorDeElemento(Elemento *e, Sitio *st)
  {
    switch(e->tpElem)
      {
        case TipoDeElemento_Bolota:
          return CalculaValorDeElementoBolota(e, st);

        case TipoDeElemento_Radial:
          return CalculaValorDeElementoRadial(e, st);
          
        case TipoDeElemento_Shepard:
          return CalculaValorDeElementoShepard(e, st);

        case TipoDeElemento_RadialShepard:
          return CalculaValorDeElementoRadialShepard(e, st);

        default:
          assert(FALSE);
      }
  }

#define TOL_RAIO_SUPORTE 1.0e-6
  /* Fator de ajuste usado para garantir que todo elemento é zero na
    fronteira de seu suporte nominal, apesar de erros de
    arredondamento no cálculo da distância ao centróide. O elemento é
    nulo se a distância do centróide for maior que
    {1-TOL_RAIO_SUPORTE} vezes o raio nominal de suporte. */

double CalculaValorDeElementoBolota(Elemento *e, Sitio *st)
  {
    demand(e->tpElem == TipoDeElemento_Bolota, "tipo de elemento inválido");
    double dist = DistanciaEntreSitios(&(e->centro), st, e->tpDist);
    return (dist*(1 + TOL_RAIO_SUPORTE) >= e->raio ? 0.0 : 1.0);
  }

double CalculaValorDeElementoRadial(Elemento *e, Sitio *st)
  {
    demand(e->tpElem == TipoDeElemento_Radial, "tipo de elemento inválido");
    double dist = DistanciaEntreSitios(&(e->centro), st, e->tpDist);
    if (dist*(1 + TOL_RAIO_SUPORTE) >= e->raio)
      { return 0.0; }
    else 
      { //a função-mãe é o sino gaussiano com {sigma = e->raio/3}:
        double z = dist/(e->raio/3.0);
        return exp(-z*z/2);
      }
  }

double CalculaValorDeElementoShepard(Elemento *e, Sitio *st)
  {
    demand(e->tpElem == TipoDeElemento_Shepard, "tipo de elemento inválido");
    double dist = DistanciaEntreSitios(&(e->centro), st, e->tpDist);
    if (dist*(1 + TOL_RAIO_SUPORTE) >= e->raio)
      { return 0.0; }
    else if (dist == 0.0)
      { return +INF; }
    else
      { //a função-mãe é {1/dist^2} vezes o sino gaussiano com {sigma = e->raio/3}:
        double z = dist/(e->raio/3.0);
        return exp(-z*z/2)/(dist*dist);
      }
  }

double CalculaValorDeElementoRadialShepard(Elemento *e, Sitio *st)
  {
    demand(e->tpElem == TipoDeElemento_RadialShepard, "tipo de elemento inválido");
    double dist = DistanciaEntreSitios(&(e->centro), st, e->tpDist);
    if (dist*(1 + TOL_RAIO_SUPORTE) >= e->raio)
      { return 0.0; }
    else if (dist == 0.0)
      { return 1.0; }
    else
      { //a função-mãe é o sino gaussiano com {sigma = e->raio/3}:
        double z = dist/(e->raio/3.0);
        return exp(-z*z/2);
      }
  }

#define MAX_NUM_VIZINHOS 100
  /* Número máximo de vizinhos que deve ser tolerado dentro do suporte nominal de
    cada elemento, na escolha deste último. */

OpcoesDeBase *AnalisaOpcoesDeBase(argparser_t *pp)
  {
    OpcoesDeBase *opBas = NOTNULL(malloc(sizeof(OpcoesDeBase)));
    //----------------------------------------------------------------------
    if (argparser_keyword_present(pp, "-numVizinhos"))
      { opBas->numVizinhos = argparser_get_next_int(pp, 1, MAX_NUM_VIZINHOS); }
    else
      { opBas->numVizinhos = 10; }
    //----------------------------------------------------------------------
    return opBas;
  }

double EscolheRaioDoElemento(TipoDeDistancia tpDist, Sitio st[], int nst, int ind, int m) 
  {
    int maxCand = m+1; //número de sítios a manter na fila.
    
    //candidatos aos {maxCand} sítios mais próximos a {st[ind]}, sem contar o próprio: 
    int indCand[maxCand];     //indices dos candidatos relevantes.
    double distCand[maxCand]; //suas distâncias a {st[ind]}, em ordem crescente.

    int nCand = 0; //número de candidatos válidos na fila.
    
    auto void InsereNaFila(int ind_c, double dist_c);
    /* Insere {st[indc]} com distância {dist_c} na fila dos {m} melhores candidatos. */

    void InsereNaFila(int ind_c, double dist_c)
      { //procura o lugar {j} para inserir este candidato:
        int j = nCand;
        while ((j > 0) && (dist_c < distCand[j-1]))
          { //empurra o candidato {j-1} para posicao {j};
            //se {j} não existe, descarta o candidato {j-1}:
            if (j < maxCand)
              { indCand[j] = indCand[j-1];
                distCand[j] = distCand[j-1];
              }
            //o lugar vago agora é {j-1}:
            j--;
          }
        if (j < maxCand) 
          { //insere o novo candidato {st[i]} na posição {j}: 
            indCand[j] = ind_c;
            distCand[j] = dist_c;
          }
        //se a fila não estava cheia, ganahamos mais um candidato:
        if (nCand < maxCand) { nCand++; }
      }

    //percorre todos os sítios:
    Sitio *stA = &(st[ind]);
    int i;
    for(i = 0; i < nst; i++)
      { if (i != ind)
          { double d = DistanciaEntreSitios(stA, &(st[i]), tpDist);
            if ((d < +INF) && (d != 0.0)) { InsereNaFila(i, d); }
          }
      }
    return (nCand == 0 ? +INF: distCand[nCand-1]);
  }