/* SOComputeMatrices - matrices of scalar products from function bases. */
/* Last edited on 2003-05-07 18:35:58 by ra014852 */

/* Computes the scalar product matrices {M[i,j] = <F[i],G[j]>}
  for two given bases {F[*]} and {G[*]} of functions.
  The dot product {<f,g>} is the integral over the root cell of 
  {f.eval(p)*geval(p)}, or {dot(f.grad(p),g.grad(p))}. */

#include <SOFunction.h> 
#include <SOBasic.h> 
#include <SOIntegral.h> 
#include <SOMatrix.h>
#include <SOBasisMatrices.h>
#include <SOParams.h>
#include <vec.h>
#include <js.h>
#include <math.h>
#include <string.h>
#include <limits.h>
#include <values.h>
#include <stdio.h>
#include <stdlib.h>

typedef struct Options
  { char *basisNameF;    /* File of basis {F} */
    char *basisNameG;    /* File of basis {G} */
    char *outName;       /* Prefix for output matrix files. */
    bool ignoreTrees;    /* TRUE ignores the grid when computing the dot product. */
    nat gaussOrder;      /* Gaussian quadrature order for dot products. */
    bool eval;           /* TRUE to write the {<eval|eval>} matrix. */
    bool grad;           /* TRUE to write the {<grad|grad>} matrix. */
    bool grad1;          /* TRUE to write the {<grad(gradaxis)|eval>} matrix. */
    nat gradaxis;        /* Axis where the differential shall be taken for a grad1 matrix. */
    double minVal;       /* Cleanup entries less than {minVal} times elem norm. */ 
    bool cholesky;       /* TRUE to write the Cholesky factor of {eval}. */
    bool checkAngles;    /* TRUE checks for close pairs in basis. */
    bool checkTriang;    /* TRUE checks the triangle inequality.  */
    bool verbose;        /* TRUE prints nonzero elements as they are computed. */
  } Options;

Options *GetOptions(int argn, char **argc);

Basis ReadBasis(char *name);
  /* Reads a function basis from the file {name} plus extension ".bas". */
  
void WriteMatrix(SOMatrix M, char *name); 
  /* Writes matrix {M} to the file {name} plus extension ".mat". */
  
void PrintSummary(SOMatrix M);
  /* Prints basic data about matrix {M}. */
  
void CheckAngles(SOMatrix M);
  /* Given the matrix {M} of dot products of a basis {B} against itself,
    checks whether {cos(B[i],B[j])} is close to 1.0 for
    any two distinct basis elements {B[i],B[j]}. */
     
void CheckTriangleInequality(SOMatrix M);
  /* Given the matrix {M} of dot products of a basis {B} against
    itself, checks whether it satisfies the triangle inequality
    {d(B[i],B[j]) \leq d(B[i],B[k])+d(B[k],B[j])}, where {d(f,g)} is
    computed from {M} by the formula 
    {d(f,g) = sqrt(<f|f> + <g|g> - 2*<f|g>)}. */
     
int main(int argn, char **argc)
  { Options *o = GetOptions(argn, argc);
    SOMatrix Meval, Mgrad, Mgrad1;
    Basis F = ReadBasis(o->basisNameF);
    Basis G = (strcmp(o->basisNameG, o->basisNameF) == 0 ? F : ReadBasis(o->basisNameG));
    
    SOIntegral_GaussOrder = o->gaussOrder;

    if (o->eval)
      { Meval = SOBasisMatrices_BuildMatrixEval
          ( F, G, NULL, o->ignoreTrees, o->minVal, o->verbose );
        PrintSummary(Meval);
        WriteMatrix(Meval, txtcat(o->outName, "-ev"));
        if (o->cholesky)
          { SOMatrix Mevch = SOMatrix_Cholesky(Meval, o->minVal);
            PrintSummary(Mevch);
            WriteMatrix(Mevch, txtcat(o->outName, "-ev-ch"));
          }
        if (o->checkAngles) { CheckAngles(Meval); }
        if (o->checkTriang) { CheckTriangleInequality(Meval); }
        free(Meval.ents.el);
      }
    if (o->grad)
      { Mgrad = SOBasisMatrices_BuildMatrixGrad
          ( F, G, NULL, o->ignoreTrees, o->minVal, o->verbose );
        PrintSummary(Mgrad);
        WriteMatrix(Mgrad, txtcat(o->outName, "-gr"));
        free(Mgrad.ents.el);
      }
    if (o->grad1)
      { Mgrad1 = SOBasisMatrices_BuildMatrixGrad1
          ( F, G, NULL, o->ignoreTrees, o->minVal, o->verbose, o->gradaxis );
        PrintSummary(Mgrad1);
        WriteMatrix(Mgrad1, txtcat(o->outName, "-gr1"));
        free(Mgrad1.ents.el);
      }
    return 0;
  }
  
Basis ReadBasis(char *name)
  { FILE *rd = open_read(txtcat(name, ".bas"));
    Basis F = SOFunction_ReadBasis(rd);
    fclose(rd);
    return F;
  }
  
void WriteMatrix(SOMatrix M, char *name)
  { FILE *wr = open_write(txtcat(name, ".mat"));
    SOMatrix_Write(wr, M);
    fclose(wr);
  }
  
void PrintSummary(SOMatrix M)
  { double tel = (double)(M.rows*M.cols);
    double nzel = (double)M.ents.nel;
    fprintf(stderr, "%5d rows\n", M.rows);
    fprintf(stderr, "%5d cols\n", M.cols);
    fprintf(stderr, "%8d entries\n", M.rows*M.cols);
    fprintf(stderr, "%8d nonzero entries\n", M.ents.nel);
    fprintf(stderr, "%8.6f occupancy fraction\n", nzel/tel);
    fprintf(stderr, "%8.6f mean entries per row\n", nzel/((double)M.rows));
    fprintf(stderr, "%8.6f mean entries per col\n", nzel/((double)M.cols));
  }

#define MAXCOS 0.9999
  /* Maximum value of {cos(B[i],B[j])} allowed in {checkAngles}. */

void CheckAngles(SOMatrix M)
  { MatEntry_vec Ments = M.ents;
    double_vec d = double_vec_new(M.rows);
    int i, k;
    
    assert(M.rows = M.cols, "matrix is not square");
    
    for (i = 0; i < M.rows; i++) { d.el[i] = 0.0; }
    for (k = 0; k < Ments.nel; k++)
      { MatEntry *Mij = &(Ments.el[k]);
        if (Mij->row == Mij->col) 
          { d.el[Mij->row] = sqrt(Mij->va); }
      }
    for (i = 0; i < M.rows; i++) 
      { assert(d.el[i] != 0.0, "diagonal element is zero"); }
    for (k = 0; k < Ments.nel; k++)
      { MatEntry *Mij = &(Ments.el[k]);
        nat i = Mij->row, j = Mij->col;
        if (i != j)
          { double c = Mij->va/d.el[i]/d.el[j];
            if (fabs(c) > MAXCOS)
              { fprintf(stderr, 
                  "basis elements %04d and %04d cos = %16.12f\n",
                  Mij->row, Mij->col, c);
              }
          }
      }
  }

void CheckTriangleInequality(SOMatrix M)
  { /* Dumb {O(N^3) log(N)} check - just to make sure. */
    int i, j, k;
    double eps = 1.0e-6;
    assert(M.rows = M.cols, "matrix is not square");
    
    for (i = 0; i < M.rows; i++)
      { for (j = 0; j < M.cols; j++)
          { for (k = 0; k < M.cols; k++)
              { double mii = SOMatrix_GetValue(M, i,i); 
                double mij = SOMatrix_GetValue(M, i,j);
                double mjj = SOMatrix_GetValue(M, j,j);
                double mik = SOMatrix_GetValue(M, i,k);
                double mjk = SOMatrix_GetValue(M, j,k);
                double mkk = SOMatrix_GetValue(M, k,k);
                double dij2 = mii - 2*mij + mjj;
                double dik2 = mii - 2*mik + mkk;
                double dkj2 = mkk - 2*mjk + mjj;
                double hij = sqrt(fabs(mii)+2*fabs(mij)+fabs(mjj));
                double hik = sqrt(fabs(mii)+2*fabs(mik)+fabs(mkk));
                double hkj = sqrt(fabs(mkk)+2*fabs(mjk)+fabs(mjj));
                double dij = sqrt(fabs(dij2));
                double dik = sqrt(fabs(dik2));
                double dkj = sqrt(fabs(dkj2));
                
                assert(dij2 >= -eps*hij, "distance dij2 negative");
                assert(dik2 >= -eps*hik, "distance dik2 negative");
                assert(dkj2 >= -eps*hkj, "distance dkj2 negative");
                assert(dij <= dik + dkj, "triangle cond failure");
              }
          }
      }
  }

#define PPUSAGE SOParams_SetUsage

Options *GetOptions(int argn, char **argc)
  {
    Options *o = (Options *)notnull(malloc(sizeof(Options)), "no mem");
    SOParams_T *pp = SOParams_NewT(stderr, argn, argc);

    PPUSAGE(pp, "SOComputeMatrices \\\n");
    PPUSAGE(pp, "  { -basis NAME | -bases FNAME GNAME } \\\n");
    PPUSAGE(pp, "  [ -eval [ -cholesky ] ] [ -grad ] [ -grad1 ] \\\n");
    PPUSAGE(pp, "  -gaussOrder NUM [ -minValue NUM ] [ -ignoreTrees ] \\\n");
    PPUSAGE(pp, "  [ -gradaxis (1:X, 2:Y, etc.) ] \\\n");
    PPUSAGE(pp, "  [ -checkAngles ] [ -checkTriang ] \\\n");
    PPUSAGE(pp, "  -outName NAME \\\n");
    PPUSAGE(pp, "  [ -verbose ] \n");

    if (SOParams_KeywordPresent(pp, "-bases"))
      { o->basisNameF = SOParams_GetNext(pp);  
        o->basisNameG = SOParams_GetNext(pp); 
        o->checkAngles = FALSE;
        o->checkTriang = FALSE;
       }
    else if (SOParams_KeywordPresent(pp, "-basis"))
      { o->basisNameF = SOParams_GetNext(pp);  
        o->basisNameG = o->basisNameF; 
        o->checkAngles = SOParams_KeywordPresent(pp, "-checkAngles");
        o->checkTriang = SOParams_KeywordPresent(pp, "-checkTriang");
      }
    else
      { SOParams_Error(pp, "cannot find neither \"-basis\" nor \"-bases\""); }
        
    o->eval = SOParams_KeywordPresent(pp, "-eval");
    o->grad = SOParams_KeywordPresent(pp, "-grad");
    o->grad1 = SOParams_KeywordPresent(pp, "-grad1");
    if (!(o->eval || o->grad || o->grad1))
      { SOParams_Error(pp, "cannot find neither \"-eval\", \"-grad\" nor \"-grad1\""); }
    if (o->eval)
      { o->cholesky = SOParams_KeywordPresent(pp, "-cholesky"); }
    else
      { o->cholesky = FALSE; }

    if (o->grad1)
      { SOParams_GetKeyword(pp, "-gradaxis");
        o->gradaxis = SOParams_GetNextInt(pp, 1, INT_MAX); o->gradaxis--;
      }
   
    if (SOParams_KeywordPresent(pp, "-minVal"))
      { o->minVal = SOParams_GetNextDouble(pp, 0.0, MAXDOUBLE); }
    else 
      { o->minVal = 0.0; }
    
    SOParams_GetKeyword(pp, "-outName");                               
    o->outName = SOParams_GetNext(pp);  

    SOParams_GetKeyword(pp, "-gaussOrder");
    o->gaussOrder = SOParams_GetNextInt(pp, 1, INT_MAX);

    o->ignoreTrees = SOParams_KeywordPresent(pp, "-ignoreTrees");
    
    o->verbose = SOParams_KeywordPresent(pp, "-verbose");

    SOParams_Finish(pp);
    return o;
  }