Changeset 7db6dd2 in tol

Timestamp:

Oct 22, 2008, 2:29:52 PM (18 years ago)

Author:

Víctor de Buen <vdebuen@…>

Branches:

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

Children:

53020ad

Parents:

baa9d61

Message:

New function GetBoundsInPolytope: Calculates minimum and maximum values of j-th component of vector z inside the polytope D*z <= d

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

Location:

tol/btol/vmatrix_type

Files:

• tol_bvmat.h (modified) (1 diff)
• vmat_stats.cpp (modified) (3 diffs)
• vmatgra.cpp (modified) (1 diff)

tol/btol/vmatrix_type/tol_bvmat.h

@@ -498 +498 @@
             ECode c=ESC_blasRdense, int nzmax=-1);
   void Gaussian(int nrow, int ncol, double nu, double sigma);
+  static int GetBoundsInPolytope(const BVMat&  D, 
+                                 const BVMat&  d, 
+                                       BVMat&  z,
+                                       int     j,
+                                       double& lower,
+                                       double& upper);
   void TruncStdGaussian(const BVMat& D, const BVMat& d_, 
                         const BVMat& z0_, int ncol, int burnin);

tol/btol/vmatrix_type/vmat_stats.cpp

@@ -431 +431 @@
 #define USE_BTruncatedNormalDist
 
+
+////////////////////////////////////////////////////////////////////////////////
+  int BVMat::GetBoundsInPolytope(const BVMat&  D, 
+                                 const BVMat&  d_, 
+                                       BVMat&  z_,
+                                       int     j,
+                                       double& lower,
+                                       double& upper)
+// Calculates minimum and maximum values of z[j] matching
+//   D*z <= d  
+////////////////////////////////////////////////////////////////////////////////
+{
+  static char* fName = "GetBoundsInPolytope";
+  if(!D .CheckDefined(fName)) { return(-1); }
+  if(!d_.CheckDefined(fName)) { return(-1); }
+  if(!z_.CheckDefined(fName)) { return(-1); }
+  int Dr = D.Rows();
+  int Dc = D.Columns();
+  int dr = d_.Rows();
+  int dc = d_.Columns();
+  int zr = z_.Rows();
+  int zc = z_.Columns();
+  if((Dr!=dr)||(dc!=1)||
+     (Dc!=zr)||(zc!=1)|| 
+     !Dc||!Dr||!dc||!dr||!zc||!zr)
+  {
+    err_invalid_dimensions(fName,D,d_,z_);
+    return(-1);
+  }
+  BVMat d(d_,ESC_blasRdense);
+  BVMat z(z_,ESC_blasRdense);
+  BVMat  b = d-D*z;
+  int i, k, h;
+  double  aux;
+  bool    sp = D.code_ == ESC_chlmRsparse;
+  int     n  = z.s_.blasRdense_->nrow;
+  int     r  = d.s_.blasRdense_->nrow;
+  int     k0 = 0;
+  int     k1 = r;
+  double* zx = (double*)z.s_.blasRdense_->x; 
+  double* dx = (double*)d.s_.blasRdense_->x;
+  double* bx = (double*)b.s_.blasRdense_->x;
+  double* Dx = (double*)(sp?D.s_.chlmRsparse_->x:D.s_.blasRdense_->x);
+  int*    Di = (int*   )(sp?D.s_.chlmRsparse_->i:NULL);
+  int*    Dp = (int*   )(sp?D.s_.chlmRsparse_->p:NULL);
+  double  negInf = BDat::NegInf();
+  double  posInf = BDat::PosInf();
+  int     nnz    = D.NonNullCells();
+  double  drop   = nnz*DEpsilon();
+  
+  double& zj = zx[j];
+//Std(BText("\nj=")<<j);
+  lower=negInf; 
+  upper=posInf;
+  if(sp)
+  {
+    k0 = Dp[j];
+    k1 = Dp[j+1];
+  }
+  else
+  {
+    h = j*Dr;
+  }
+  for(k=k0; k<k1; k++)
+  {
+    i = sp?Di[k]:k;
+    const double& Dij = sp?Dx[k]:Dx[h++];
+    int Dij_sign = gsl_fcmp(Dij+1.0,1.0,DBL_EPSILON);
+    if(Dij_sign!=0)
+    {
+      bx[i] += Dij*zj;
+    //Std(BText("\nD[")+i+","+j+"]="+Dij+"\tz["+j+"]="+zj+"\tb["+i+"]="+bx[i]+"\td["+i+"]="+dx[i]);
+      if(Dij_sign<0) 
+      {
+        aux = bx[i]/Dij;
+        if(lower<aux) { lower = aux; }
+      //Std(BText("\naux=")+aux+"\tlower="<<lower);
+      }
+      else //if(Dij_sign>0) 
+      {
+        aux = bx[i]/Dij;
+        if(upper>aux) { upper = aux; }
+      //Std(BText("\naux=")+aux+"\tupper="<<upper);
+      }
+    }
+  }
+//Std(BText("\nlower=")<<lower<<"\tupper="<<upper);
+  double length = lower-upper;
+  int length_sign = gsl_fcmp(length+1.0,1.0,DBL_EPSILON);
+  if(length_sign>=0)
+  {
+    if(length_sign==0)
+    {
+      lower = upper;
+    }
+    else
+    {
+      err_cannot_apply("GetBoundsInPolytope",
+        I2("(empty interval",
+           "(intervalo vacio")+" ["+lower+","+upper+"])",z);
+      return(-1);
+    }
+  }
+  return(0);
+}
+
 ////////////////////////////////////////////////////////////////////////////////
   int BVMat::trunc_std_gaussian(const BVMat& D, 
@@ -436 +542 @@
                                       BVMat& b, 
                                       BVMat& z)
-// Generates a set of standarized multibormal z constrained to
+// Generates a set of standarized multinormal z constrained to
 //   D*z <= d  
 // Value od d-D*z is stored in b and must be initialized by the caller function
@@ -563 +669 @@
 ////////////////////////////////////////////////////////////////////////////////
 {
-  if(!D .CheckDefined("TruncStdGaussian")) { return; }
-  if(!d_.CheckDefined("TruncStdGaussian")) { return; }
-  if(!z0.CheckDefined("TruncStdGaussian")) { return; }
   static char* fName = "TruncStdGaussian";
+  if(!D .CheckDefined(fName)) { return; }
+  if(!d_.CheckDefined(fName)) { return; }
+  if(!z0.CheckDefined(fName)) { return; }
   int Dr = D.Rows();
   int Dc = D.Columns();

tol/btol/vmatrix_type/vmatgra.cpp

@@ -2034 +2034 @@
   assert(contens_.Check());
 }
+
+//--------------------------------------------------------------------
+DeclareContensClass(BSet, BSetTemporary, BSetGetBoundsInPolytope);
+DefExtOpr(1, BSetGetBoundsInPolytope, "GetBoundsInPolytope", 4, 4, 
+  "VMatrix VMatrix VMatrix Real",
+  "(VMatrix D, VMatrix d, VMatrix z, Real j)",
+  I2("Calculates minimum and maximum values of j-th component of "
+     "vector z inside the polytope D*z <= d",
+     "Calcula los valores mínimo y máximo de la j-ésima componente "
+     "del vector z dentro del politopo D*z <= d"),
+BOperClassify::MatrixAlgebra_);
+//--------------------------------------------------------------------
+void BSetGetBoundsInPolytope::CalcContens()
+//--------------------------------------------------------------------
+{
+  BVMat& D = VMat(Arg(1));
+  BVMat& d = VMat(Arg(2));
+  BVMat& z = VMat(Arg(3));
+  int j =  (int)Real(Arg(4))-1;
+  double lower;
+  double upper;
+  int ok = BVMat::GetBoundsInPolytope(D, d, z, j, lower, upper)==0;
+  if(ok)
+  {
+    contens_.PrepareStore(2);
+    BGrammar::IncLevel();
+    contens_.AddElement(BContensDat ::New("Lower",lower,""));
+    contens_.AddElement(BContensDat ::New("Upper",upper,""));
+    contens_.SetIndexByName();
+    BGrammar::DecLevel();
+  }
+}
+
 /*
 //--------------------------------------------------------------------