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


//////////////////////////////////////////////////////////////////////////////
Text _.autodoc.member.DesignC2 = 
"";
Class @DesignC2
//////////////////////////////////////////////////////////////////////////////
{
//Nodal points
  VMatrix _.X;
  Real _.n;
  Real _.J;

//////////////////////////////////////////////////////////////////////////////
//Configuration parameters
//////////////////////////////////////////////////////////////////////////////
  Real prior.Z.order = 2;
  Real min.contrast.surface = 3;

//////////////////////////////////////////////////////////////////////////////
//Auxiliar members
//////////////////////////////////////////////////////////////////////////////

  Real _.X.maxNeighbour = ?;
  Real _.Z.maxNeighbour = ?;
  VMatrix _.Y = Rand(0,2,0,0);
  VMatrix _.Z = Rand(0,2,0,0);
  VMatrix _.Z.h = Rand(0,1,0,0);
  VMatrix _.B = Rand(0,2,0,0);
  VMatrix _.B.h = Rand(0,0,0,0);
  Real _.K = 0;
  Real _.H = 0;
  Real _.Tz = 0;
  Real _.Mz = 0;
  Real _.Tb = 0;
  Real _.Mb = 0;
  Real _.M = 0;
  Real _.N = 0;
  Set _.base = Copy(Empty);
  Real _.base.sel = ?;
  VMatrix _.QX = Rand(0,0,0,0);
  VMatrix _.QB = Rand(0,0,0,0);
  VMatrix _.PZ = Rand(0,0,0,0);
  VMatrix _.Q = Rand(0,0,0,0);
  VMatrix _.LQ = Rand(0,0,0,0);
  VMatrix _.Q.nonNullCellsByColumns = Rand(0,0,0,0);
  //////////////////////////////////////////////////////////////////////////////
  NameBlock get.base(Real sel)
  //////////////////////////////////////////////////////////////////////////////
  {
    Real sel_ = If(Or(IsUnknown(sel),sel<=0,sel>Card(_.base)),_.base.sel);
    _.base[sel_]
  };

  //////////////////////////////////////////////////////////////////////////////
  //Virtual methods
  //////////////////////////////////////////////////////////////////////////////
  Real build.base(Real void);
  Real build.reference.points(Real void);
  Real build.stress.points(Real void) { False };
  Real build.border.points(Real void) { False };
  VMatrix closest.reference.index(VMatrix X);
  Set get.X.weight(Set X.dif);
  VMatrix match.evaluable(VMatrix X);

  //////////////////////////////////////////////////////////////////////////////
    Real _build.model.design(Real void)
  //////////////////////////////////////////////////////////////////////////////
  {
    Real _.X.maxNeighbour := 1+2*_.n;
    Real _.Z.maxNeighbour := If(prior.Z.order>=0,2,0);
    Real build.reference.points(0);
    Real build.stress.points(0);
    Real build.border.points(0);
    Real _.M := _.J+_.Mb+_.Mz;
    Real build.base(0);
    Real _.N := _.H*_.K;
    _.N<_.M
  };

  //////////////////////////////////////////////////////////////////////////////
    Real build.model.design(Real void)
  //////////////////////////////////////////////////////////////////////////////
  {
    _build.model.design(void)
  };

  //////////////////////////////////////////////////////////////////////////////
    VMatrix build.obs.eval(VMatrix X, VMatrix X.h, FunRn2R::@BaseC2.Eval eval)
  //////////////////////////////////////////////////////////////////////////////
  {
    Text _MID = "[@DesignC2::build.obs] ";
    Real m = VRows(X);
    WriteLn("TRACE "+_MID+"1 m="<<m);
    If(!m, Rand(0,_.N,0,0), {
    WriteLn("TRACE "+_MID+"2");
    WriteLn("TRACE "+_MID+"2");
    VMatrix X.j=DifEq(1/(1-B),Rand(m,1,1,1));
    WriteLn("TRACE "+_MID+"2");
    Set X.dif = For(1,_.X.maxNeighbour,VMatrix(Real i)
    {
      Set X.h.i = VCol2Set(SubCol(X.h,[[i]]));
      X - SubRow(_.Y,X.h.i)
    });
    WriteLn("TRACE "+_MID+"2");
    Set X.weight = get.X.weight(X.dif);
    WriteLn("TRACE "+_MID+"2");
    Matrix QX.triplet = Group("ConcatRows", For(1,_.X.maxNeighbour,Matrix(Real i)
    {
      VMatrix w.i = X.weight[i];
      Set sel.rows = MatQuery::SelectMatch(NE(Drop(Abs(w.i),1.E-15),w.i*0));
      VMatrix wi = SubRow(w.i,sel.rows);
      VMatrix X.hiK = SubRow((SubCol(X.h,[[i]])-1)*_.K,sel.rows);
      VMatrix Xdif.i = SubRow(X.dif[i],sel.rows);
      VMatrix X.j.i = SubRow(X.j,sel.rows);
      Real eval::upgrade.cache(Xdif.i);
      Group("ConcatRows", For(1,_.K,Matrix(Real k)
      {
        VMatrix pk = eval::evaluate.v(k,Xdif.i);
        VMatrix qk = wi$*pk; 
        VMat2Mat(X.j.i | (X.hiK+k) | qk)
      }))
    }));
    Pack(Convert(Triplet(QX.triplet, m, _.N),"Cholmod.R.Sparse"))
  })};

  //////////////////////////////////////////////////////////////////////////////
    VMatrix build.obs(VMatrix X, VMatrix X.h)
  //////////////////////////////////////////////////////////////////////////////
  {
    FunRn2R::@BaseC2.Eval.Function eval = 
     FunRn2R::@BaseC2.Eval.Function::New(get.base(0));
    build.obs.eval(X,X.h,eval)
  };
  //////////////////////////////////////////////////////////////////////////////
    VMatrix build.obs.partial.derivative.1(VMatrix X, VMatrix X.h, Real s)
  //////////////////////////////////////////////////////////////////////////////
  {
    FunRn2R::@BaseC2.Eval.Partial.Derivative.1 eval = 
     FunRn2R::@BaseC2.Eval.Partial.Derivative.1::New(get.base(0), s);
    build.obs.eval(X,X.h,eval)
  };
  //////////////////////////////////////////////////////////////////////////////
    VMatrix build.obs.partial.derivative.2(VMatrix X, VMatrix X.h, 
                                           Real s1, Real s2)
  //////////////////////////////////////////////////////////////////////////////
  {
    FunRn2R::@BaseC2.Eval.Partial.Derivative.2 eval = 
     FunRn2R::@BaseC2.Eval.Partial.Derivative.2::New(get.base(0), s1, s2);
    build.obs.eval(X,X.h,eval)
  };
  //////////////////////////////////////////////////////////////////////////////
    VMatrix build.obs.partial.integrate.1(VMatrix X, VMatrix X.h, Real s)
  //////////////////////////////////////////////////////////////////////////////
  {
    FunRn2R::@BaseC2.Eval.Partial.Integrate.1 eval = 
     FunRn2R::@BaseC2.Eval.Partial.Integrate.1::New(get.base(0), s);
    build.obs.eval(X,X.h,eval)
  };
  //////////////////////////////////////////////////////////////////////////////
    VMatrix build.border.prior(Real void)
  //////////////////////////////////////////////////////////////////////////////
  {
    If(!_.Tb, Rand(0,_.N,0,0), {
    VMatrix B.t=DifEq(1/(1-B),Rand(_.Tb,1,1,1));
    VMatrix Y = SubRow(_.Y, VCol2Set(_.B.h));
    VMatrix B.dif = _.B - Y;
    FunRn2R::@BaseC2 b = get.base(0);    
    Matrix PB.triplet = 
    {
      VMatrix B.hK = (_.B.h-1)*_.K;
      Real b::upgrade.cache(B.dif,prior.Z.order);
      Real lastRow = 0;
      Group("ConcatRows", For(1,_.K,Matrix(Real k)
      {
        Real cumRows = lastRow;
        VMatrix pk = b::function.v(k,B.dif);
        VMatrix col.idx = (B.hK+k);
        {
          VMatrix row.idx = B.t+cumRows;
          Real cumRows := cumRows + _.Tb;
          VMat2Mat(row.idx | col.idx | pk)
        }  << If(prior.Z.order<1,Rand(0,3,0,0),
        Group("ConcatRows",For(1,_.n,Matrix(Real d)
        {
          VMatrix pkd = b::partial.derivative.1.v(k,B.dif,d);
          VMatrix row.idx = B.t+cumRows;
          Real cumRows := cumRows + _.Tb;
          VMat2Mat(row.idx | col.idx | pkd)
        })))<< If(prior.Z.order<2,Rand(0,3,0,0),
        Group("ConcatRows",For(1,_.n,Matrix(Real d1)
        {
          Group("ConcatRows",For(1,d1,Matrix(Real d2)
          {
            VMatrix pkd2 = b::partial.derivative.2.v(k,B.dif,d1,d2);
            VMatrix row.idx = B.t+cumRows;
            Real cumRows := cumRows + _.Tb;
            VMat2Mat(row.idx | col.idx | pkd2)
          }))
        })))
      }))
    };
    Convert(Triplet(PB.triplet, _.Mb, _.N),"Cholmod.R.Sparse")
  })};
  //////////////////////////////////////////////////////////////////////////////
    VMatrix build.stress.prior(Real void)
  //////////////////////////////////////////////////////////////////////////////
  {
    If(!_.Tz, Rand(0,_.N,0,0), {
    VMatrix Z.t=DifEq(1/(1-B),Rand(_.Tz,1,1,1));
    
    Set Z.dif = For(1,_.Z.maxNeighbour,VMatrix(Real i)
    {
      _.Z - SubRow(_.Y,VCol2Set(SubCol(_.Z.h,[[i]])))
    });
    FunRn2R::@BaseC2 b = get.base(0);    
    Matrix PZ.triplet = Group("ConcatRows", For(1,_.Z.maxNeighbour,Matrix(Real i)
    {
      VMatrix Z.hiK = (SubCol(_.Z.h,[[i]])-1)*_.K;
      VMatrix Zdif.i = Z.dif[i];
      VMatrix Z.t.i = Z.t;
      Real b::upgrade.cache(Zdif.i,prior.Z.order);
      Real lastRow = 0;
      Group("ConcatRows", For(1,_.K,Matrix(Real k)
      {
        Real cumRows = lastRow;
        VMatrix pk = b::function.v(k,Zdif.i)*If(i==1,+1,-1);
        VMatrix col.idx = (Z.hiK+k);
        {
          VMatrix row.idx = Z.t.i+cumRows;
          Real cumRows := cumRows + _.Tz;
          VMat2Mat(row.idx | col.idx | pk)
        }  << If(prior.Z.order<1,Rand(0,3,0,0),
        Group("ConcatRows",For(1,_.n,Matrix(Real d)
        {
          Real sign = If(i==1,+1/2,-1/2);
          VMatrix pkd = b::partial.derivative.1.v(k,Zdif.i,d)*sign;
          VMatrix row.idx = Z.t.i+cumRows;
          Real cumRows := cumRows + _.Tz;
          VMat2Mat(row.idx | col.idx | pkd)
        })))<< If(prior.Z.order<2,Rand(0,3,0,0),
        Group("ConcatRows",For(1,_.n,Matrix(Real d1)
        {
          Group("ConcatRows",For(1,d1,Matrix(Real d2)
          {
            Real sign = If(i==1,+1/4,-1/4);
            VMatrix pkd2 = b::partial.derivative.2.v(k,Zdif.i,d1,d2)*sign;
            VMatrix row.idx = Z.t.i+cumRows;
            Real cumRows := cumRows + _.Tz;
            VMat2Mat(row.idx | col.idx | pkd2)
          }))
        })))
      }))
    }));
    Convert(Triplet(PZ.triplet, _.Mz, _.N),"Cholmod.R.Sparse")
  })};
  //////////////////////////////////////////////////////////////////////////////
    Real build.decomposition(Real void)
  //////////////////////////////////////////////////////////////////////////////
  {
    Text _MID = "[@DesignC2::build.decomposition] ";
    WriteLn("TRACE "+_MID+"1");
    VMatrix X.h=closest.reference.index(_.X);
    VMatrix _.QX := build.obs(_.X, X.h);
    WriteLn("TRACE "+_MID+"2");
    VMatrix _.QB := build.border.prior(0);
    WriteLn("TRACE "+_MID+"3");
    VMatrix _.PZ := build.stress.prior(0);
    WriteLn("TRACE "+_MID+"4");
    VMatrix _.Q := _.QX << _.QB << _.PZ;
    WriteLn("TRACE "+_MID+"7");
    VMatrix _.Q.nonNullCellsByColumns := 
    {
      VMatrix aux = Tra(Rand(1,_.M,1,1)*Or(_.Q,Zeros(_.M,_.N)));
      Group("ConcatColumns",For(1,_.K,VMatrix(Real k)
      {
        SubRow(aux,Range(k,_.N,_.K)) 
      }))
    };
    WriteLn("TRACE "+_MID+"8");
    VMatrix _.LQ := CholeskiFactor(_.Q,"XtX");
    WriteLn("TRACE "+_MID+"9");
    True
  }
};