/* SOMultiLevelMatrices -- Finds all minimal or maximal basis for each
  level from a regular dyadic grid of rank A to a regular dyadic grid
  of rank B, and writes dot product matrices with all combinations of
  those basis.  */

/* Last edited on 2003-05-07 18:16:57 by ra014852 */

#include <SOGrid.h>
#include <SOFunction.h>
#include <SOTentFunction.h>
#include <SOMatrix.h>
#include <SOParams.h>
#include <SOBasic.h>
#include <js.h>
#include <stdio.h>
#include <values.h>
 
typedef struct Options
  { char *outName;    /* Name of output file (minus the ".bas" extension) */
    SOGrid_Dim pDim;  /* Grid dimension. */
    bool maximal;     /* TRUE for maximal-support tents, FALSE for minimal ones. */
    int minRank;      /* Minimum rank of the grids. */
    int maxRank;      /* Maximum rank of the grids. */
  } Options;
  
/* INTERNAL PROTOTYPES */

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

void WriteMatrix(SOMatrix M, char *name);

/* IMPLEMENTATIONS */

int main(int argn, char **argc)
  {
    Options *o = GetOptions(argn, argc);
    int i, j, k, levels = o->maxRank - o->minRank + 1;
    double val;
    SOGrid_Tree *trees[levels];
    SOTent_vec tents[levels];
    Basis allbasis;
    SOMatrix Meval, Mgrad;
    SOGrid_Dim pDim = o->pDim;

    auto void ShatterNode( SOGrid_Node *p, SOGrid_Rank r, SOGrid_Rank maxr );

    void ShatterNode( SOGrid_Node *p, SOGrid_Rank r, SOGrid_Rank maxr )
      { if( r < maxr )
	  { 
            SOGrid_Node_split(p);
            ShatterNode( p->c[LO], r+1, maxr );
            ShatterNode( p->c[HI], r+1, maxr );
	  }
      } 

    for(i = 0; i < levels; i++) 
      { 
        trees[i] = SOGrid_Tree_new(pDim);
        ShatterNode(trees[i]->root, 0, o->minRank + i);
      }

    j = 0;
    if (o->maximal)
      for(i = 0; i < levels; i++){ tents[i] = SOTent_maximal_basis(trees[i]); j += tents[i].nel; }
    else 
      for(i = 0; i < levels; i++){ tents[i] = SOTent_minimal_basis(trees[i]); j += tents[i].nel; }

    allbasis = SOFunctionRef_vec_new(j);

    k = 0; 
    for(i = 0; i < levels; i++) 
      for(j = 0; j < tents[i].nel; j++)
        { SOTent t = tents[i].el[j];

	  allbasis.el[k] = (SOFunction *)SOTentFunction_Make(pDim, t.cx, t.r);
          k++;
        }

    Meval.ents = MatEntry_vec_new((allbasis.nel * allbasis.nel + allbasis.nel) / 2);
    Meval.ents.nel = 0;

    Mgrad.ents = MatEntry_vec_new((allbasis.nel * allbasis.nel + allbasis.nel) / 2);
    Mgrad.ents.nel = 0;

    int maxcolev = 0, maxcolgr = 0, maxcolevi, maxcolgri;

    for(i = 0; i < allbasis.nel; i++)
      { 
        maxcolevi = 0;
        maxcolgri = 0;

        for(j = 0; j < allbasis.nel; j++)
          {
            SOTentFunction *itf = SOTentFunction_Cast(allbasis.el[i]);
            SOTentFunction *jtf = SOTentFunction_Cast(allbasis.el[j]);

            if( itf->d->index >= jtf->d->index )
              {                   
                val = SOTentFunction_DotBoth(itf, jtf);
                if( val != 0.0 )
                  {
                    Meval.ents.el[Meval.ents.nel].va = val;
                    Meval.ents.el[Meval.ents.nel].row = itf->d->index;
                    Meval.ents.el[Meval.ents.nel].col = jtf->d->index;
                    Meval.ents.nel++;
                    maxcolevi++;
	          }

                val = SOTentFunction_GradDotBoth(itf, jtf);
                if( val != 0.0 )
                  {
                    Mgrad.ents.el[Mgrad.ents.nel].va = val;
                    Mgrad.ents.el[Mgrad.ents.nel].row = itf->d->index;
                    Mgrad.ents.el[Mgrad.ents.nel].col = jtf->d->index;
                    Mgrad.ents.nel++;
                    maxcolgri++;
  	          }
              }   
          }

        if(maxcolevi > maxcolev) maxcolev = maxcolevi;
        if(maxcolgri > maxcolgr) maxcolgr = maxcolgri;

      }

    Meval.rows = allbasis.nel; Meval.cols = maxcolev;
    Mgrad.rows = allbasis.nel; Mgrad.cols = maxcolgr;

    WriteMatrix(Meval, txtcat(o->outName, "-ev"));
    WriteMatrix(Mgrad, txtcat(o->outName, "-gr"));

    return 0;
  }
  

#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, "SOMultiLevelMatrices \\\n");
    PPUSAGE(pp, "  -outName MATRIXNAME\n");
    PPUSAGE(pp, "  -pDim D \\\n");
    PPUSAGE(pp, "  [ -maximal | -minimal ] \\\n");
    PPUSAGE(pp, "  -minRank NUM \\\n");
    PPUSAGE(pp, "  -maxRank NUM \\\n");

    SOParams_GetKeyword(pp, "-outName");                               
    o->outName = SOParams_GetNext(pp);  

    SOParams_GetKeyword(pp, "-pDim");                               
    o->pDim = SOParams_GetNextInt(pp, 1, 4);  

    if (SOParams_KeywordPresent(pp, "-maximal"))
      { o->maximal = TRUE; }
    else if (SOParams_KeywordPresent(pp, "-minimal"))
      { o->maximal = FALSE; }
    else
      { SOParams_Error(pp, "must specify \"-maximal\" or \"-minimal\""); }

    SOParams_GetKeyword(pp, "-maxRank");                               
    o->maxRank = SOParams_GetNextInt(pp, 0, SOGRID_MAX_RANK);  

    SOParams_GetKeyword(pp, "-minRank");                               
    o->minRank = SOParams_GetNextInt(pp, 0, SOGRID_MAX_RANK);  

    SOParams_Finish(pp);
    return o;
  }


void WriteMatrix(SOMatrix M, char *name)
  { FILE *wr = open_write(txtcat(name, ".mat"));
    SOMatrix_Write(wr, M);
    fclose(wr);
  }