Changeset c209b7c in tol

Timestamp:

Oct 26, 2008, 10:51:11 AM (18 years ago)

Author:

Víctor de Buen <vdebuen@…>

Branches:

1941-docker-image-to-test, devel, master

Children:

eaba07d

Parents:

64e44a5

Message:

Trying to avoid points that do not match constrains

git-svn-id: ​https://www.tol-project.org/svn/tolp/trunk@87 be4846ae-f051-483a-831f-b64e7d08527f

File:

• tol/stdlib/math/stat/models/bayesian/bysMcmc/bsr/gibbs/_stdLinearBlock.tol (modified) (3 diffs)

tol/stdlib/math/stat/models/bayesian/bysMcmc/bsr/gibbs/_stdLinearBlock.tol

@@ -41 +41 @@
   //explicitely stored
   VMatrix _.L  = Rand(0,0,0,0);
+  VMatrix _.Ls = Rand(0,0,0,0);
+  //Last known valid initial point for constrained simulation
+  VMatrix _.beta0 = Rand(0,0,0,0);
   //If you know that design matrix will be the same for all simulations you
   //should to set this member as TRUE to avoid unneeded Cholesky 
@@ -123 +126 @@
         CholeskiFactor(XtX,"X")
       });
+      VMatrix If(_.R, _.Ls := Convert(_.L, "Cholmod.R.Sparse"));
       1
     }) 
@@ -187 +191 @@
   ////////////////////////////////////////////////////////////////////////////
   {
-    VMatrix   Ls  = Convert(_.L, "Cholmod.R.Sparse");
+  //WriteLn("TRACE "+_MID+" draw.constrained 1.1 Max(A*b0-a)="<<VMatMax(_.A*_.blockDef::_.store - a));
     VMatrix   d   = _.a - _.A*sr::nu;
     VMatrix   D   = Tra(CholeskiSolve(_.L,_.At,"PtL"));
-  //WriteLn("TRACE "+_MID+" draw.constrained 2 check D="<<Real VMatMax(Abs(_.At-CholeskiSolve(_.L,Ls*Tra(D),"Pt"))));
-    VMatrix   dif = _.blockDef::_.store-sr::nu;
-    VMatrix   b0  = Tra(Tra(CholeskiSolve(_.L,dif,"P"))*Ls);
-  //WriteLn("TRACE "+_MID+" draw.constrained 3 check b0="<<Real VMatMax(Abs(dif-CholeskiSolve(_.L,b0,"LtP"))));
-    VMatrix   b1  = TruncStdGaussian(D,d,b0,1,
+    VMatrix   z0  = Tra(Tra(_.blockDef::_.store-sr::nu)*_.Ls);
+  //WriteLn("TRACE "+_MID+" draw.constrained 1.2 Max(D*z0-d)="<<VMatMax(D*z0 - d));
+    VMatrix   z   = TruncStdGaussian(D,d,z0,1,
                      _.blockDef::_.config::bsr.truncMNormal.gibbsNumIter);
-    VMatrix   b2  = CholeskiSolve(_.L,b1,"LtP");
-  //WriteLn("TRACE "+_MID+" draw.constrained 4 check b2="<<Real VMatMax(Abs(b0-Tra(Tra(CholeskiSolve(_.L,b1,"P"))*Ls))));
-    VMatrix  (_.blockDef::_.store) := sr::nu+b2;
+    Real distBorder = VMatMax(D*z-d);
+    Real If(distBorder<=0,
+    {
+      VMatrix _.beta0 := _.blockDef::_.store;
+      1
+    },
+    {
+      WriteLn(_MID+" Found a point that doesn't match constrains Max(D*z-d)="<<
+        VMatMax(D*z-d)+"\n Trying again");
+      If(!VRows(_.beta0),0,
+      {
+        VMatrix z0 := Tra(Tra(_.beta0-sr::nu)*_.Ls);
+        VMatrix z  := TruncStdGaussian(D,d,z0,1,
+                     _.blockDef::_.config::bsr.truncMNormal.gibbsNumIter);
+        1
+      });
+      1
+    });
+    VMatrix  (_.blockDef::_.store) := sr::nu + CholeskiSolve(_.L,z,"LtP");
     Copy(_.blockDef::_.store)
   };