Changeset 82ef353 in tol


Ignore:
Timestamp:
Oct 26, 2008, 10:30:40 AM (18 years ago)
Author:
Víctor de Buen <vdebuen@…>
Branches:
1941-docker-image-to-test, devel, master
Children:
de5a7dd
Parents:
6558498
Message:

Using warn_cannot_apply instead of err_cannot_apply in trunc_std_gaussian

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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • tol/btol/vmatrix_type/vmat_stats.cpp

    r6558498 r82ef353  
    537537}
    538538
    539 ////////////////////////////////////////////////////////////////////////////////
    540   int BVMat::trunc_std_gaussian(const BVMat& D,
    541                                 const BVMat& d,
    542                                       BVMat& b,
    543                                       BVMat& z)
    544 // Generates a set of standarized multinormal z constrained to
     539/*
     540////////////////////////////////////////////////////////////////////////////////
     541  int BVMat::trunc_std_kick(const BVMat& D,
     542                            const BVMat& d,
     543                                  BVMat& b,
     544                                  BVMat& z)
     545// Given a vector z matching
    545546//   D*z <= d 
     547// moves it to another point constrained these inequations.
     548// For each component builds bounds conditioning to rest of components
     549//   a[j] <= z[j] <= b[j]
     550// If a[j]>-Inf and b[j]<+Inf then z[j] is generated as a uniform
     551// If a[j]=-Inf and b[j]=+Inf then z[j] is generated as a  N(0,1)
     552// If a[j]=-Inf and b[j]<+Inf then z[j] is generated as a TN(b[j],1,-Inf,b[j])
     553// If a[j]>-Inf and b[j]=+Inf then z[j] is generated as a TN(a[j],1,a[j],+Inf)
    546554// Value od d-D*z is stored in b and must be initialized by the caller function
    547555////////////////////////////////////////////////////////////////////////////////
    … …  
    569577  {
    570578    double& zj = zx[j];
    571   //Std(BText("\nj=")<<j);
     579  //Std(BText("\ntrunc_std_gaussian column j=")<<j);
    572580    lower=negInf;
    573581    upper=posInf;
    … …  
    589597      {
    590598        bx[i] += Dij*zj;
    591       //Std(BText("\nD[")+i+","+j+"]="+Dij+"\tz["+j+"]="+zj+"\tb["+i+"]="+bx[i]+"\td["+i+"]="+dx[i]);
     599      //Std(BText("\n  D[")+i+","+j+"]="+BDat(Dij)+"\tz["+j+"]="+BDat(zj)+"\tb["+i+"]="+BDat(bx[i])+"\td["+i+"]="+BDat(dx[i]));
    592600        if(Dij_sign<0)
    593601        {
    594602          aux = bx[i]/Dij;
    595603          if(lower<aux) { lower = aux; }
    596         //Std(BText("\naux=")+aux+"\tlower="<<lower);
     604        //Std(BText("\n  aux=")+BDat(aux)+"\tlower="<<BDat(lower));
    597605        }
    598606        else //if(Dij_sign>0)
    … …  
    600608          aux = bx[i]/Dij;
    601609          if(upper>aux) { upper = aux; }
    602         //Std(BText("\naux=")+aux+"\tupper="<<upper);
     610        //Std(BText("\n  aux=")+BDat(aux)+"\tupper="<<BDat(upper));
    603611        }
    604612      }
    605613    }
    606   //Std(BText("\nlower=")<<lower<<"\tupper="<<upper);
     614  //Std(BText("\nlower=")<<BDat(lower)<<"\tupper="<<BDat(upper));
    607615    double length = lower-upper;
    608616    int length_sign = gsl_fcmp(length+1.0,1.0,DBL_EPSILON);
    … …  
    615623      else
    616624      {
    617         err_cannot_apply("TruncStdGaussian",
     625        warn_cannot_apply("TruncStdGaussian",
    618626          I2("(empty interval",
    619627          "(intervalo vacio")+" ["+BDat(lower).Format("%.16lg")+","+
    620628                                   BDat(upper).Format("%.16lg")+"])",z);
    621       //return(-1);
     629        return(-1);
    622630      }
    623631    }
    … …  
    662670  return(0);
    663671};
    664 
    665 
    666 ////////////////////////////////////////////////////////////////////////////////
    667   void BVMat::TruncStdGaussian(const BVMat& D, const BVMat& d_,
    668                                const BVMat& z0, int ncol, int burnin)
    669 //Matrix instances
    670 ////////////////////////////////////////////////////////////////////////////////
    671 {
    672   static char* fName = "TruncStdGaussian";
    673   if(!D .CheckDefined(fName)) { return; }
    674   if(!d_.CheckDefined(fName)) { return; }
    675   if(!z0.CheckDefined(fName)) { return; }
    676   int Dr = D.Rows();
    677   int Dc = D.Columns();
    678   int dr = d_.Rows();
    679   int dc = d_.Columns();
    680   int zr = z0.Rows();
    681   int zc = z0.Columns();
    682   if((Dr!=dr)||(dc!=1)||
    683      (Dc!=zr)||(zc!=1)||
    684      !Dc||!Dr||!dc||!dr||!zc||!zr)
    685   {
    686     err_invalid_dimensions(fName,D,d_,z0);
    687     return;
    688   }
    689   double* aux;
    690   int k, nrow = Dc;
    691   BVMat d(d_,ESC_blasRdense);
    692   BVMat z(z0,ESC_blasRdense);
    693   BVMat Dz = D*z;
    694   BVMat b  = d-Dz;
    695   double bMin = b.Min();
    696   int    nnz  = D.NonNullCells();
    697   double drop = nnz*DEpsilon();
    698   if(bMin<-drop)
    699   {
    700     err_cannot_apply(fName,
    701       I2("(non feasible initial point)",
    702          "(valor inicial no factible)"),z0);
    703   //return;
    704   }
    705   if(burnin<=0 ) { burnin = 1; }
    706   BVMat burn;
    707   burn.Zeros(nrow, burnin, ESC_blasRdense);
    708   Zeros(nrow, ncol, ESC_blasRdense);
    709   int size_row = nrow*sizeof(double);
    710   switch(D.code_) {
    711   case(ESC_blasRdense ) :
    712   case(ESC_chlmRsparse ) :
    713     aux = (double*)burn.s_.blasRdense_->x;
    714     for(k=0; k<burnin; k++, aux+=nrow)
    715     {
    716       if(trunc_std_gaussian(D, d, b, z)!=0) { return; } 
    717       memcpy(aux, z.s_.blasRdense_->x, size_row);
    718     }
    719     aux = (double*)s_.blasRdense_->x;
    720     for(k=0; k<ncol; k++, aux+=nrow)
    721     {
    722       if(trunc_std_gaussian(D, d, b, z)!=0) { return; } 
    723       memcpy(aux, z.s_.blasRdense_->x, size_row);
    724     }
    725     break;
    726   default:
    727     err_cannot_create(fName,D.code_);  }
    728 };
    729 
    730 
    731 /*
    732 ////////////////////////////////////////////////////////////////////////////////
    733   int BVMat::constrained_multivariate(const BVMat& D,
    734                                       const BVMat& d,
    735                                             BVMat& b,
    736                                             BVMat& z,
    737                                       const BArray<BTruncatedDist*>& dist)
    738 // Generates a set of independent multivariate z constrained to
     672*/
     673
     674////////////////////////////////////////////////////////////////////////////////
     675  int BVMat::trunc_std_gaussian(const BVMat& D,
     676                                const BVMat& d,
     677                                      BVMat& b,
     678                                      BVMat& z)
     679// Generates a set of standarized multinormal z constrained to
    739680//   D*z <= d 
    740 // given a truncated univariate distribution
    741681// Value od d-D*z is stored in b and must be initialized by the caller function
    742682////////////////////////////////////////////////////////////////////////////////
    … …  
    764704  {
    765705    double& zj = zx[j];
     706  //Std(BText("\ntrunc_std_gaussian column j=")<<j);
     707    lower=negInf;
     708    upper=posInf;
     709    if(sp)
     710    {
     711      k0 = Dp[j];
     712      k1 = Dp[j+1];
     713    }
     714    else
     715    {
     716      hj = h;
     717    }
     718    for(k=k0; k<k1; k++)
     719    {
     720      i = sp?Di[k]:k;
     721      const double& Dij = sp?Dx[k]:Dx[h++];
     722      int Dij_sign = gsl_fcmp(Dij+1.0,1.0,DBL_EPSILON);
     723      if(Dij_sign!=0)
     724      {
     725        bx[i] += Dij*zj;
     726      //Std(BText("\n  D[")+i+","+j+"]="+BDat(Dij)+"\tz["+j+"]="+BDat(zj)+"\tb["+i+"]="+BDat(bx[i])+"\td["+i+"]="+BDat(dx[i]));
     727        if(Dij_sign<0)
     728        {
     729          aux = bx[i]/Dij;
     730          if(lower<aux) { lower = aux; }
     731        //Std(BText("\n  aux=")+BDat(aux)+"\tlower="<<BDat(lower));
     732        }
     733        else //if(Dij_sign>0)
     734        {
     735          aux = bx[i]/Dij;
     736          if(upper>aux) { upper = aux; }
     737        //Std(BText("\n  aux=")+BDat(aux)+"\tupper="<<BDat(upper));
     738        }
     739      }
     740    }
     741  //Std(BText("\nlower=")<<BDat(lower)<<"\tupper="<<BDat(upper));
     742    double length = lower-upper;
     743    int length_sign = gsl_fcmp(length+1.0,1.0,DBL_EPSILON);
     744    if(length_sign>=0)
     745    {
     746      if(length_sign==0)
     747      {
     748        lower = upper;
     749      }
     750      else
     751      {
     752        warn_cannot_apply("TruncStdGaussian",
     753          I2("(empty interval",
     754          "(intervalo vacio")+" ["+BDat(lower).Format("%.16lg")+","+
     755                                   BDat(upper).Format("%.16lg")+"])",z);
     756        return(-1);
     757      }
     758    }
     759# ifdef USE_BTruncatedNormalDist
     760    if(lower>=upper)
     761    {
     762      zj = lower;
     763    }
     764    else if((lower>negInf)||(upper<posInf))
     765    {
     766      BTruncatedNormalDist tn(lower,upper);
     767      zj = tn.Random().Value();
     768      assert((lower<=zj) && (zj<=upper));
     769    //Std(BText("\nBTruncatedNormalDist z[")+j+"]="+zj);
     770    }
     771    else
     772    {
     773      zx[j] = BNormalDist::Random01().Value();
     774    //Std(BText("\nBNormalDist z[")+j+"]="+zj);
     775    }
     776# else
     777    zj = gsl_rtabnorm_combo(BProbDist::rng(), 0.0, 1.0,
     778                     lower, upper, 10);
     779    assert((lower<=zj) && (zj<=upper));
     780    if(!BDat(zj).IsFinite())
     781    {
     782      err_cannot_apply("TruncStdGaussian",
     783        I2("(empty interval",
     784           "(intervalo vacio")+" ["+lower+","+upper+"])",z);
     785      return(-2);
     786    }
     787# endif
     788    if(!sp) { h=hj; }
     789    for(k=k0; k<k1; k++)
     790    {
     791      i = sp?Di[k]:k;
     792      const double& Dij = sp?Dx[k]:Dx[h++];
     793      bx[i] -= Dij*zj;
     794    //Std(BText("\nD[")+i+","+j+"]="+Dij+"\tz["+j+"]="+zj+"\tb["+i+"]="+bx[i]);
     795    }
     796  }
     797  return(0);
     798};
     799
     800
     801////////////////////////////////////////////////////////////////////////////////
     802  void BVMat::TruncStdGaussian(const BVMat& D, const BVMat& d_,
     803                               const BVMat& z0, int ncol, int burnin)
     804//Matrix instances
     805////////////////////////////////////////////////////////////////////////////////
     806{
     807  static char* fName = "TruncStdGaussian";
     808  if(!D .CheckDefined(fName)) { return; }
     809  if(!d_.CheckDefined(fName)) { return; }
     810  if(!z0.CheckDefined(fName)) { return; }
     811  int Dr = D.Rows();
     812  int Dc = D.Columns();
     813  int dr = d_.Rows();
     814  int dc = d_.Columns();
     815  int zr = z0.Rows();
     816  int zc = z0.Columns();
     817  if((Dr!=dr)||(dc!=1)||
     818     (Dc!=zr)||(zc!=1)||
     819     !Dc||!Dr||!dc||!dr||!zc||!zr)
     820  {
     821    err_invalid_dimensions(fName,D,d_,z0);
     822    return;
     823  }
     824  double* aux;
     825  int k, nrow = Dc;
     826  BVMat d(d_,ESC_blasRdense);
     827  BVMat z(z0,ESC_blasRdense);
     828  BVMat Dz = D*z;
     829  BVMat b  = d-Dz;
     830  double bMin = b.Min();
     831  int    nnz  = D.NonNullCells();
     832  double drop = nnz*DEpsilon();
     833  if(bMin<-drop)
     834  {
     835    err_cannot_apply(fName,
     836      I2("(non feasible initial point)",
     837         "(valor inicial no factible)"),z0);
     838    return;
     839  }
     840  if(burnin<=0 ) { burnin = 1; }
     841  BVMat burn;
     842  burn.Zeros(nrow, burnin, ESC_blasRdense);
     843  Zeros(nrow, ncol, ESC_blasRdense);
     844  int size_row = nrow*sizeof(double);
     845  switch(D.code_) {
     846  case(ESC_blasRdense ) :
     847  case(ESC_chlmRsparse ) :
     848    aux = (double*)burn.s_.blasRdense_->x;
     849    for(k=0; k<burnin; k++, aux+=nrow)
     850    {
     851      if(trunc_std_gaussian(D, d, b, z)!=0) { return; } 
     852      memcpy(aux, z.s_.blasRdense_->x, size_row);
     853    }
     854    aux = (double*)s_.blasRdense_->x;
     855    for(k=0; k<ncol; k++, aux+=nrow)
     856    {
     857      if(trunc_std_gaussian(D, d, b, z)!=0) { return; } 
     858      memcpy(aux, z.s_.blasRdense_->x, size_row);
     859    }
     860    break;
     861  default:
     862    err_cannot_create(fName,D.code_);  }
     863};
     864
     865
     866/*
     867////////////////////////////////////////////////////////////////////////////////
     868  int BVMat::constrained_multivariate(const BVMat& D,
     869                                      const BVMat& d,
     870                                            BVMat& b,
     871                                            BVMat& z,
     872                                      const BArray<BTruncatedDist*>& dist)
     873// Generates a set of independent multivariate z constrained to
     874//   D*z <= d 
     875// given a truncated univariate distribution
     876// Value od d-D*z is stored in b and must be initialized by the caller function
     877////////////////////////////////////////////////////////////////////////////////
     878{
     879  int i, j, k, h, hj;
     880  double  aux;
     881  bool    sp = D.code_ == ESC_chlmRsparse;
     882  int     n  = z.s_.blasRdense_->nrow;
     883  int     r  = d.s_.blasRdense_->nrow;
     884  int     k0 = 0;
     885  int     k1 = r;
     886  double* zx = (double*)z.s_.blasRdense_->x;
     887  double* dx = (double*)d.s_.blasRdense_->x;
     888  double* bx = (double*)b.s_.blasRdense_->x;
     889  double* Dx = (double*)(sp?D.s_.chlmRsparse_->x:D.s_.blasRdense_->x);
     890  int*    Di = (int*   )(sp?D.s_.chlmRsparse_->i:NULL);
     891  int*    Dp = (int*   )(sp?D.s_.chlmRsparse_->p:NULL);
     892  double  negInf = BDat::NegInf();
     893  double  posInf = BDat::PosInf();
     894  double  lower  = BDat::NegInf();
     895  double  upper  = BDat::PosInf();
     896  int     nnz    = D.NonNullCells();
     897  double  drop   = nnz*DEpsilon();
     898  for(h=hj=j=0; j<n; j++)
     899  {
     900    double& zj = zx[j];
    766901  //Std(BText("\nj=")<<j);
    767902    lower=negInf;
Note: See TracChangeset for help on using the changeset viewer.