Index: tol/btol/vmatrix_type/vmat_stats.cpp
===================================================================
--- tol/btol/vmatrix_type/vmat_stats.cpp	(revision 205d23b82b035076fef19041a3002946eb2054b9)
+++ tol/btol/vmatrix_type/vmat_stats.cpp	(revision 08b6ebcb6ca9d2c0199edbb985f3cac27d51cfa7)
@@ -585,138 +585,4 @@
 }
 
-/*
-////////////////////////////////////////////////////////////////////////////////
-  int BVMat::trunc_std_kick(const BVMat& D, 
-                            const BVMat& d, 
-                                  BVMat& b, 
-                                  BVMat& z)
-// Given a vector z matching
-//   D*z <= d  
-// moves it to another point constrained these inequations.
-// For each component builds bounds conditioning to rest of components 
-//   a[j] <= z[j] <= b[j]
-// If a[j]>-Inf and b[j]<+Inf then z[j] is generated as a uniform 
-// If a[j]=-Inf and b[j]=+Inf then z[j] is generated as a  N(0,1)
-// If a[j]=-Inf and b[j]<+Inf then z[j] is generated as a TN(b[j],1,-Inf,b[j])
-// If a[j]>-Inf and b[j]=+Inf then z[j] is generated as a TN(a[j],1,a[j],+Inf)
-// 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("\ntrunc_std_gaussian column j=")<<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("\n  D[")+i+","+j+"]="+BDat(Dij)+"\tz["+j+"]="+BDat(zj)+"\tb["+i+"]="+BDat(bx[i])+"\td["+i+"]="+BDat(dx[i]));
-        if(Dij_sign<0) 
-        {
-          aux = bx[i]/Dij;
-          if(lower<aux) { lower = aux; }
-        //Std(BText("\n  aux=")+BDat(aux)+"\tlower="<<BDat(lower));
-        }
-        else //if(Dij_sign>0) 
-        {
-          aux = bx[i]/Dij;
-          if(upper>aux) { upper = aux; }
-        //Std(BText("\n  aux=")+BDat(aux)+"\tupper="<<BDat(upper));
-        }
-      }
-    }
-  //Std(BText("\nlower=")<<BDat(lower)<<"\tupper="<<BDat(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
-      {
-        warn_cannot_apply("TruncStdGaussian",
-          I2("(empty interval",
-          "(intervalo vacio")+" ["+BDat(lower).Format("%.16lg")+","+
-                                   BDat(upper).Format("%.16lg")+"])",z);
-        return(-1);
-      }
-    }
-# ifdef USE_BTruncatedNormalDist
-    if(lower>=upper)
-    {
-      zj = lower;
-    }
-    else if((lower>negInf)||(upper<posInf))
-    {
-      BTruncatedNormalDist tn(lower,upper);
-      zj = tn.Random().Value();
-      assert((lower<=zj) && (zj<=upper));
-    //Std(BText("\nBTruncatedNormalDist z[")+j+"]="+zj);
-    }
-    else
-    {
-      zx[j] = BNormalDist::Random01().Value();
-    //Std(BText("\nBNormalDist z[")+j+"]="+zj);
-    }
-# else
-    zj = gsl_rtabnorm_combo(BProbDist::rng(), 0.0, 1.0,
-                     lower, upper, 10);
-    assert((lower<=zj) && (zj<=upper));
-    if(!BDat(zj).IsFinite())
-    {
-      err_cannot_apply("TruncStdGaussian",
-        I2("(empty interval",
-           "(intervalo vacio")+" ["+lower+","+upper+"])",z);
-      return(-2);
-    }
-# endif
-    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);
-};
-*/
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -913,239 +779,2 @@
 
 
-/*
-////////////////////////////////////////////////////////////////////////////////
-  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*vmat_trunc_eps;
-  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,vmat_trunc_eps);
-      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,vmat_trunc_eps);
-    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*vmat_trunc_eps;
-  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);
-};
-
-/* */