Changeset 43274223 in tol

Timestamp:

Oct 20, 2008, 8:11:38 AM (18 years ago)

Author:

Víctor de Buen <vdebuen@…>

Branches:

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

Children:

ada82ec

Parents:

f22ae37

Message:

Just commentaries

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

File:

• tol/btol/vmatrix_type/vmat_stats.cpp (modified) (2 diffs)

tol/btol/vmatrix_type/vmat_stats.cpp

@@ -436 +436 @@
                                       BVMat& b, 
                                       BVMat& z)
-// Generates z as a truncated standard multinormal
-// b <= d  
-// b = D*z 
+// Generates a set of standarized multibormal z constrained to
+//   D*z <= d  
+// Value od d-D*z is stored in b and must be initialized by the caller function
 ////////////////////////////////////////////////////////////////////////////////
 {
@@ -621 +621 @@
 };
 
+
+/*
+////////////////////////////////////////////////////////////////////////////////
+  int BVMat::constrained_multivariate(const BVMat& D, 
+                                      const BVMat& d, 
+                                            BVMat& b, 
+                                            BVMat& z,
+                                      const BArray<BTruncatedDist*>& dist)
+// Generates a set of independent multivariate z constrained to
+//   D*z <= d  
+// given a truncated univariate distribution
+// Value od d-D*z is stored in b and must be initialized by the caller function
+////////////////////////////////////////////////////////////////////////////////
+{
+  int i, j, k, h, hj;
+  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();
+  double  lower  = BDat::NegInf();
+  double  upper  = BDat::PosInf();
+  int     nnz    = D.NonNullCells();
+  double  drop   = nnz*DEpsilon();
+  for(h=hj=j=0; j<n; j++)
+  {
+    double& zj = zx[j];
+  //Std(BText("\nj=")<<j);
+    lower=negInf; 
+    upper=posInf;
+    if(sp)
+    {
+      k0 = Dp[j];
+      k1 = Dp[j+1];
+    }
+    else
+    {
+      hj = h;
+    }
+    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("TruncStdGaussian",
+          I2("(empty interval",
+             "(intervalo vacio")+" ["+lower+","+upper+"])",z);
+      //return(-1);
+      }
+    }
+  //if((j>=2)&&(lower<0)&&(upper>0))      printf("");
+    dist[j]->SetBounds(lower, upper);
+    zj = dist[j]->Random().Value();
+    Std(BText("\nTRACE constrained_multivariate ")+lower+"<=z["+j+"]="+zj+"<="+upper);
+    assert((lower<=zj) && (zj<=upper));
+    if(!sp) { h=hj; }
+    for(k=k0; k<k1; k++)
+    {
+      i = sp?Di[k]:k;
+      const double& Dij = sp?Dx[k]:Dx[h++];
+      bx[i] -= Dij*zj;
+    //Std(BText("\nD[")+i+","+j+"]="+Dij+"\tz["+j+"]="+zj+"\tb["+i+"]="+bx[i]);
+    }
+  }
+  return(0);
+};
+
+
+////////////////////////////////////////////////////////////////////////////////
+  int BVMat::FindFeasiblePoint(const BVMat& D, const BVMat& d_, BVMat& z0)
+////////////////////////////////////////////////////////////////////////////////
+{
+  z0.Delete();
+  static char* fName = "FindFeasiblePoint";
+  if(!D .CheckDefined(fName)) { return(-1); }
+  if(!d_.CheckDefined(fName)) { return(-2); }
+  BVMat d(d_,ESC_blasRdense);
+  int Dr = D.Rows();
+  int Dc = D.Columns();
+  int dr = d.Rows();
+  int dc = d.Columns();
+  if((Dr!=dr)||(dc!=1)||
+     !Dc||!Dr||!dc||!dr)
+  {
+    err_invalid_dimensions(fName,D,d_,z0);
+    return(-3);
+  }
+  BVMat Ir, DI, u, b;
+  BArray<BTruncatedDist*> dist(Dc+Dr);
+  int i,j,k,iter;
+  int cm = 0;
+  bool isFeasible = false;
+  switch(D.code_) {
+  case(ESC_blasRdense ) :
+  case(ESC_chlmRsparse ) :
+    z0.Zeros(Dc,1,ESC_blasRdense);
+    Ir.Eye(Dr,Dr);
+    DI = D | Ir;
+    u = z0 << d;
+  //for(k=0; k<Dc; k++) { dist[k] = new BClosestInsideDist(-1,1,0); }
+  //for(k=Dc; k<Dc+Dr; k++) { dist[k] = new BClosestInsideDist(-1,1,0); }
+    for(k=0;  k<Dc;    k++) { dist[k] = new BTruncatedNormalDist(-1,1,0,1); }
+    for(k=Dc; k<Dc+Dr; k++) { dist[k] = new BTruncatedNormalDist(-1,1,0,1); }
+    b = d-DI*u;
+    for(iter=0; !isFeasible && (iter<10*(Dr+Dc)); iter++)
+    {
+      Std(BText("\nTRACE constrained_multivariate iteration ")+iter);
+      cm = constrained_multivariate(DI,d,b,u,dist);
+      if(cm!=0) 
+      { return(-4); }
+      isFeasible = true;
+      k=Dc;
+      for(i=0; isFeasible && (i<Dr); i++, k++)
+      {
+        BDat uk = u.GetCell(k,0);
+        isFeasible = (uk>=0.0);
+      }
+      for(k=Dc; k<Dc+Dr; k++) 
+      { 
+        BDat s = BDat(0.5)^iter;
+        delete dist[k];
+        dist[k] = new BTruncatedNormalDist(-1,1,0,s); 
+      }
+    }
+    for(k=0; k<dist.Size(); k++) { delete dist[k]; }
+    break;
+  default: 
+    err_cannot_create(fName,D.code_);  return(-5); }
+  if(!isFeasible)
+  { return(-6); }
+  u.Sub(0, 0, Dc, 1, z0);
+  return(0);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+  int BVMat::TruncStdGaussian(const BVMat& D, const BVMat& d_, 
+                               const BVMat& z0, int ncol, int burnin)
+//Matrix instances
+////////////////////////////////////////////////////////////////////////////////
+{
+  static char* fName = "TruncStdGaussian";
+  if(!D .CheckDefined(fName)) { return(-1); }
+  if(!d_.CheckDefined(fName)) { return(-2); }
+  if(!z0.CheckDefined(fName)) { return(-3); }
+  int Dr = D.Rows();
+  int Dc = D.Columns();
+  int dr = d_.Rows();
+  int dc = d_.Columns();
+  int zr = z0.Rows();
+  int zc = z0.Columns();
+  if((Dr!=dr)||(dc!=1)||
+     (Dc!=zr)||(zc!=1)|| 
+     !Dc||!Dr||!dc||!dr||!zc||!zr)
+  {
+    err_invalid_dimensions(fName,D,d_,z0);
+    return(-4);
+  }
+  double* aux;
+  int k, nrow = Dc;
+  BVMat d(d_,ESC_blasRdense);
+  BVMat z(z0,ESC_blasRdense);
+  BVMat Dz = D*z;
+  BVMat b  = d-Dz;
+  double bMin = b.Min();
+  int    nnz  = D.NonNullCells();
+  double drop = nnz*DEpsilon();
+  if(bMin<-drop)
+  {
+    err_cannot_apply(fName,
+      I2("(non feasible initial point)",
+         "(valor inicial no factible)"),z0);
+    return(-5);
+  }
+  if(burnin<=0 ) { burnin = 1; }
+  BVMat burn;
+  burn.Zeros(nrow, burnin, ESC_blasRdense);
+  Zeros(nrow, ncol, ESC_blasRdense);
+  int size_row = nrow*sizeof(double);
+  switch(D.code_) {
+  case(ESC_blasRdense ) :
+  case(ESC_chlmRsparse ) :
+    aux = (double*)burn.s_.blasRdense_->x;
+    for(k=0; k<burnin; k++, aux+=nrow)
+    {
+      if(trunc_std_gaussian(D, d, b, z)!=0) { return(-6); }  
+      memcpy(aux, z.s_.blasRdense_->x, size_row);
+    }
+    aux = (double*)s_.blasRdense_->x;
+    for(k=0; k<ncol; k++, aux+=nrow)
+    {
+      if(trunc_std_gaussian(D, d, b, z)!=0) { return(-7); }  
+      memcpy(aux, z.s_.blasRdense_->x, size_row);
+    }
+    break;
+  default: 
+    err_cannot_create(fName,D.code_);  return(-8); }
+  return(0);
+};
+
 /* */