Changeset 08b6ebc in tol


Ignore:
Timestamp:
Sep 21, 2009, 7:59:12 PM (17 years ago)
Author:
Víctor de Buen <vdebuen@…>
Branches:
1941-docker-image-to-test, devel, master
Children:
7154de3
Parents:
205d23b
Message:

Cleaning unused code

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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • tol/btol/vmatrix_type/vmat_stats.cpp

    r205d23b r08b6ebc  
    585585}
    586586
    587 /*
    588 ////////////////////////////////////////////////////////////////////////////////
    589   int BVMat::trunc_std_kick(const BVMat& D,
    590                             const BVMat& d,
    591                                   BVMat& b,
    592                                   BVMat& z)
    593 // Given a vector z matching
    594 //   D*z <= d 
    595 // moves it to another point constrained these inequations.
    596 // For each component builds bounds conditioning to rest of components
    597 //   a[j] <= z[j] <= b[j]
    598 // If a[j]>-Inf and b[j]<+Inf then z[j] is generated as a uniform
    599 // If a[j]=-Inf and b[j]=+Inf then z[j] is generated as a  N(0,1)
    600 // If a[j]=-Inf and b[j]<+Inf then z[j] is generated as a TN(b[j],1,-Inf,b[j])
    601 // If a[j]>-Inf and b[j]=+Inf then z[j] is generated as a TN(a[j],1,a[j],+Inf)
    602 // Value od d-D*z is stored in b and must be initialized by the caller function
    603 ////////////////////////////////////////////////////////////////////////////////
    604 {
    605   int i, j, k, h, hj;
    606   double  aux;
    607   bool    sp = D.code_ == ESC_chlmRsparse;
    608   int     n  = z.s_.blasRdense_->nrow;
    609   int     r  = d.s_.blasRdense_->nrow;
    610   int     k0 = 0;
    611   int     k1 = r;
    612   double* zx = (double*)z.s_.blasRdense_->x;
    613   double* dx = (double*)d.s_.blasRdense_->x;
    614   double* bx = (double*)b.s_.blasRdense_->x;
    615   double* Dx = (double*)(sp?D.s_.chlmRsparse_->x:D.s_.blasRdense_->x);
    616   int*    Di = (int*   )(sp?D.s_.chlmRsparse_->i:NULL);
    617   int*    Dp = (int*   )(sp?D.s_.chlmRsparse_->p:NULL);
    618   double  negInf = BDat::NegInf();
    619   double  posInf = BDat::PosInf();
    620   double  lower  = BDat::NegInf();
    621   double  upper  = BDat::PosInf();
    622   int     nnz    = D.NonNullCells();
    623   double  drop   = nnz*DEpsilon();
    624   for(h=hj=j=0; j<n; j++)
    625   {
    626     double& zj = zx[j];
    627   //Std(BText("\ntrunc_std_gaussian column j=")<<j);
    628     lower=negInf;
    629     upper=posInf;
    630     if(sp)
    631     {
    632       k0 = Dp[j];
    633       k1 = Dp[j+1];
    634     }
    635     else
    636     {
    637       hj = h;
    638     }
    639     for(k=k0; k<k1; k++)
    640     {
    641       i = sp?Di[k]:k;
    642       const double& Dij = sp?Dx[k]:Dx[h++];
    643       int Dij_sign = gsl_fcmp(Dij+1.0,1.0,DBL_EPSILON);
    644       if(Dij_sign!=0)
    645       {
    646         bx[i] += Dij*zj;
    647       //Std(BText("\n  D[")+i+","+j+"]="+BDat(Dij)+"\tz["+j+"]="+BDat(zj)+"\tb["+i+"]="+BDat(bx[i])+"\td["+i+"]="+BDat(dx[i]));
    648         if(Dij_sign<0)
    649         {
    650           aux = bx[i]/Dij;
    651           if(lower<aux) { lower = aux; }
    652         //Std(BText("\n  aux=")+BDat(aux)+"\tlower="<<BDat(lower));
    653         }
    654         else //if(Dij_sign>0)
    655         {
    656           aux = bx[i]/Dij;
    657           if(upper>aux) { upper = aux; }
    658         //Std(BText("\n  aux=")+BDat(aux)+"\tupper="<<BDat(upper));
    659         }
    660       }
    661     }
    662   //Std(BText("\nlower=")<<BDat(lower)<<"\tupper="<<BDat(upper));
    663     double length = lower-upper;
    664     int length_sign = gsl_fcmp(length+1.0,1.0,DBL_EPSILON);
    665     if(length_sign>=0)
    666     {
    667       if(length_sign==0)
    668       {
    669         lower = upper;
    670       }
    671       else
    672       {
    673         warn_cannot_apply("TruncStdGaussian",
    674           I2("(empty interval",
    675           "(intervalo vacio")+" ["+BDat(lower).Format("%.16lg")+","+
    676                                    BDat(upper).Format("%.16lg")+"])",z);
    677         return(-1);
    678       }
    679     }
    680 # ifdef USE_BTruncatedNormalDist
    681     if(lower>=upper)
    682     {
    683       zj = lower;
    684     }
    685     else if((lower>negInf)||(upper<posInf))
    686     {
    687       BTruncatedNormalDist tn(lower,upper);
    688       zj = tn.Random().Value();
    689       assert((lower<=zj) && (zj<=upper));
    690     //Std(BText("\nBTruncatedNormalDist z[")+j+"]="+zj);
    691     }
    692     else
    693     {
    694       zx[j] = BNormalDist::Random01().Value();
    695     //Std(BText("\nBNormalDist z[")+j+"]="+zj);
    696     }
    697 # else
    698     zj = gsl_rtabnorm_combo(BProbDist::rng(), 0.0, 1.0,
    699                      lower, upper, 10);
    700     assert((lower<=zj) && (zj<=upper));
    701     if(!BDat(zj).IsFinite())
    702     {
    703       err_cannot_apply("TruncStdGaussian",
    704         I2("(empty interval",
    705            "(intervalo vacio")+" ["+lower+","+upper+"])",z);
    706       return(-2);
    707     }
    708 # endif
    709     if(!sp) { h=hj; }
    710     for(k=k0; k<k1; k++)
    711     {
    712       i = sp?Di[k]:k;
    713       const double& Dij = sp?Dx[k]:Dx[h++];
    714       bx[i] -= Dij*zj;
    715     //Std(BText("\nD[")+i+","+j+"]="+Dij+"\tz["+j+"]="+zj+"\tb["+i+"]="+bx[i]);
    716     }
    717   }
    718   return(0);
    719 };
    720 */
    721587
    722588////////////////////////////////////////////////////////////////////////////////
    … …  
    913779
    914780
    915 /*
    916 ////////////////////////////////////////////////////////////////////////////////
    917   int BVMat::constrained_multivariate(const BVMat& D,
    918                                       const BVMat& d,
    919                                             BVMat& b,
    920                                             BVMat& z,
    921                                       const BArray<BTruncatedDist*>& dist)
    922 // Generates a set of independent multivariate z constrained to
    923 //   D*z <= d 
    924 // given a truncated univariate distribution
    925 // Value od d-D*z is stored in b and must be initialized by the caller function
    926 ////////////////////////////////////////////////////////////////////////////////
    927 {
    928   int i, j, k, h, hj;
    929   double  aux;
    930   bool    sp = D.code_ == ESC_chlmRsparse;
    931   int     n  = z.s_.blasRdense_->nrow;
    932   int     r  = d.s_.blasRdense_->nrow;
    933   int     k0 = 0;
    934   int     k1 = r;
    935   double* zx = (double*)z.s_.blasRdense_->x;
    936   double* dx = (double*)d.s_.blasRdense_->x;
    937   double* bx = (double*)b.s_.blasRdense_->x;
    938   double* Dx = (double*)(sp?D.s_.chlmRsparse_->x:D.s_.blasRdense_->x);
    939   int*    Di = (int*   )(sp?D.s_.chlmRsparse_->i:NULL);
    940   int*    Dp = (int*   )(sp?D.s_.chlmRsparse_->p:NULL);
    941   double  negInf = BDat::NegInf();
    942   double  posInf = BDat::PosInf();
    943   double  lower  = BDat::NegInf();
    944   double  upper  = BDat::PosInf();
    945   int     nnz    = D.NonNullCells();
    946   double  drop   = nnz*vmat_trunc_eps;
    947   for(h=hj=j=0; j<n; j++)
    948   {
    949     double& zj = zx[j];
    950   //Std(BText("\nj=")<<j);
    951     lower=negInf;
    952     upper=posInf;
    953     if(sp)
    954     {
    955       k0 = Dp[j];
    956       k1 = Dp[j+1];
    957     }
    958     else
    959     {
    960       hj = h;
    961     }
    962     for(k=k0; k<k1; k++)
    963     {
    964       i = sp?Di[k]:k;
    965       const double& Dij = sp?Dx[k]:Dx[h++];
    966       int Dij_sign = gsl_fcmp(Dij+1.0,1.0,vmat_trunc_eps);
    967       if(Dij_sign!=0)
    968       {
    969         bx[i] += Dij*zj;
    970       //Std(BText("\nD[")+i+","+j+"]="+Dij+"\tz["+j+"]="+zj+"\tb["+i+"]="+bx[i]+"\td["+i+"]="+dx[i]);
    971         if(Dij_sign<0)
    972         {
    973           aux = bx[i]/Dij;
    974           if(lower<aux) { lower = aux; }
    975         //Std(BText("\naux=")+aux+"\tlower="<<lower);
    976         }
    977         else //if(Dij_sign>0)
    978         {
    979           aux = bx[i]/Dij;
    980           if(upper>aux) { upper = aux; }
    981         //Std(BText("\naux=")+aux+"\tupper="<<upper);
    982         }
    983       }
    984     }
    985   //Std(BText("\nlower=")<<lower<<"\tupper="<<upper);
    986     double length = lower-upper;
    987     int length_sign = gsl_fcmp(length+1.0,1.0,vmat_trunc_eps);
    988     if(length_sign>=0)
    989     {
    990       if(length_sign==0)
    991       {
    992         lower = upper;
    993       }
    994       else
    995       {
    996         err_cannot_apply("TruncStdGaussian",
    997           I2("(empty interval",
    998              "(intervalo vacio")+" ["+lower+","+upper+"])",z);
    999       //return(-1);
    1000       }
    1001     }
    1002   //if((j>=2)&&(lower<0)&&(upper>0)) { printf(""); }
    1003     dist[j]->SetBounds(lower, upper);
    1004     zj = dist[j]->Random().Value();
    1005     Std(BText("\nTRACE constrained_multivariate ")+lower+"<=z["+j+"]="+zj+"<="+upper);
    1006     assert((lower<=zj) && (zj<=upper));
    1007     if(!sp) { h=hj; }
    1008     for(k=k0; k<k1; k++)
    1009     {
    1010       i = sp?Di[k]:k;
    1011       const double& Dij = sp?Dx[k]:Dx[h++];
    1012       bx[i] -= Dij*zj;
    1013     //Std(BText("\nD[")+i+","+j+"]="+Dij+"\tz["+j+"]="+zj+"\tb["+i+"]="+bx[i]);
    1014     }
    1015   }
    1016   return(0);
    1017 };
    1018 
    1019 
    1020 ////////////////////////////////////////////////////////////////////////////////
    1021   int BVMat::FindFeasiblePoint(const BVMat& D, const BVMat& d_, BVMat& z0)
    1022 ////////////////////////////////////////////////////////////////////////////////
    1023 {
    1024   z0.Delete();
    1025   static char* fName = "FindFeasiblePoint";
    1026   if(!D .CheckDefined(fName)) { return(-1); }
    1027   if(!d_.CheckDefined(fName)) { return(-2); }
    1028   BVMat d(d_,ESC_blasRdense);
    1029   int Dr = D.Rows();
    1030   int Dc = D.Columns();
    1031   int dr = d.Rows();
    1032   int dc = d.Columns();
    1033   if((Dr!=dr)||(dc!=1)||
    1034      !Dc||!Dr||!dc||!dr)
    1035   {
    1036     err_invalid_dimensions(fName,D,d_,z0);
    1037     return(-3);
    1038   }
    1039   BVMat Ir, DI, u, b;
    1040   BArray<BTruncatedDist*> dist(Dc+Dr);
    1041   int i,j,k,iter;
    1042   int cm = 0;
    1043   bool isFeasible = false;
    1044   switch(D.code_) {
    1045   case(ESC_blasRdense ) :
    1046   case(ESC_chlmRsparse ) :
    1047     z0.Zeros(Dc,1,ESC_blasRdense);
    1048     Ir.Eye(Dr,Dr);
    1049     DI = D | Ir;
    1050     u = z0 << d;
    1051   //for(k=0; k<Dc; k++) { dist[k] = new BClosestInsideDist(-1,1,0); }
    1052   //for(k=Dc; k<Dc+Dr; k++) { dist[k] = new BClosestInsideDist(-1,1,0); }
    1053     for(k=0;  k<Dc;    k++) { dist[k] = new BTruncatedNormalDist(-1,1,0,1); }
    1054     for(k=Dc; k<Dc+Dr; k++) { dist[k] = new BTruncatedNormalDist(-1,1,0,1); }
    1055     b = d-DI*u;
    1056     for(iter=0; !isFeasible && (iter<10*(Dr+Dc)); iter++)
    1057     {
    1058       Std(BText("\nTRACE constrained_multivariate iteration ")+iter);
    1059       cm = constrained_multivariate(DI,d,b,u,dist);
    1060       if(cm!=0)
    1061       { return(-4); }
    1062       isFeasible = true;
    1063       k=Dc;
    1064       for(i=0; isFeasible && (i<Dr); i++, k++)
    1065       {
    1066         BDat uk = u.GetCell(k,0);
    1067         isFeasible = (uk>=0.0);
    1068       }
    1069       for(k=Dc; k<Dc+Dr; k++)
    1070       {
    1071         BDat s = BDat(0.5)^iter;
    1072         delete dist[k];
    1073         dist[k] = new BTruncatedNormalDist(-1,1,0,s);
    1074       }
    1075     }
    1076     for(k=0; k<dist.Size(); k++) { delete dist[k]; }
    1077     break;
    1078   default:
    1079     err_cannot_create(fName,D.code_);  return(-5); }
    1080   if(!isFeasible)
    1081   { return(-6); }
    1082   u.Sub(0, 0, Dc, 1, z0);
    1083   return(0);
    1084 }
    1085 
    1086 ////////////////////////////////////////////////////////////////////////////////
    1087   int BVMat::TruncStdGaussian(const BVMat& D, const BVMat& d_,
    1088                                const BVMat& z0, int ncol, int burnin)
    1089 //Matrix instances
    1090 ////////////////////////////////////////////////////////////////////////////////
    1091 {
    1092   static char* fName = "TruncStdGaussian";
    1093   if(!D .CheckDefined(fName)) { return(-1); }
    1094   if(!d_.CheckDefined(fName)) { return(-2); }
    1095   if(!z0.CheckDefined(fName)) { return(-3); }
    1096   int Dr = D.Rows();
    1097   int Dc = D.Columns();
    1098   int dr = d_.Rows();
    1099   int dc = d_.Columns();
    1100   int zr = z0.Rows();
    1101   int zc = z0.Columns();
    1102   if((Dr!=dr)||(dc!=1)||
    1103      (Dc!=zr)||(zc!=1)||
    1104      !Dc||!Dr||!dc||!dr||!zc||!zr)
    1105   {
    1106     err_invalid_dimensions(fName,D,d_,z0);
    1107     return(-4);
    1108   }
    1109   double* aux;
    1110   int k, nrow = Dc;
    1111   BVMat d(d_,ESC_blasRdense);
    1112   BVMat z(z0,ESC_blasRdense);
    1113   BVMat Dz = D*z;
    1114   BVMat b  = d-Dz;
    1115   double bMin = b.Min();
    1116   int    nnz  = D.NonNullCells();
    1117   double drop = nnz*vmat_trunc_eps;
    1118   if(bMin<-drop)
    1119   {
    1120     err_cannot_apply(fName,
    1121       I2("(non feasible initial point)",
    1122          "(valor inicial no factible)"),z0);
    1123     return(-5);
    1124   }
    1125   if(burnin<=0 ) { burnin = 1; }
    1126   BVMat burn;
    1127   burn.Zeros(nrow, burnin, ESC_blasRdense);
    1128   Zeros(nrow, ncol, ESC_blasRdense);
    1129   int size_row = nrow*sizeof(double);
    1130   switch(D.code_) {
    1131   case(ESC_blasRdense ) :
    1132   case(ESC_chlmRsparse ) :
    1133     aux = (double*)burn.s_.blasRdense_->x;
    1134     for(k=0; k<burnin; k++, aux+=nrow)
    1135     {
    1136       if(trunc_std_gaussian(D, d, b, z)!=0) { return(-6); } 
    1137       memcpy(aux, z.s_.blasRdense_->x, size_row);
    1138     }
    1139     aux = (double*)s_.blasRdense_->x;
    1140     for(k=0; k<ncol; k++, aux+=nrow)
    1141     {
    1142       if(trunc_std_gaussian(D, d, b, z)!=0) { return(-7); } 
    1143       memcpy(aux, z.s_.blasRdense_->x, size_row);
    1144     }
    1145     break;
    1146   default:
    1147     err_cannot_create(fName,D.code_);  return(-8); }
    1148   return(0);
    1149 };
    1150 
    1151 /* */
Note: See TracChangeset for help on using the changeset viewer.