Index: tol_pkg/MulVarFunAprox/test/test_0001/test.tol
===================================================================
--- tol_pkg/MulVarFunAprox/test/test_0001/test.tol	(revision 29c08825cd97ebe7b20dcda2d1ecacd993735ecd)
+++ tol_pkg/MulVarFunAprox/test/test_0001/test.tol	(revision 29c08825cd97ebe7b20dcda2d1ecacd993735ecd)
@@ -0,0 +1,480 @@
+/* -*- mode: C++ -*- */
+//////////////////////////////////////////////////////////////////////////////
+// FILE   : test.tol
+// PURPOSE: Package MulVarFunAprox test
+//////////////////////////////////////////////////////////////////////////////
+
+Text email_ = "vdebuen@tol-project.org"; 
+Text link_ = ""; 
+Text summary_ = "Test for mutivariant function aproximation";
+
+//Write here the initial test TOL code if needed
+Real numErr0 = Copy(NError);
+Real numWar0 = Copy(NWarning);
+
+//Random seed setting
+Real PutRandomSeed(0);
+Real rndSeed = GetRandomSeed(0);
+WriteLn("Current Random Seed = "<<rndSeed);
+
+//#Require MulVarFunAprox;
+#Embed "../../../_embed_all.tol";
+
+Real n = 2;
+
+
+Set Base.VFun = {SetOfCode(
+  VMatrix p00 (VMatrix x) { Rand(VRows(x),1,1,1) },
+  VMatrix p10 (VMatrix x) { SubCol(x,[[1]]) },
+  VMatrix p01 (VMatrix x) { SubCol(x,[[2]]) },
+  VMatrix p20 (VMatrix x) { RPow(SubCol(x,[[1]]),2) },
+  VMatrix p02 (VMatrix x) { RPow(SubCol(x,[[2]]),2) },
+  VMatrix p11 (VMatrix x) { SubCol(x,[[1]])$*SubCol(x,[[2]]) }
+)};
+
+Set Base.VDerFst = {[[
+{SetOfCode(
+  VMatrix p00 (VMatrix x) { Rand(VRows(x),1,0,0) },
+  VMatrix p10 (VMatrix x) { Rand(VRows(x),1,1,1) },
+  VMatrix p01 (VMatrix x) { Rand(VRows(x),1,0,0) },
+  VMatrix p20 (VMatrix x) { SubCol(x,[[1]])*2 },
+  VMatrix p02 (VMatrix x) { Rand(VRows(x),1,0,0) },
+  VMatrix p11 (VMatrix x) { SubCol(x,[[2]]) }
+)};
+{SetOfCode(
+  VMatrix p00 (VMatrix x) { Rand(VRows(x),1,0,0) },
+  VMatrix p10 (VMatrix x) { Rand(VRows(x),1,0,0) },
+  VMatrix p01 (VMatrix x) { Rand(VRows(x),1,1,1) },
+  VMatrix p20 (VMatrix x) { Rand(VRows(x),1,0,0) },
+  VMatrix p02 (VMatrix x) { SubCol(x,[[2]])*2 },
+  VMatrix p11 (VMatrix x) { SubCol(x,[[1]]) }
+)}
+]]};
+
+Set Base.VDerSnd = {[[
+{[[
+  {SetOfCode(
+    VMatrix p00 (VMatrix x) { Rand(VRows(x),1,0,0) },
+    VMatrix p10 (VMatrix x) { Rand(VRows(x),1,0,0) },
+    VMatrix p01 (VMatrix x) { Rand(VRows(x),1,0,0) },
+    VMatrix p20 (VMatrix x) { Rand(VRows(x),1,2,2) },
+    VMatrix p02 (VMatrix x) { Rand(VRows(x),1,0,0) },
+    VMatrix p11 (VMatrix x) { Rand(VRows(x),1,0,0) }
+  )};
+  {SetOfCode(
+    VMatrix p00 (VMatrix x) { Rand(VRows(x),1,0,0) },
+    VMatrix p10 (VMatrix x) { Rand(VRows(x),1,0,0) },
+    VMatrix p01 (VMatrix x) { Rand(VRows(x),1,0,0) },
+    VMatrix p20 (VMatrix x) { Rand(VRows(x),1,0,0) },
+    VMatrix p02 (VMatrix x) { Rand(VRows(x),1,0,0) },
+    VMatrix p11 (VMatrix x) { Rand(VRows(x),1,1,1) }
+  )}
+]]},
+{[[
+  {SetOfCode(
+    VMatrix p00 (VMatrix x) { Rand(VRows(x),1,0,0) },
+    VMatrix p10 (VMatrix x) { Rand(VRows(x),1,0,0) },
+    VMatrix p01 (VMatrix x) { Rand(VRows(x),1,0,0) },
+    VMatrix p20 (VMatrix x) { Rand(VRows(x),1,0,0) },
+    VMatrix p02 (VMatrix x) { Rand(VRows(x),1,0,0) },
+    VMatrix p11 (VMatrix x) { Rand(VRows(x),1,1,1) }
+  )};
+  {SetOfCode(
+    VMatrix p00 (VMatrix x) { Rand(VRows(x),1,0,0) },
+    VMatrix p10 (VMatrix x) { Rand(VRows(x),1,0,0) },
+    VMatrix p01 (VMatrix x) { Rand(VRows(x),1,0,0) },
+    VMatrix p20 (VMatrix x) { Rand(VRows(x),1,0,0) },
+    VMatrix p02 (VMatrix x) { Rand(VRows(x),1,2,2) },
+    VMatrix p11 (VMatrix x) { Rand(VRows(x),1,0,0) }
+  )}
+]]} ]]};
+
+Real K = Card(Base.VFun);
+
+/* */
+
+Real J = K*10000;
+
+VMatrix X = Group("ConcatColumns",For(1,n,VMatrix(Real i)
+{
+  Real min = Rand(0,0);
+  Real len = Rand(5,10)*(i^(1/n));
+  Rand(J,1,min,min+len)
+}));
+
+
+VMatrix fun.eval.intra.base(VMatrix X) 
+{
+  Set x = For(1,n,VMatrix(Real i) { SubCol(X,[[i]]) });
+  (x[1]*0+1)*Rand(-10,10)+
+  x[1]*Rand(-10,10)+
+  x[2]*Rand(-10,10)+
+  x[1]^2*Rand(-10,10)+
+  x[2]^2*Rand(-10,10)+
+  x[1]$*x[2]*Rand(-10,10)
+};
+
+VMatrix fun.eval.extra.base(VMatrix X) 
+{
+  Set x = For(1,n,VMatrix(Real i) { SubCol(X,[[i]]) });
+  Sin(x[1])+Cos(x[1])
+};
+
+VMatrix Fx = fun.eval.extra.base(X);
+
+Real lengthRatio = 1*(J/(K*n))^(1/n);
+
+Set X.col.info = For(1,n,Set(Real i)
+{[[
+  VMatrix col = SubCol(X,[[i]]);
+  Real min = VMatMin(col);
+  Real max = VMatMax(col);
+  Real length = max-min;
+  Real grid.size = 1*length/lengthRatio
+]]});
+
+Real grid.size = SetAvr(Traspose(Extract(X.col.info,5))[1]);
+Real radius = grid.size / Sqrt(2);
+
+Set range = For(1,n,Set(Real i)
+{
+  Set info = X.col.info[i];
+  Real points = 1+Ceil(info::length/grid.size);
+  Real exceed = (points-1)*grid.size-info::length;
+  Matrix rnk = Tra(DifEq(1/(1-B),Rand(points,1,1,1))-Rand(points,1,1,1));   
+  Set MatSet((rnk*grid.size)+info::min-exceed/2)[1]
+});
+
+VMatrix Y = Mat2VMat(SetMat(Group("CartProd",range)));
+Real H = VRows(Y);
+Real N = H*K;
+
+VMatrix Y.min = Group("ConcatColumns", For(1,n,VMatrix(Real i)
+{
+  Real aux = range[i][1];
+  Rand(H,1,aux,aux)
+}));
+
+VMatrix Y.max = Group("ConcatColumns", For(1,n,VMatrix(Real i)
+{
+  Real aux = range[i][Card(range[i])];
+  Rand(H,1,aux,aux)
+}));
+
+
+VMatrix Y.index = Mat2VMat(SetMat(Group("CartProd",For(1,n,Set(Real i)
+{
+  Range(1,Card(range[i]),1)
+}))));
+
+VMatrix show.X.Y = 
+  ((X*0) <<  Y   ) |
+  ( X    << (Y*0)); 
+
+VMatrix closest.Y.index(VMatrix A)
+{
+  Real m = VRows(A);
+  VMatrix cumulated = Rand(m,1,1,1);
+  VMatrix A.position = Rand(m,1,1,1)+SetSum(For(1,n,VMatrix(Real i_)
+  {
+    Real i = n-i_+1;
+    Real cRi = Card(range[i]);
+    Set info = X.col.info[i];
+    VMatrix Ai = SubCol(A,[[i]]);
+    VMatrix Yi.min = Rand(m,1,range[i][1],range[i][1]);
+    VMatrix Ri = Round((Ai-Yi.min)/grid.size);
+    VMatrix minRi = Rand(m,1,0,0);
+    VMatrix maxRi = Rand(m,1,cRi-1,cRi-1);
+    VMatrix aux = Min(maxRi,Max(minRi,Ri));
+    VMatrix pos = aux $* cumulated;
+    VMatrix cumulated := cumulated*cRi;
+    pos
+  }))
+};
+
+VMatrix Y.closest.Y.index = closest.Y.index(Y);
+VMatrix Y.closest.Y = SubRow(Y,VCol2Set(Y.closest.Y.index));
+Real ok.1 = VMatMax(Abs(Y-Y.closest.Y))==0;
+
+VMatrix Y.neighbour.pos = Min(Y.max, Y + Rand(H,n,grid.size,grid.size));
+VMatrix Y.neighbour.neg = Max(Y.min, Y - Rand(H,n,grid.size,grid.size));
+
+VMatrix Y.neighbourhood = DifEq(1/(1-B),Rand(H,1,1,1)) |
+{
+  Set aux.1 = BinGroup("<<", For(1,n,Set(Real i)
+  {
+    VMatrix pre  = SubCol(Y,Range(1,  i-1,1));
+    VMatrix post = SubCol(Y,Range(i+1,n,  1));
+    SetOfVMatrix(pre | SubCol(Y.neighbour.neg,[[i]]) | post,
+                 pre | SubCol(Y.neighbour.pos,[[i]]) | post)
+  }));
+  Group("ConcatColumns",EvalSet(aux.1,closest.Y.index))
+};
+
+Set Y.neighbourhood.extrictly = For(1,H,Set(Real h) 
+{ 
+  VCol2Set(SubRow(Y.neighbourhood,[[h]]))  - [[h]]
+});
+
+Real X.maxNeighbour = 1+2*n;
+
+VMatrix X.closest.Y.index = closest.Y.index(X);
+VMatrix X.closest.Y = SubRow(Y,VCol2Set(X.closest.Y.index));
+
+VMatrix X.j=DifEq(1/(1-B),Rand(J,1,1,1));
+
+VMatrix X.h=SubRow(Y.neighbourhood, VCol2Set(X.closest.Y.index));
+
+Set X.dif = For(1,X.maxNeighbour,VMatrix(Real i)
+{
+  X - SubRow(Y,VCol2Set(SubCol(X.h,[[i]])))
+});
+
+Set X.tao = For(1,X.maxNeighbour,VMatrix(Real i)
+{
+  VMatrix X.norm2 = RPow(X.dif[i],2)*Rand(n,1,1,1);
+  RPow(X.norm2,0.5)/radius
+});
+
+Set X.tao2 = For(1,X.maxNeighbour,VMatrix(Real i)
+{
+  RPow(X.tao[i],2)
+});
+
+Set X.match = For(1,X.maxNeighbour,VMatrix(Real i)
+{
+  LE(X.tao[i],Rand(J,1,1,1))
+});
+
+Set X.lambda = For(1,X.maxNeighbour,VMatrix(Real i)
+{
+  RPow(X.match[i] $* (Rand(J,1,1,1)-X.tao2[i]),4)
+});
+
+VMatrix X.lambda.sum = SetSum(X.lambda);
+
+Set X.weight = For(1,X.maxNeighbour,VMatrix(Real i)
+{
+  X.lambda[i] $/ X.lambda.sum
+});
+
+Matrix QX.triplet = Group("ConcatRows", For(1,X.maxNeighbour,Matrix(Real i)
+{
+  VMatrix X.hiK = (SubCol(X.h,[[i]])-1)*K;
+  Group("ConcatRows", For(1,K,Matrix(Real k)
+  {
+    Code bvf = Base.VFun[k];
+    VMatrix pk = bvf(X.dif[i]);
+    VMatrix qk = X.weight[i]$*pk; 
+    Set sel.rows = MatQuery::SelectMatch(NE(Drop(Abs(qk),1.E-15),qk*0));
+    VMat2Mat(
+      SubRow(X.j,sel.rows) |
+     (SubRow(X.hiK,sel.rows)+k) |
+      SubRow(qk,sel.rows))
+  }))
+}));
+
+VMatrix QX = Pack(Convert(Triplet(QX.triplet, J, N),"Cholmod.R.Sparse"));
+
+VMatrix displ.Z.1 = 
+  Rand(H,1,-0.1,0.1)*grid.size/2 |
+  Rand(H,1, 0.9,1.1)*grid.size/2;
+VMatrix displ.Z.2 = 
+  Rand(H,1, 0.9,1.1)*grid.size/2 |
+  Rand(H,1,-0.1,0.1)*grid.size/2;
+VMatrix displ.Z.3 = 
+  Rand(H,1,-0.9,-1.1)*grid.size/2 |
+  Rand(H,1,-0.9,-1.1)*grid.size/2;
+
+VMatrix Z.all = (Y+displ.Z.1) << (Y+displ.Z.2) << (Y+displ.Z.3);
+VMatrix Z.Y.1.all = Y << Y << Y;
+VMatrix Z.Y.2.all = (Y+displ.Z.1*2) << (Y+displ.Z.2*2) << (Y+displ.Z.3*2);
+
+VMatrix Z = Z.all;
+VMatrix Z.Y.1 = Z.Y.1.all;
+VMatrix Z.Y.2 = Z.Y.2.all;
+VMatrix Z.h.1 = closest.Y.index(Z.Y.1);
+VMatrix Z.h.2 = closest.Y.index(Z.Y.2);
+
+
+VMatrix show.X.Y.Z = 
+  ( X    << (Y*0) << (Z*0)) |
+  ((X*0) <<  Y    << (Z*0)) | 
+  ((X*0) << (Y*0) <<  Z   ); 
+
+VMatrix show.Y.Z = 
+  ( Y    << (Z*0)) | 
+  ((Y*0) <<  Z   ); 
+
+
+Real T = VRows(Z);
+Real prior.Z.order = 2;
+Real Mz = Case(
+  prior.Z.order< 0, 0,
+  prior.Z.order==0, T,
+  prior.Z.order==1, T*(1+n*K),
+  prior.Z.order==2, T*(1+(n^2)*K));
+
+Real Z.maxNeighbour = If(Mz,2,0);
+VMatrix Z.j=DifEq(1/(1-B),Rand(T,1,1,1));
+
+VMatrix Z.h=Z.h.1 | Z.h.2;
+
+Set Z.dif = For(1,Z.maxNeighbour,VMatrix(Real i)
+{
+  Z - SubRow(Y,VCol2Set(SubCol(Z.h,[[i]])))
+});
+
+Set Z.tao = For(1,Z.maxNeighbour,VMatrix(Real i)
+{
+  VMatrix Z.norm2 = RPow(Z.dif[i],2)*Rand(n,1,1,1);
+  RPow(Z.norm2,0.5)/radius
+});
+
+Set Z.match = For(1,Z.maxNeighbour,VMatrix(Real i)
+{
+  LE(Z.tao[i],Rand(T,1,1,1))
+});
+
+Matrix PZ.triplet = Group("ConcatRows", For(1,Z.maxNeighbour,Matrix(Real i)
+{
+  VMatrix Z.hiK = (SubCol(Z.h,[[i]])-1)*K;
+  Group("ConcatRows", For(1,K,Matrix(Real k)
+  {
+    Real cumRows = 0;
+    Code bvf = Base.VFun[k];
+    VMatrix pk = bvf(Z.dif[i])*If(i==1,+1,-1);
+    Set sel.rows = MatQuery::SelectMatch(NE(Drop(Abs(pk),1.E-15),pk*0));
+    VMatrix col.idx = (SubRow(Z.hiK,sel.rows)+k);
+    {
+      VMatrix row.idx = SubRow(Z.j,sel.rows)+cumRows;
+      Real cumRows := cumRows + T;
+      VMat2Mat(row.idx | col.idx | SubRow(pk,sel.rows))
+    }  << If(prior.Z.order<1,Rand(0,3,0,0),
+    Group("ConcatRows",For(1,n,Matrix(Real d)
+    {
+      Code bvfd = Base.VDerFst[d][k];
+      VMatrix pkd = bvfd(Z.dif[i])*If(i==1,+1/2,-1/2);
+      VMatrix row.idx = SubRow(Z.j,sel.rows)+cumRows;
+      Real cumRows := cumRows + T;
+      VMat2Mat(row.idx | col.idx | SubRow(pkd,sel.rows))
+    })))<< If(prior.Z.order<2,Rand(0,3,0,0),
+    Group("ConcatRows",For(1,n,Matrix(Real d1)
+    {
+      Group("ConcatRows",For(d1,n,Matrix(Real d2)
+      {
+        Code bvfd2 = Base.VDerSnd[d1][d2][k];
+        VMatrix pkd2 = bvfd2(Z.dif[i])*If(i==1,+1/4,-1/4);
+        VMatrix row.idx = SubRow(Z.j,sel.rows)+cumRows;
+        Real cumRows := cumRows + T;
+        VMat2Mat(row.idx | col.idx | SubRow(pkd2,sel.rows))
+      }))
+    })))
+  }))
+}));
+
+
+Real M = J+Mz;
+
+VMatrix PZ = Pack(Convert(Triplet(PZ.triplet, Mz, N),"Cholmod.R.Sparse"));
+
+VMatrix Fz = Rand(Mz,1,0,0);
+VMatrix F = Fx<<Fz;
+VMatrix Q.all = QX << PZ;
+Set Q.active = MatQuery::SelectMatch(
+  NE(Tra(Rand(1,M,1,1)*Abs(Q.all)),Rand(N,1,0,0)));
+
+VMatrix Q = SubCol(Q.all,Q.active);
+
+VMatrix LQ = CholeskiFactor(Q,"XtX");
+
+VMatrix alpha = CholeskiSolve(LQ,Tra(Tra(F)*Q),"PtLLtP");
+VMatrix F.frc = SubRow(Q*alpha,Range(1,J,1));
+VMatrix F.cmp = Fx | F.frc;
+VMatrix E = Fx - F.frc;
+Real E.avr = VMatAvr(E);
+Real E.ste = VMatStDs(E);
+
+Set View([[
+Real n, K, J, H, 
+X.maxNeighbour, Z.maxNeighbour,
+T, N, Mz,M, E.avr, E.ste, "\n"]],"");
+
+/*
+Matrix Q_ = VMat2Mat(Q);
+
+Set svd = SVD(Q_);
+Matrix svd.D = Tra(SubDiag(svd[2],0));
+
+/*
+Real numErr1 = Copy(NError);
+Real numWar1 = Copy(NWarning);
+
+Set partialResults_ = [[numErr0, numErr1, 
+                        ok.1, ok.2]];
+
+//This is a messure of the success of the test 
+Real quality_ = And(numErr1 == numErr0, 
+                    numWar1 == numWar0, 
+                    ok.1, ok.2);
+
+//Return the results 
+Set resultStr_ = @strTestStatus(summary_, link_, quality_,
+                  "Partial results = "<<partialResults_,
+                  "NO DBServerType", "NO DBServerHost", "NO DBAlias",
+                  email_);
+WriteLn(""<<resultStr_);
+resultStr_;
+
+/* */
+
+/*
+VMatrix Z.min = Group("ConcatColumns", For(1,n,VMatrix(Real i)
+{
+  Real aux = range[i][1];
+  Rand(2*H,1,aux,aux)
+}));
+
+VMatrix Z.max = Group("ConcatColumns", For(1,n,VMatrix(Real i)
+{
+  Real aux = range[i][Card(range[i])];
+  Rand(2*H,1,aux,aux)
+}));
+
+VMatrix Z.interior = EQ(
+  And(GE(Z.all,Z.min), LE(Z.all,Z.max))*Rand(n,1,1,1),
+  Rand(2*H,1,n,n));
+
+VMatrix Z = SubRow(Z.all,MatQuery::SelectMatch(Z.interior));
+VMatrix Z.Y.1 = SubRow(Z.Y.1.all,MatQuery::SelectMatch(Z.interior));
+VMatrix Z.Y.2 = SubRow(Z.Y.2.all,MatQuery::SelectMatch(Z.interior));
+VMatrix Z.h.1 = closest.Y.index(Z.Y.1);
+VMatrix Z.h.2 = closest.Y.index(Z.Y.2);
+*/
+
+/*
+Real deg = 3;
+
+Real dim = 2;
+Real maxDeg = 3;
+VMatrix candidates = Mat2VMat(SetMat(Range(0,maxDeg,1)^dim));
+VMatrix match = LE(candidates*Rand(dim,1,1,1),Rand(VRows(candidates),1,maxDeg,maxDeg));
+VMatrix degrees = SubRow(candidates, MatQuery::SelectMatch(match));
+Real K = VRows(degrees);
+Class @PolynBasic
+{
+  Real k;
+  VMatrix fun (VMatrix x) 
+  {
+    SetProd(For(1,dim,VMatrix(Real i)
+    {
+      
+    })
+  };
+  
+};
+Set base = For(1,K,NameBlock(Real k)
+{
+});
+/* */
+