/* See msg_tent.h */
/* Last edited on 2009-08-23 19:31:21 by stolfi */

#include <math.h>
#include <assert.h>

#include <affirm.h>
#include <vec.h>

#include <bz_basic.h>
/* #include <bz_patch.h> */

#include <msg_pulse.h>
#include <msg_tent.h>

/* IMPLEMENTATIONS */

void msg_tent_eval
  ( msg_dim_t d, 
    msg_pulse_family_t fam[], 
    msg_pulse_t p[],
    interval_t R[],
    double x[],
    msg_degree_t ord[],
    double f[]
  )
  { double pf[msg_pulse_MAX_DEGREE + 1];
    f[0] = 1.0;
    int nf = 1;
    int i;
    for (i = 0; i < d; i++) 
      { interval_t *Ri = (R == NULL ? NULL : &(R[i]));
        msg_pulse_eval(&(fam[i]), &(p[i]), Ri, x[i], ord[i], pf);
        int j;
        for (j = 0; j < nf; j++)
          { int k; 
            for (k = 0; k <= ord[i]; k++)
              { f[k*nf + j] = f[j]*pf[k]; }
          }
        nf *= (ord[i] + 1);
      }
  }

void msg_tent_eval_total
  ( msg_dim_t d, 
    msg_pulse_family_t fam[], 
    msg_pulse_t p[],
    interval_t R[],
    double x[],
    msg_degree_t ord,
    double f[]
  )
  { double pf[msg_pulse_MAX_DEGREE + 1];
    f[0] = 1.0;
    int i;
    for (i = 0; i < d; i++) 
      { interval_t *Ri = (R == NULL ? NULL : &(R[i]));
        msg_pulse_eval(&(fam[i]), &(p[i]), Ri, x[i], ord, pf);
        affirm(FALSE, "not implemented");
      }
  }

msg_tent_vec_t msg_tents_from_cells
  ( msg_dim_t d, 
    msg_tent_mother_index_t tix,
    msg_cell_index_vec_t cell
  )
  { msg_tent_vec_t tv = msg_tent_vec_new(cell.ne);
    int nt = 0;
    int i; 
    for (i = 0; i < cell.ne; i++)
      { msg_cell_index_t ck = cell.e[i];
        tv.e[nt] = (msg_tent_t){tix,ck};
        nt++;
      }
    assert(nt == tv.ne);
    return tv;
  }

msg_tent_mother_index_t msg_tent_mother_index_max(msg_dim_t d, msg_pulse_family_t fam[])
  { msg_tent_mother_index_t tix = 0;
    int i;
    for (i = d-1; i >= 0; i--) { tix = tix * fam[i].nmp + (fam[i].nmp - 1); }
    return tix; 
  }

msg_tent_mother_index_t msg_tent_mother_index_pack
  ( msg_dim_t d,
    msg_pulse_family_t fam[],
    msg_pulse_mother_index_t pix[]
  )
  { msg_tent_mother_index_t tix = 0;
    int i;
    for (i = d-1; i >= 0; i--) { tix = tix * fam[i].nmp + pix[i]; }
    return tix;
  }

void msg_tent_mother_index_unpack
  ( msg_dim_t d, 
    msg_pulse_family_t fam[],
    msg_tent_mother_index_t tix,
    msg_pulse_mother_index_t pix[]
  )
  { int i;
    for (i = 0; i < d; i++) { pix[i] = tix % fam[i].nmp; tix /= fam[i].nmp; }
  }

/* PRINTOUT: */

void msg_tent_print
  ( FILE *wr, 
    msg_dim_t d,
    msg_pulse_family_t fam[], 
    msg_pulse_t p[]
  )
  { int i;
    /* Print pulse indices: */
    for (i = 0; i < d; i++) 
      { fprintf(wr, "[");
        msg_pulse_print(wr, &(fam[i]), &(p[i]));
        fprintf(wr, "]");
      }
  }

void msg_tent_packed_print
  ( FILE *wr, 
    msg_dim_t d,
    msg_pulse_family_t fam[], 
    msg_tent_t *t
  )
  { /* Unpack puse indices from tent: */
    msg_pulse_t p[d];
    msg_tent_unpack(d, fam, t, p);
    msg_tent_print(wr, d, fam, p);
  }

/* MANIPULATION OF TENT VECTORS */

vec_typeimpl(msg_tent_vec_t,msg_tent_vec,msg_tent_t);