Index: tol_pkg/MhCnstrLinRegArima/test/test_0001/test.tol
===================================================================
--- tol_pkg/MhCnstrLinRegArima/test/test_0001/test.tol	(revision f77a89a8df763d85aec12b8d8bd8c4b823dd5bf8)
+++ tol_pkg/MhCnstrLinRegArima/test/test_0001/test.tol	(revision 77c938e484e61813a8546d9f0f9294fb16f7bd69)
@@ -32,5 +32,6 @@
 Real MLHSteps = 0;
 //Length of MCMC
-Real SampleLength = 500;
+Real Burnin = 500;
+Real SampleLength = 2000;
 
 /////////////////////////////////////////////////////////////////////////////
@@ -38,11 +39,11 @@
 /////////////////////////////////////////////////////////////////////////////
 //Number of variables of each observational node
-Real n = 2;
+Real n = 200;
 //Data length of each observational node
-Real m  = IntRand(100, 400);
+Real m  = n*IntRand(2, 10)+IntRand(100,1000);
 //Number of inequations
-Real r = IntRand(0,2*n);
+Real r = IntRand(0,4*n);
 //Enables/disables full ARIMA components
-Real enableARMA = IntRand(0,1);
+Real enableARMA = IntRand(1,1);
 //Seasonal periodicity
 Real S = 7;
@@ -73,5 +74,5 @@
   {
     Real avr = Rand(-10, 10);
-    Real sig = Rand(2, 4);
+    Real sig = Rand(1, 2);
     Gaussian(m,1,avr,sig)
   }));
@@ -87,5 +88,5 @@
 
 //Standard deviation
-Real s = Rand(0.1,0.4);
+Real s = Rand(0.1,0.2);
 //Variance
 Real s2 = s^2;
@@ -159,9 +160,10 @@
 {
   WriteLn("Initializing Markov Chain ...");
+  Real t0 = Copy(Time);
   Real mclra::_.gInitial := s/2;
   Real mclra::initialize(b0_);
   Matrix mcmc  = Rand(0,mclra::_.n,0,0);
   Matrix cache = Rand(0,mclra::_.n,0,0);
-  While(Rows(mcmc)<SampleLength,
+  While(Rows(mcmc)<Burnin+SampleLength,
   {
     Real If(Rows(mcmc) < MLHSteps,
@@ -177,12 +179,9 @@
     True 
   });
+  Real t1 = Copy(Time);
+  Real ts = (t1-t0)/(Burnin+SampleLength);
+  WriteLn("  Time by simulation:"<<ts+" seconds");
   mcmc
 });
-
-Set Include("test_circle.tol");
-Set Include("test_bsr.tol");
-
-//Chart circle on real target beta and MCMCŽs all together odd pairs against
-Matrix mcmc = circle | bsr.mcmc | Group("ConcatColumns",MCMC);
 
 Set gelman.diag = RCODA::gelman.diag(MCMC,{[[
@@ -192,4 +191,29 @@
   Real transform=FALSE, 
   Real autoburnin=TRUE]]});
+
+Set MCMC.b = EvalSet(MCMC,Matrix(Matrix mcmc)
+{
+  SubRow(mcmc,Range(Burnin+1,Burnin+SampleLength,1))
+});
+Set gelman.diag.b = RCODA::gelman.diag(MCMC.b,{[[
+  Real verbose=TRUE, 
+  Real varByCol=TRUE, 
+  Real confidence=0.95, 
+  Real transform=FALSE, 
+  Real autoburnin=TRUE]]});
+
+Set Include("test_bsr.tol");
+
+//Chart circle on real target beta and MCMCŽs all together odd pairs against
+Matrix mcmc  = Group("ConcatColumns",MCMC);
+Matrix mcmc.b = Matrix If(n==2,
+{
+  Set Include("test_circle.tol");
+  circle | bsr.mcmc | Group("ConcatColumns",MCMC.b)  
+},
+{
+  Set Include("test_circle.tol");
+  circle | bsr.mcmc | Group("ConcatColumns",MCMC.b)
+});
 
 
Index: tol_pkg/MhCnstrLinRegArima/test/test_0001/test_bsr.tol
===================================================================
--- tol_pkg/MhCnstrLinRegArima/test/test_0001/test_bsr.tol	(revision f77a89a8df763d85aec12b8d8bd8c4b823dd5bf8)
+++ tol_pkg/MhCnstrLinRegArima/test/test_0001/test_bsr.tol	(revision 77c938e484e61813a8546d9f0f9294fb16f7bd69)
@@ -55,8 +55,8 @@
 BysMcmc::@Config bsr.config = [[
 //MCMC dimensions
-  Real mcmc.burnin        =  100; 
+  Real mcmc.burnin        =  Burnin; 
   Real mcmc.sampleLength  = SampleLength;
   Real mcmc.cacheLength   =  100;
-  Real bsr.linBlk.MH.useAfterIter   = 10
+  Real bsr.linBlk.MH.useAfterIter   = 1/0
 ]];
 
Index: tol_pkg/MhCnstrLinRegArima/test/test_0001/test_circle.tol
===================================================================
--- tol_pkg/MhCnstrLinRegArima/test/test_0001/test_circle.tol	(revision f77a89a8df763d85aec12b8d8bd8c4b823dd5bf8)
+++ tol_pkg/MhCnstrLinRegArima/test/test_0001/test_circle.tol	(revision 77c938e484e61813a8546d9f0f9294fb16f7bd69)
@@ -8,5 +8,5 @@
 Real mclra::initialize(b);
 Real realMaxLikelihood = mclra::_.log.P;
-Real mclra::_.ratio := s/100;
+Real mclra::_.ratio := s/500;
 Matrix circle = Group("ConcatRows", For(1,SampleLength, Matrix(Real numSim)
 {
Index: tol_pkg/MulVarFunAprox/DesignC2.Monomial.Spheric.MeshFree.tol
===================================================================
--- tol_pkg/MulVarFunAprox/DesignC2.Monomial.Spheric.MeshFree.tol	(revision f77a89a8df763d85aec12b8d8bd8c4b823dd5bf8)
+++ tol_pkg/MulVarFunAprox/DesignC2.Monomial.Spheric.MeshFree.tol	(revision 77c938e484e61813a8546d9f0f9294fb16f7bd69)
@@ -52,5 +52,5 @@
   //////////////////////////////////////////////////////////////////////////////
   {
-    Text _MID = "[@DesignC2.Monomial.Spheric.MeshFree::build.stress.points] ";
+    Text _MID = "[@DesignC2.Monomial.Spheric.Delaunay2D::build.stress.points] ";
     Real _.Z.maxNeighbour := If(prior.Z.order>=0,2,0);
     If(!_.Z.maxNeighbour, 
@@ -74,5 +74,5 @@
     {
       And(
-        LE(distance_[k], Constant(_.J,1,stress.maxDistance)),
+        LE(distance_[k], Constant(_.H,1,stress.maxDistance)),
         LT(Y1_[k],Y2_[k]))
     });
@@ -115,4 +115,5 @@
   })};
 
+
   //////////////////////////////////////////////////////////////////////////////
   Real build.border.points(Real void)
Index: tol_pkg/MulVarFunAprox/DesignC2.Spheric.tol
===================================================================
--- tol_pkg/MulVarFunAprox/DesignC2.Spheric.tol	(revision f77a89a8df763d85aec12b8d8bd8c4b823dd5bf8)
+++ tol_pkg/MulVarFunAprox/DesignC2.Spheric.tol	(revision 77c938e484e61813a8546d9f0f9294fb16f7bd69)
@@ -72,4 +72,68 @@
     SubRow(_.R,VCol2Set(Y.idx))
   };
+
+  //////////////////////////////////////////////////////////////////////////////
+    VMatrix build.obs(VMatrix X)
+  //////////////////////////////////////////////////////////////////////////////
+  {
+    FunRn2R::@BaseC2.Eval.Function eval = 
+       FunRn2R::@BaseC2.Eval.Function::New(get.base(0));
+    MatQuery::@KNN.Neighbourhood ngbh = 
+      _.ann::search(VMat2Mat(X),_.X.maxNeighbour);
+    Real m = VRows(X);
+    VMatrix dist = Mat2VMat(Sqrt(ngbh::_.distance2));
+    VMatrix X.h = Mat2VMat(ngbh::_.neighbour);
+    VMatrix X.idx = DifEq(1/(1-B),Constant(m,1,1));
+    VMatrix triplet.all = Group("ConcatRows", For(1,_.X.maxNeighbour,VMatrix(Real k)
+    {
+      VMatrix X.h.k = SubCol(X.h,[[k]]);
+      VMatrix dist.k = SubCol(dist,[[k]]);
+      VMatrix R = get.radius(X.h.k);
+      VMatrix tao = dist.k $/ R;
+      VMatrix lambda = (-tao^2+1)^4;
+      VMatrix match = LT(dist.k,R);
+      SubRow(X.idx | X.h.k | dist.k | tao| lambda, MatQuery::SelectMatch(match))
+    }));
+    VMatrix x.idx = SubCol(triplet.all,[[1]]);
+    VMatrix y.idx = SubCol(triplet.all,[[2]]);
+    VMatrix xy.idx = x.idx | y.idx;
+    VMatrix lambda.col = SubCol(triplet.all,[[5]]);
+    VMatrix triplet.lambda = xy.idx | lambda.col;
+    VMatrix lambda = Triplet(VMat2Mat(triplet.lambda),m,_.H);
+    VMatrix sum.lambda = lambda*Constant(_.H,1,1);
+    VMatrix sum.lambda.inv = Convert(sum.lambda,"Blas.R.Dense")^-1;
+    VMatrix weight = Eye(m,m,0,sum.lambda.inv)*lambda;
+    VMatrix dif = 
+      SubRow(  X,VCol2Set(SubCol(triplet.all,[[1]]))) -
+      SubRow(_.Y,VCol2Set(SubCol(triplet.all,[[2]]))) ;
+    VMatrix P.row = Group("ConcatRows", For(1,_.K,VMatrix(Real k)
+    {
+      x.idx
+    }));
+    VMatrix P.y = Group("ConcatRows", For(1,_.K,VMatrix(Real k)
+    {
+      y.idx
+    }));
+    VMatrix YnCum = DifEq(1/(1-B),_.Yn);
+    VMatrix YnCumB = Sub(SetCol([[0]])<<YnCum,1,1,_.H,1);
+    VMatrix P.YnCum  = SubRow(YnCum, VCol2Set(P.y));
+    VMatrix P.YnCumB = SubRow(YnCumB,VCol2Set(P.y));
+    VMatrix P.k = Group("ConcatRows", For(1,_.K,VMatrix(Real k)
+    {
+      y.idx*0+k
+    }));
+    VMatrix P.col = P.YnCumB + P.k;
+    VMatrix P.match = LE(P.col,P.YnCum);
+    Real eval::upgrade.cache(dif);
+    VMatrix P.ev = Group("ConcatRows", For(1,_.K,VMatrix(Real k)
+    {
+      lambda.col $* eval::evaluate.v(k,dif)
+    }));
+    Matrix P.triplet = VMat2Mat(
+      SubRow(P.row | P.col | P.ev, MatQuery::SelectMatch(P.match)));
+    VMatrix Px = Pack(Convert(Triplet(P.triplet, m, _.N),"Cholmod.R.Sparse"));
+    VMatrix Qx = Eye(m,m,0,sum.lambda.inv)*Px
+  };
+
 
   ///////////////////////////////////////////////////////////////////////////
Index: tol_pkg/MulVarFunAprox/MulVarFunAprox.tol
===================================================================
--- tol_pkg/MulVarFunAprox/MulVarFunAprox.tol	(revision f77a89a8df763d85aec12b8d8bd8c4b823dd5bf8)
+++ tol_pkg/MulVarFunAprox/MulVarFunAprox.tol	(revision 77c938e484e61813a8546d9f0f9294fb16f7bd69)
@@ -22,4 +22,5 @@
 #Embed "DesignC2.Monomial.Spheric.Grid.tol";
 #Embed "DesignC2.Monomial.Spheric.MeshFree.tol";
+#Embed "DesignC2.Monomial.Spheric.Delaunay2D.tol";
 #Embed "ModelC2.tol"