/* See SOTentFunction.h */
/* Last edited on 2003-05-07 17:52:57 by ra014852 */

#include <SOTentFunction.h>
#include <SOFunction.h>

#include <filefmt.h>
#include <fget.h>
#include <nget.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
    
#define T SOTentFunction


void SOTentFunction_EvalMth(T *f, double *p, double *fp);
void SOTentFunction_GradMth(T *f, double *p, double *dp);
void SOTentFunction_HessMth(T *f, double *p, double *ddp);
void SOTentFunction_WriteMth(T *f, FILE *wr);
SOTentFunction *SOTentFunction_CopyMth(T *f);

SOTentFunction_Methods *SOTentFunction_Methods_New(void);
SOTentFunction_Data *SOTentFunction_Data_New(void);
SOTentFunction *SOTentFunction_New(void);

SOTentFunction *SOTentFunction_FullNew(void);
  /* Allocates a new {SOTentFunction} object {f}
    (and its data record), with the proper {type}
    fields and methods record.*/

SOTentFunction *SOTentFunction_Cast(OBJ *f)
  { SOFunction *ff = (SOFunction *)f;
    if ((f != NULL) && isprefix(SOTentFunction_TypeId, ff->type))
      { return (SOTentFunction *)f; }
    else
      { return NULL; }
  }

#define SOTentFunction_FileFormat "2003-04-01"

/* OVERRIDES FOR PARENT CLASS METHODS */
  
void SOTentFunction_EvalMth(T *f, double *p, double *fp) 
  { double min[f->d->fn.pDim], max[f->d->fn.pDim];
    double tpoint[f->d->fn.pDim];
    SOGrid_Dim d = (SOGrid_Dim)f->d->fn.pDim;
    int i;
  
    assert(isprefix(SOTentFunction_TypeId, f->type), "type mismatch");

    SOGrid_cell_coords(f->d->index, f->d->rank, d, min, max);

    for(i = 0; i < d; i++) 
      { 
        tpoint[i] = p[i];
        tpoint[i] = (double)tpoint[i] * pow(2.0, (double)i/d); //Canonical geometry

        tpoint[i] = (double)((tpoint[i] - max[i]) / (max[i] - min[i]));
      }

    *fp = SOTent_eval(d, f->d->rank, tpoint);
  }
  /* Returns the SOTentFunction value relative to the (root cell) 
     considering its coordinates [0..2^{-i/d}].*/

void SOTentFunction_GradMth(T *f, double *p, double *dp) 
  { double min[f->d->fn.pDim], max[f->d->fn.pDim];
    double tpoint[f->d->fn.pDim];
    SOGrid_Dim d = (SOGrid_Dim)f->d->fn.pDim;
    int i;

    assert(isprefix(SOTentFunction_TypeId, f->type), "type mismatch");

    SOGrid_cell_coords(f->d->index, f->d->rank, d, min, max);

    for(i = 0; i < d; i++) 
      { 
        tpoint[i] = p[i];
        tpoint[i] = (double)tpoint[i] * pow(2.0, (double)i/d); //Canonical geometry

        tpoint[i] = (double)((tpoint[i] - max[i]) / (max[i] - min[i]));

        dp[i] = 1 / (max[i] - min[i]);

        dp[i] = (double)dp[i] * pow(2.0, (double)i/d); //Canonical geometry
      }

    SOTent_grad_eval(d, f->d->rank, tpoint, dp);
  }
  /* Returns the SOTentFunction gradient values relative to the (root cell) 
     considering its coordinates [0..2^{-i/d}].*/
  
void SOTentFunction_HessMth(T *f, double *p, double *ddp) 
  { assert(FALSE, "Hess not implemented"); }
  
T *SOTentFunction_CopyMth(T *f)
  { assert(isprefix(SOTentFunction_TypeId, f->type), "type/method bug");
    SOTentFunction *g = SOTentFunction_FullNew();
    g->d->fn.pDim = f->d->fn.pDim;
    g->d->fn.fDim = f->d->fn.fDim;
    assert(g->d->fn.fDim == 1, "tent functions should be scalar");
    g->d->index = g->d->index;
    g->d->rank = g->d->rank;
    return g;
  }

/* CLASS-SPECIFIC METHODS */
  
void SOTentFunction_WriteMth(T *f, FILE *wr)
  { assert(isprefix(SOTentFunction_TypeId, f->type), "type mismatch");
    
    filefmt_write_header(wr, "SOTentFunction", SOTentFunction_FileFormat);

    fprintf(wr, "index = %ld", f->d->index); fputc('\n', wr);
    fprintf(wr, "rank = %d", f->d->rank); fputc('\n', wr);

    filefmt_write_footer(wr, "SOTentFunction");
  }
    
/* OTHER PROCS */
  
SOTentFunction_Methods *SOTentFunction_Methods_New(void)
  { void *v = malloc(sizeof(SOTentFunction_Methods));
    return (SOTentFunction_Methods *)notnull(v, "no mem for SOTentFunction_Methods");
  }

SOTentFunction_Data *SOTentFunction_Data_New(void)
  { void *v = malloc(sizeof(SOTentFunction_Data));
    return (SOTentFunction_Data *)notnull(v, "no mem for SOTentFunction_Data");
  }

SOTentFunction *SOTentFunction_New(void)
  { void *v = malloc(sizeof(SOTentFunction));
    return (SOTentFunction *)notnull(v, "no mem for SOTentFunction");
  }

static SOTentFunction_Methods *TentMths = NULL;

SOTentFunction *SOTentFunction_FullNew(void)
  { SOTentFunction *f = SOTentFunction_New();
    f->type = SOTentFunction_TypeId;
    f->d = SOTentFunction_Data_New();
    if (TentMths == NULL)
      { TentMths = SOTentFunction_Methods_New();
        TentMths->fn.eval = (EvalMth *)&SOTentFunction_EvalMth;
        TentMths->fn.grad = (GradMth *)&SOTentFunction_GradMth;
        TentMths->fn.hess = (HessMth *)&SOTentFunction_HessMth;
	TentMths->fn.write = (WriteMth *)&SOFunction_WriteMth; /* NOT SOTent! */
        TentMths->fn.copy = (CopyMth *)&SOTentFunction_CopyMth;
        /* Class-specific methods */
        TentMths->write = (WriteMth *)&SOTentFunction_WriteMth;
      }
    f->m = TentMths;
    return f;
  }

SOTentFunction *SOTentFunction_Read(FILE *rd, SOGrid_Dim dDim)
  { 
    SOTentFunction *f = SOTentFunction_FullNew();
    
    filefmt_read_header(rd, "SOTentFunction", SOTentFunction_FileFormat);
    
    f->d->fn.pDim = dDim;
    f->d->fn.fDim = 1;
  
    f->d->index = nget_int(rd, "index"); fget_eol(rd);
    f->d->rank = nget_int(rd, "rank"); fget_eol(rd);

    filefmt_read_footer(rd, "SOTentFunction");
    return f;
  }
  
SOTentFunction *SOTentFunction_Make(SOGrid_Dim pDim, SOGrid_Index k, SOGrid_Rank r)
  { 
    SOTentFunction *f = SOTentFunction_FullNew();
    f->d->fn.pDim = pDim;
    f->d->fn.fDim = 1;
    f->d->index = k;
    f->d->rank = r;
    return f;
  }


double SOTentFunction_DotBoth(SOTentFunction *ftf, SOTentFunction *gtf)
  { SOTent_Pair tp;   
    SOGrid_Dim d;         /* Dimension of domain. */
    double sum = 0;       /* Accumulator for dot product. */
    double corr = 0;      /* Low-order bits of {sum}. */

    tp.cx[0] = ftf->d->index; 
    tp.cx[1] = gtf->d->index;  

    tp.r[0] = ftf->d->rank;
    tp.r[1] = gtf->d->rank;

    assert(ftf->d->fn.pDim == gtf->d->fn.pDim, "Tent Functions domain dimension must be equal");

    d = ftf->d->fn.pDim;

    SOTent_tt_dot(&tp, d, &sum, &corr);

    return(sum);
  }


double *SOTentFunction_DotListBoth(SOTentFunction *ftf, SOTentFunction *gtf)
  { SOTent_Pair tp;   
    SOGrid_Dim d = ftf->d->fn.pDim;         /* Dimension of domain. */
    int i, max_cells = 1 << d;
    double *sum;                            /* Accumulator for dot product. */

    sum = (double *)malloc(sizeof(double)*max_cells);
    for(i = 0; i < max_cells; i++) sum[i] = 0.0;

    tp.cx[0] = ftf->d->index; 
    tp.cx[1] = gtf->d->index;  

    tp.r[0] = ftf->d->rank;
    tp.r[1] = gtf->d->rank;

    assert(ftf->d->fn.pDim == gtf->d->fn.pDim, "Tent Functions domain dimension must be equal");

    SOTent_tt_dot_list(&tp, d, sum);

    /* Returns Null if all dot products are 0.0. */
    bool isNull = TRUE;
    for(i = 0; i < max_cells; i++) if(sum[i] != 0.0) isNull = FALSE;
    if(isNull){free(sum); sum = (double *)NULL;}

    return(sum);
  }
   

void SOTentFunction_DotSingle(SOTentFunction *tf, SOFunction *f, FuncMap *FMap, double *sum)
  { SOGrid_Dim pd, fd;          /* Dimension of domain and (f) general function range */
    double corr[f->d->fDim];    /* Low-order bits of {*sum}. */
    double min[f->d->pDim];     /* Lower coordinates of (tf) {tf's support} on root-cell. */
    double max[f->d->pDim];     /* Higher coordinates of (tf) {tf's support} on root-cell. */
    double mink[f->d->pDim], maxk[f->d->pDim], diff[f->d->pDim];
    int i, k, max_cells = (1 << f->d->pDim);

    pd = f->d->pDim;
    if(FMap->map == NULL) fd = f->d->fDim; else fd = FMap->vDim;
    for(i = 0; i < fd; i++){ sum[i] = 0.0; corr[i] = 0.0; }

    assert(pd == tf->d->fn.pDim, "Tent and General Functions domain dimension must be equal");
   
    auto void integrand(double *var, double *fvar);
    /* Integrand of dot product: computes {h(x) = f(x)*t(x)},
       where {t} is the tent function, and {x} is a point within
       the root cell (defined as [0..1] along each coordinate). */

    void integrand(double *x, double *fx)
    { double tx;
      int j;
    
      /* Evaluate the function {f} at {var}: */
      if(FMap->map == NULL)
        f->m->eval((SOFunction *)f, x, fx);
      else
        { double val[fd = f->d->fDim];
          f->m->eval((SOFunction *)f, x, val);
          FMap->map(val, x, fx);  
        }

      /* Evaluate the tent function {tf} at {var}: */
      tf->m->fn.eval((SOTentFunction *)tf, x, &tx);

      /* Multiply fvar by the tent value: */
      for(j = 0; j < fd; j++) { fx[j] *= tx; }
    }

    SOGrid_cell_coords(tf->d->index, tf->d->rank, pd, min, max); /* Find limits for the lower cell */
    for(i = 0; i < pd; i++) diff[i] = max[i] - min[i];
  
    for(k = 0; k < max_cells; k++)
      {
        for(i = 0; i < pd; i++)
          {
            mink[i] = min[i] + (double)((k >> i)&1)*diff[i];    
            mink[i] = (double)mink[i]/pow(2.0, (double)i/pd); //Canonical geometry

            maxk[i] = max[i] + (double)((k >> i)&1)*diff[i];
            maxk[i] = (double)maxk[i]/pow(2.0, (double)i/pd); //Canonical geometry
          }
        SOIntegral_Gauss_Limits(integrand, pd, fd, sum, corr, mink, maxk, 0);
      }
    
  }


double SOTentFunction_GradDotBoth(SOTentFunction *ftf, SOTentFunction *gtf)
  { SOTent_Pair tp;   
    SOGrid_Dim d;         /* Dimension of domain. */
    double sum = 0;       /* Accumulator for dot product. */
    double corr = 0;      /* Low-order bits of {sum}. */

    tp.cx[0] = ftf->d->index; 
    tp.cx[1] = gtf->d->index;  

    tp.r[0] = ftf->d->rank;
    tp.r[1] = gtf->d->rank;

    assert(ftf->d->fn.pDim == gtf->d->fn.pDim, "Tent Functions domain dimension must be equal");

    d = ftf->d->fn.pDim;

    SOTent_tt_grad_dot(&tp, d, &sum, &corr);

    return(sum);
  }


double *SOTentFunction_GradDotListBoth(SOTentFunction *ftf, SOTentFunction *gtf)
  { SOTent_Pair tp;   
    SOGrid_Dim d = ftf->d->fn.pDim;         /* Dimension of domain. */
    int i, max_cells = 1 << d;
    double *sum;                            /* Accumulator for dot product. */

    sum = (double *)malloc(sizeof(double)*max_cells);
    for(i = 0; i < max_cells; i++) sum[i] = 0.0;

    tp.cx[0] = ftf->d->index; 
    tp.cx[1] = gtf->d->index;  

    tp.r[0] = ftf->d->rank;
    tp.r[1] = gtf->d->rank;

    assert(ftf->d->fn.pDim == gtf->d->fn.pDim, "Tent Functions domain dimension must be equal");

    SOTent_tt_grad_dot_list(&tp, d, sum);

    /* Returns Null if all dot products are 0.0. */
    bool isNull = TRUE;
    for(i = 0; i < max_cells; i++) if(sum[i] != 0.0) isNull = FALSE;
    if(isNull){free(sum); sum = (double *)NULL;}

    return(sum);
  }


void SOTentFunction_GradDotSingle(SOTentFunction *tf, SOFunction *f, double *sum)
  { SOGrid_Dim pd, fd;          /* Dimension of domain and (f) general function range */
    double corr[f->d->fDim];    /* Low-order bits of {*sum}. */
    double min[f->d->pDim];     /* Lower coordinates of (tf) {tf's support} on root-cell. */
    double max[f->d->pDim];     /* Higher coordinates of (tf) {tf's support} on root-cell. */
    double mink[f->d->pDim], maxk[f->d->pDim], diff[f->d->pDim];
    int i, k, max_cells = (1 << f->d->pDim);

    pd = f->d->pDim;
    fd = f->d->fDim;
    for(i = 0; i < fd; i++){ sum[i] = 0.0; corr[i] = 0.0; }

    assert(pd == tf->d->fn.pDim, "Tent and General Functions domain dimension must be equal");

    /* DEBUG -- just in case we're testing 2 tent functions. */
    /* should be kept as comments.*/
       /*     SOTentFunction *dbf = SOTentFunction_Cast(f); */
       /*     if(dbf->d->index > tf->d->index){f = SOFunction_Cast(tf); tf = dbf;} */
    /* DEBUG */
   
    auto void integrand(double *var, double *fvar);
    /* Integrand of the gradient dot product: computes 
       {h(x) = h0(x) + h1(x) + ... + h(pd-1)(x), and hi(x) = gradi(f)*gradi(t)},
       where {t} is the tent function, and {x} is a point within
       the root cell. */

    void integrand(double *x, double *fx)
    { double dt[pd], df[pd*fd];
      int i, j;

      /* Evaluate the function {f} gradient {df} at {x}: */
      f->m->grad((SOFunction *)f, x, df);

      /* Evaluate the tent function {tf} gradient {dt} at {x}: */
      tf->m->fn.grad((SOTentFunction *)tf, x, dt);

      for(i = 0; i < fd; i++)
        {  
          fx[i] = 0;
 
          for(j = 0; j < pd; j++) 
            { fx[i] += dt[j]*df[(i*pd)+j]; }
        }
    }

    SOGrid_cell_coords(tf->d->index, tf->d->rank, pd, min, max); /* Find limits for the lower cell */
    for(i = 0; i < pd; i++) diff[i] = max[i] - min[i];

    for(k = 0; k < max_cells; k++)
      {
        for(i = 0; i < pd; i++)
          {
            mink[i] = min[i] + (double)((k >> i)&1)*diff[i];      
            mink[i] = (double)mink[i]/pow(2.0, (double)i/pd); //Canonical geometry
    
            maxk[i] = max[i] + (double)((k >> i)&1)*diff[i];
            maxk[i] = (double)maxk[i]/pow(2.0, (double)i/pd); //Canonical geometry
          }

        /* IMPORTANT: to run this with can. geometry this procedure must NOT do can. geo. correction! */ 
	SOIntegral_Gauss_Limits(integrand, pd, fd, sum, corr, mink, maxk, 0);
      }
    
  }


double SOTentFunction_Grad1DotBoth(SOTentFunction *ftf, SOTentFunction *gtf, int gradaxis)
  { SOTent_Pair tp;   
    SOGrid_Dim d;         /* Dimension of domain. */
    double sum = 0;       /* Accumulator for dot product. */
    double corr = 0;      /* Low-order bits of {sum}. */

    tp.cx[0] = ftf->d->index; 
    tp.cx[1] = gtf->d->index;  

    tp.r[0] = ftf->d->rank;
    tp.r[1] = gtf->d->rank;

    assert(ftf->d->fn.pDim == gtf->d->fn.pDim, "Tent Functions domain dimension must be equal");

    d = ftf->d->fn.pDim;

    SOTent_tt_grad1_dot(&tp, d, &sum, &corr, gradaxis);

    return(sum);
  }


double SOTentFunction_Integral(SOTentFunction *ftf)
  { double sum;
    int i, volcorr=0, fulldivs; /* Volume correction of integral relative to root cell */ 

    SOGrid_Rank rank = ftf->d->rank;
    SOGrid_Dim d = ftf->d->fn.pDim;

    fulldivs = rank/d;      
    for(i = 0; i < d; i++)
      volcorr += fulldivs + ((rank % d) > i ? 1 : 0);

    sum = (double) 1 / (1 << volcorr);
    return(sum);

  }