/* -*- mode: C++ -*- */
//////////////////////////////////////////////////////////////////////////////
// FILE   : DesignC2.tol
// PURPOSE: Defines Class @DesignC2
//////////////////////////////////////////////////////////////////////////////


//////////////////////////////////////////////////////////////////////////////
Text _.autodoc.member.DesignC2.Monomial.Spheric.Grid = 
"";
Class @DesignC2.Monomial.Spheric.Grid : 
  @DesignC2.Monomial, 
  @DesignC2.Spheric
//////////////////////////////////////////////////////////////////////////////
{
  @IndexedGrid.Border _.grid;
//////////////////////////////////////////////////////////////////////////////
//Configuration parameters
//////////////////////////////////////////////////////////////////////////////
  Real _.radius;
  Real _.radius.margin;
  Real _.grid.size;

//////////////////////////////////////////////////////////////////////////////
//Auxiliar members
//////////////////////////////////////////////////////////////////////////////
  VMatrix _.Y.h = Rand(0,0,0,0);
  VMatrix _.Y.neighbourhood = Rand(0,0,0,0);

  ////////////////////////////////////////////////////////////////////////////
  Static @DesignC2.Monomial.Spheric.Grid New(
    VMatrix X, 
    Real radius.margin,
    Real radius,
    Real do.build.border)
  ////////////////////////////////////////////////////////////////////////////
  {
    @DesignC2.Monomial.Spheric.Grid new =
    [[ 
      VMatrix _.X = X;
      Real    _.n = VColumns(X);
      Real    _.J = VRows(X);
      Real    _.radius = radius;
      Real    _.radius.margin = radius.margin;
      Real    _.grid.size = radius/(1+radius.margin);
      @IndexedGrid.Border _.grid = @IndexedGrid.Border::New(
         X,_.grid.size,True,do.build.border,do.build.border)
    ]] 
  };

  //////////////////////////////////////////////////////////////////////////////
  VMatrix get.radius(VMatrix Y.idx)
  //////////////////////////////////////////////////////////////////////////////
  {
    Constant(VRows(Y.idx),1,_.radius)
  };

  //////////////////////////////////////////////////////////////////////////////
  VMatrix closest.reference.index(VMatrix X)
  //////////////////////////////////////////////////////////////////////////////
  {
    WriteLn("TRACE [closest.reference.index] 1");
    SubRow(_.Y.neighbourhood, Col2Set(_.grid::point2index(X)))
  };

  //////////////////////////////////////////////////////////////////////////////
  Real build.reference.points(Real void)
  //////////////////////////////////////////////////////////////////////////////
  {
    Text _MID = "[@DesignC2.Monomial.Spheric.Grid::Grid.reference.points] ";
    VMatrix _.Y := _.grid::_.Y;
    VMatrix _.Y.h := Mat2VMat(_.grid::_.Y.h);
    VMatrix _.Y.neighbourhood := _.Y.h | Group("ConcatColumns",
    For(1,2*_.n, VMatrix(Real i)
    {
      VMatrix aux = SubCol(_.grid::_.Y.neighbourhood.index,[[i]]);
      IfVMat(aux,aux,_.Y.h)
    }));
    Set X.class = Classify(Col2Set(_.grid::point2index(_.X)),Compare);
    VMatrix _.Yj := Mat2VMat(SetMat([[EvalSet(X.class,Card)]]),1);
    Real _.H := _.grid::_.H
  };

  //////////////////////////////////////////////////////////////////////////////
  Real build.stress.points(Real void)
  //////////////////////////////////////////////////////////////////////////////
  {
    Text _MID = "[@DesignC2.Spheric::build.stress.points] ";
    WriteLn("TRACE "+_MID+"1");
    Matrix Y.neighbour.triplet = Group("ConcatRows", For(1,_.H,Matrix(Real h)
    {
      Set neighbours = Unique(VRow2Set(SubRow(_.Y.neighbourhood,[[h]])));
      Set neighbours.up = Select(neighbours,Real(Real s) { s>h });
      Matrix s = SetCol(Sort(neighbours.up,Compare));
      Rand(Rows(s),1,h,h)|s|Rand(Rows(s),1,1,1)
    }));

    WriteLn("TRACE "+_MID+"2");
    VMatrix Y.neighbour.1.map = Triplet(Y.neighbour.triplet,_.H,_.H);
    VMatrix Y.neighbour.2.map = Y.neighbour.1.map*Y.neighbour.1.map;
    VMatrix Y.neighbour.12.map = Max(Y.neighbour.1.map,Y.neighbour.2.map);
    Matrix Y.neighbour.12.triplet = VMat2Triplet(Y.neighbour.12.map);
    VMatrix Z.h.a = Mat2VMat(SubCol(Y.neighbour.12.triplet,[[1]]));
    VMatrix Z.h.b = Mat2VMat(SubCol(Y.neighbour.12.triplet,[[2]]));
    VMatrix Y.neighbour.a = SubRow(_.Y,VCol2Set(Z.h.a));
    VMatrix Y.neighbour.b = SubRow(_.Y,VCol2Set(Z.h.b));
    VMatrix Z.all = (Y.neighbour.a + Y.neighbour.b)*0.5;

    WriteLn("TRACE "+_MID+"3");
    VMatrix Z.interior.a = 
    {
      VMatrix Z.norm2 = RPow(Y.neighbour.a-Z.all,2)*Rand(_.n,1,1,1);
      VMatrix Z.tao = RPow(Z.norm2,0.5)/_.radius;
      LE(Z.tao,Rand(VRows(Z.tao),1,1,1)) 
    };
    VMatrix Z.interior.b = 
    {
      VMatrix Z.norm2 = RPow(Y.neighbour.b-Z.all,2)*Rand(_.n,1,1,1);
      VMatrix Z.tao = RPow(Z.norm2,0.5)/_.radius;
      LE(Z.tao,Rand(VRows(Z.tao),1,1,1)) 
    };
    VMatrix match.interior = Z.interior.a $* Z.interior.b;

    Set sel.interior = MatQuery::SelectMatch(match.interior);
    VMatrix Z.interior = SubRow(Z.all,sel.interior);
    VMatrix Z.h.interior = SubRow(Z.h.a | Z.h.b,sel.interior);
    Set criterium = Range(1,_.n,1);
    Set sel.unique = Col2Set(Sort(VMat2Mat(Z.interior),criterium,True));
    WriteLn("TRACE "+_MID+"4");
    VMatrix _.Z := SubRow(Z.interior, sel.unique);
    VMatrix _.Z.h := SubRow(Z.h.interior, sel.unique);
    Real _.Tz := VRows(_.Z);
    VMatrix _.Z.Y.neighbour := 
    {
      Matrix Z.t = DifEq(1/(1-B),Rand(_.Tz,1,1,1));
      Matrix Z.th = (Z.t << Z.t) | 
       VMat2Mat((SubCol(_.Z.h,[[1]])<<SubCol(_.Z.h,[[2]])));
      Matrix Z.Y.idx = PivotByRows(Z.th,Sort(Z.th,[[2]]));
      VMatrix Z.Y.dif = 
       SubRow(_.Z,Col2Set(SubCol(Z.Y.idx,[[1]]))) -
       SubRow(_.Y,Col2Set(SubCol(Z.Y.idx,[[2]])));
        VMatrix Z.Y.norm2 = RPow(Z.Y.dif,2)*Rand(_.n,1,1,1);
        VMatrix Z.Y.tao = RPow(Z.Y.norm2,0.5)/_.radius;
        VMatrix Z.Y.match = LE(Z.Y.tao,Rand(VRows(Z.Y.tao),1,1,1));
      SubRow(Mat2VMat(Z.Y.idx),MatQuery::SelectMatch(Z.Y.match))
    };
    VMatrix aux = MatQuery::SortAndCount(SubCol(_.Z.Y.neighbour,[[2]]));
    VMatrix _.Yt := Triplet(VMat2Mat(SubCol(aux,[[1]])) |  
                            Rand(VRows(aux),1,1,1) | 
                            VMat2Mat(SubCol(aux,[[2]])),_.H,1);
    Real _.Mz := Case(
      prior.Z.order< 0, 0,
      prior.Z.order==0, _.Tz,
      prior.Z.order==1, _.Tz*(1+_.n),
      prior.Z.order==2, _.Tz*(1+_.n*3/2+(1/2)*_.n^2));
    _.Tz
  };

  //////////////////////////////////////////////////////////////////////////////
  Real build.border.points(Real void)
  //////////////////////////////////////////////////////////////////////////////
  {
    If(!Rows(_.grid::_.border.assign),{
      VMatrix _.Yb := Zeros(_.H,1);
      Real _.Tb := 0
    },{
    Text _MID = "[@DesignC2.Monomial.Spheric.Grid::build.border.points] ";
    VMatrix   B.t  = Mat2VMat(SubCol(_.grid::_.border.assign,[[1]]));
    VMatrix _.B.h := Mat2VMat(SubCol(_.grid::_.border.assign,[[2]]));
    VMatrix _.B := SubRow(_.grid::_.border::_.Y, VCol2Set(B.t));
    VMatrix aux = MatQuery::SortAndCount(_.B.h);
    VMatrix _.Yb := Triplet(VMat2Mat(SubCol(aux,[[1]])) |  
                            Rand(VRows(aux),1,1,1) | 
                            VMat2Mat(SubCol(aux,[[2]])),_.H,1);

    Real _.Tb := VRows(_.B);
    Real _.Mb := Case(
      prior.B.order< 0, 0,
      prior.B.order==0, _.Tb,
      prior.B.order==1, _.Tb*(1+_.n),
      prior.B.order==2, _.Tb*(1+_.n*3/2+(1/2)*_.n^2));
    _.Tb
    })
  };

  //////////////////////////////////////////////////////////////////////////////
    Real build.model.design(Real void)
  //////////////////////////////////////////////////////////////////////////////
  {
    Text _MID = "[@DesignC2.Monomial.Spheric.Grid::build.model.design] ";
    FunRn2R::@BaseC2 base.max = 
      FunRn2R::@BaseC2.Monomial::New(_.n,max.deg);
    Set Append(_.base, [[base.max]]);
    Real _.base.sel := Card(_.base);
    Real _.K0 := base.max::_.K;
    _build.model.design(void)
  };

  ///////////////////////////////////////////////////////////////////////////
    VMatrix generate.internal.points(
      Real m, 
      Real force.reference.near.nodal)
  ///////////////////////////////////////////////////////////////////////////
  {
    VMatrix h = Round(Rand(m,1,0.5,_.H+0.499999999));
    VMatrix y = SubRow(_.Y,VCol2Set(h));
    VMatrix v = Gaussian(m,_.n,0,1);
    VMatrix u = v*(VMatDat(MtMSqr(v),1,1)^-0.5);
    y+u*(Rand(0,1)*_.radius)
  };

  ////////////////////////////////////////////////////////////////////////////
    VMatrix match.evaluable(VMatrix X)
  ////////////////////////////////////////////////////////////////////////////
  {
    Real m = VRows(X);
    Mat2VMat(_.grid::point2index(X))
/*
    VMatrix X.closest.index = closest.reference.index(X);
    VMatrix X.h=SubRow(SubCol(_.Y.neighbourhood,[[1]]), 
                       VCol2Set(X.closest.index)); 
    VMatrix X.dif = X - SubRow(_.Y,VCol2Set(X.h));
    VMatrix X.norm2 = RPow(X.dif,2)*Rand(_.n,1,1,1);
    VMatrix X.tao = RPow(X.norm2,0.5)/_.radius;
    LE(X.tao,Rand(m,1,1,1))
*/
  }

};