sigma2 changed to an array to allow reestimation

PyUQTk/PyPCE/pce_tools.py

@@ -374,7 +374,7 @@ def UQTkRegression(pc_model,f_evaluations, samplepts):
 def UQTkBCS(pc_begin, xdata, ydata, eta=1.e-3, niter=1, mindex_growth=None, ntry=1,\
  eta_folds=5, eta_growth = False, eta_plot = False,\
  regparams=None, sigma2=1e-8, npccut=None, pcf_thr=None,\
- verbose=0):
+ verbose=0, return_sigma2=False):
  """
  Obtain PC coefficients by Bayesian compressive sensing
 
@@ -411,7 +411,8 @@ def UQTkBCS(pc_begin, xdata, ydata, eta=1.e-3, niter=1, mindex_growth=None, ntry
  pcf_thr: Minimum value (magnitude) for PC coefficients, for pruning low PC coefficients 'by hand'
  default is None
  verbose: Flag for optional print statements
- 
+ return_sigma2: Flag to retun reestimated sigma2
+
 
  Output:
  pc_model_final: PC object with basis expanded by the iterations
@@ -457,7 +458,7 @@ def UQTkBCS(pc_begin, xdata, ydata, eta=1.e-3, niter=1, mindex_growth=None, ntry
  full_basis_size = pc_begin.GetNumberPCTerms()
  else:
  full_basis_size = uqtkpce.PCSet("NISPnoq", pc_begin.GetOrder() + niter -1, pc_begin.GetNDim(), pc_begin.GetPCType(), pc_begin.GetAlpha(), pc_begin.GetBeta()).GetNumberPCTerms()
- 
+
  # loop through iterations with different splits of the data
  for i in range(ntry):
  # Split the data
@@ -482,7 +483,7 @@ def UQTkBCS(pc_begin, xdata, ydata, eta=1.e-3, niter=1, mindex_growth=None, ntry
  print(mindex)
 
  # One run of BCS to obtain an array of coefficients and a new multiindex
- c_k, used_mi_np = UQTkEvalBCS(pc_model, y_split, x_split, sigma2, eta_opt, regparams, verbose)
+ c_k, used_mi_np, sigma2 = UQTkEvalBCS(pc_model, y_split, x_split, sigma2, eta_opt, regparams, verbose)
 
  # Custom 'cuts' by number of PC terms or by value of PC coefficients
  npcall = c_k.shape[0] # number of PC terms
@@ -560,8 +561,12 @@ def UQTkBCS(pc_begin, xdata, ydata, eta=1.e-3, niter=1, mindex_growth=None, ntry
  print("Coefficients:")
  print(cfs_final)
  print(len(cfs_final), " terms retained out of a full basis of size", full_basis_size)
+ print("Reestimated sigma2:", sigma2)
 
- return pc_model_final, cfs_final
+ if return_sigma2:
+ return pc_model_final, cfs_final, sigma2
+ else:
+ return pc_model_final, cfs_final
 ################################################################################
 def UQTkOptimizeEta(pc_start, y, x, etas, niter, nfolds, mindex_growth, verbose, plot=False):
  """
@@ -600,7 +605,7 @@ def UQTkOptimizeEta(pc_start, y, x, etas, niter, nfolds, mindex_growth, verbose,
  full_basis_size = pc_start.GetNumberPCTerms()
  else:
  full_basis_size = uqtkpce.PCSet("NISPnoq", pc_start.GetOrder() + niter -1, pc_start.GetNDim(), pc_start.GetPCType(), pc_start.GetAlpha(), pc_start.GetBeta()).GetNumberPCTerms()
- 
+
  # loop through each fold
  for i in range(nfolds):
  # retrieve training and validation data
@@ -699,8 +704,10 @@ def UQTkEvalBCS(pc_model, f_evaluations, samplepts, sigma2, eta, regparams, verb
  verbose: Flag for optional print statements
 
  Output:
- c_k: 1D NumPy array of nonzero coefficients
- used_mi_np: NumPy array with the multiindex containing only terms selected by BCS
+ c_k: 1D NumPy array of nonzero coefficients
+ used_mi_np: NumPy array with the multiindex containing only terms
+ selected by BCS
+ sigma2_reestimated: Noise variance reestimated by BCS
  """
  # Configure BCS parameters to defaults
  adaptive = 0 # Flag for adaptive CS, using a generative basis, set to 0 or 1
@@ -721,6 +728,9 @@ def UQTkEvalBCS(pc_model, f_evaluations, samplepts, sigma2, eta, regparams, verb
  for i2 in range(lambda_init.shape[0]):
  lam_uqtk.assign(i2, lambda_init[i2])
 
+ # uqtkarray for sigma2
+ sigma2_array=uqtkarray.dblArray1D(1,sigma2)
+
  #UQTk array for the basis terms evaluated at the sample points
  psi_uqtk = uqtkarray.dblArray2D()
  pc_model.EvalBasisAtCustPts(sam_uqtk, psi_uqtk)
@@ -734,10 +744,10 @@ def UQTkEvalBCS(pc_model, f_evaluations, samplepts, sigma2, eta, regparams, verb
  #vector
  alpha = uqtkarray.dblArray1D() # inverse variance of the coefficient priors,
  # updated through the algorithm
- Sig = uqtkarray.dblArray2D() # re-estimated noise variance
+ Sig = uqtkarray.dblArray2D() # covariance matrix of the weights
 
  # Run BCS through the c++ implementation
- bcs.WBCS(psi_uqtk, y, sigma2, eta, lam_uqtk, adaptive, optimal, scale,\
+ bcs.WBCS(psi_uqtk, y, sigma2_array, eta, lam_uqtk, adaptive, optimal, scale,\
  bcs_verbose, weights, used, errbars, basis, alpha, Sig)
 
  # Print result of the BCS iteration
@@ -757,8 +767,11 @@ def UQTkEvalBCS(pc_model, f_evaluations, samplepts, sigma2, eta, regparams, verb
  uqtkarray.subMatrix_row_int(mi_uqtk, used, used_mi_uqtk)
  used_mi_np=uqtkarray.uqtk2numpy(used_mi_uqtk)
 
+ # convert the returned noise variance to a scalar
+ sigma2_reestimated=sigma2_array[0]
+
  # Return coefficients and their locations with respect to the basis terms
- return c_k, used_mi_np
+ return c_k, used_mi_np, sigma2_reestimated
 ################################################################################
 def UQTkCallBCSDirect(vdm_np, rhs_np, sigma2, eta=1.e-8, regparams_np=None, verbose=False):
  """
@@ -813,11 +826,14 @@ def UQTkCallBCSDirect(vdm_np, rhs_np, sigma2, eta=1.e-8, regparams_np=None, verb
  regparams_np = np.array([])
  elif type(regparams_np)==int or type(regparams_np)==float:
  regparams_np = regparams_np*np.ones((n_basis_terms,))
+
  # Convert to UQTk array
  lam_uqtk=uqtkarray.dblArray1D(regparams_np.shape[0])
  for i2 in range(regparams_np.shape[0]):
  lam_uqtk.assign(i2, regparams_np[i2])
 
+ sigma2_array = uqtkarray.dblArray1D(1,sigma2)
+
 
  # UQTk arrays for outputs
  weights = uqtkarray.dblArray1D() # sparse weights
@@ -831,7 +847,7 @@ def UQTkCallBCSDirect(vdm_np, rhs_np, sigma2, eta=1.e-8, regparams_np=None, verb
  Sig = uqtkarray.dblArray2D() # re-estimated noise variance
 
  # Run BCS through the c++ implementation
- bcs.WBCS(psi_uqtk, rhs_uqtk, sigma2, eta, lam_uqtk, adaptive, optimal, scale,\
+ bcs.WBCS(psi_uqtk, rhs_uqtk, sigma2_array, eta, lam_uqtk, adaptive, optimal, scale,\
  bcs_verbose, weights, used, errbars, basis, alpha, Sig)
 
  # Print result of the BCS iteration

PyUQTk/pytests/PyBCSTest1D.py

@@ -108,8 +108,7 @@
  value = 1.0/(1.0 + x.at(i,0)*x.at(i,0));
  y.assign(i,value)
 
-#sigma = uqtkarray.dblArray1D(1,1e-8)
-sigma = 1e-8
+sigma = uqtkarray.dblArray1D(1,1e-8)
 eta = 1e-8
 lambda_init = uqtkarray.dblArray1D()
 scale = 0.1
@@ -119,14 +118,13 @@
 basis = uqtkarray.dblArray1D()
 alpha = uqtkarray.dblArray1D()
 used = uqtkarray.intArray1D()
-#_lambda = uqtkarray.dblArray1D(1,0.0)
-_lambda=0.0
+Sig=uqtkarray.dblArray2D()
 
 adaptive = 1
 optimal = 1
 verbose = 0
 
-bcs.BCS(Phi,y,sigma,eta,lambda_init,adaptive,optimal,scale,verbose,weights,used,errbars,basis,alpha,_lambda)
+bcs.WBCS(Phi,y,sigma,eta,lambda_init,adaptive,optimal,scale,verbose,weights,used,errbars,basis,alpha,Sig)
 
 uqtkarray.printarray(weights)
 uqtkarray.printarray(used)

PyUQTk/pytests/PyBCSTest2D.py

@@ -130,8 +130,8 @@
 Phi = uqtkarray.dblArray2D()
 pcmodel.EvalBasisAtCustPts(x,Phi)
 
-#sigma = uqtkarray.dblArray1D(1,1e-8)
-sigma = 1e-8
+sigma = uqtkarray.dblArray1D(1,1e-8)
+
 eta = 1e-12
 lambda_init = uqtkarray.dblArray1D()
 scale = 0.1
@@ -141,14 +141,13 @@
 basis = uqtkarray.dblArray1D()
 alpha = uqtkarray.dblArray1D()
 used = uqtkarray.intArray1D()
-#_lambda = uqtkarray.dblArray1D(1,0.0)
-_lambda=0.0
+Sig=uqtkarray.dblArray2D()
 
 adaptive = 1
 optimal = 1
 verbose = 0
 
-bcs.BCS(Phi,y,sigma,eta,lambda_init,adaptive,optimal,scale,verbose,weights,used,errbars,basis,alpha,_lambda)
+bcs.WBCS(Phi,y,sigma,eta,lambda_init,adaptive,optimal,scale,verbose,weights,used,errbars,basis,alpha,Sig)
 
 uqtkarray.printarray(weights)
 uqtkarray.printarray(used)

cpp/lib/bcs/bcs.cpp

@@ -88,9 +88,9 @@
 // % errbars: one standard deviation around the sparse weights
 // % basis: if adaptive==1, then basis = the next projection vector, see [2]
 // % alpha: inverse variance of the coefficient priors, updated through the algorithm, see [1]
-// % Sig: re-estimated noise variance
+// % Sig: covariance matrix of the weights
 
-void WBCS(Array2D<double> &PHI, Array1D<double> &y, double &sigma2,
+void WBCS(Array2D<double> &PHI, Array1D<double> &y, Array1D<double> &sigma2,
  double eta, Array1D<double> &lambda_init,
  int adaptive, int optimal, double scale, int verbose,
  Array1D<double> &weights, Array1D<int> &used,
@@ -121,27 +121,27 @@ void WBCS(Array2D<double> &PHI, Array1D<double> &y, double &sigma2,
  double maxr = maxVal(ratio,&indx) ;
  Array1D<int>index(1,indx);
  alpha.Resize(1) ;
- alpha(0)= PHI2(index(0))/(maxr-sigma2);
+ alpha(0)= PHI2(index(0))/(maxr-sigma2(0));
 
  // Compute initial mu, Sig, S, Q
  Array2D<double> phi(n,1,0.0);
  for (int i = 0; i<n; i++ ) phi(i,0) = PHI(i,index(0));
  double Hessian=alpha(0);
- for (int i = 0; i<n; i++ ) Hessian += phi(i,0)*phi(i,0)/sigma2;
+ for (int i = 0; i<n; i++ ) Hessian += phi(i,0)*phi(i,0)/sigma2(0);
 
  Sig.Resize(1,1,1./(Hessian+1.e-12));
- double dtmp1=Sig(0,0)*PHIy(index(0))/sigma2;
+ double dtmp1=Sig(0,0)*PHIy(index(0))/sigma2(0);
 
  Array1D<double> mu(1,dtmp1);
  Array1D<double> left(m,0.0),phitmp(n);
  for ( int i=0; i<n; i++ )phitmp(i)=phi(i,0);
- prodAlphaMatVec(PHIT, phitmp, 1.0/sigma2, left);
+ prodAlphaMatVec(PHIT, phitmp, 1.0/sigma2(0), left);
 
  Array1D<double> S(m), Q(m);
  for (int j=0; j<m; j++)
  {
- S(j) = PHI2(j)/sigma2-Sig(0,0)*left(j)*left(j);
- Q(j) = PHIy(j)/sigma2-Sig(0,0)*PHIy(index(0))*left(j)/sigma2;
+ S(j) = PHI2(j)/sigma2(0)-Sig(0,0)*left(j)*left(j);
+ Q(j) = PHIy(j)/sigma2(0)-Sig(0,0)*PHIy(index(0))*left(j)/sigma2(0);
  }
 
  // Keep track of the positions selected during the iterations
@@ -154,23 +154,23 @@ void WBCS(Array2D<double> &PHI, Array1D<double> &y, double &sigma2,
  int count = 0;
  for ( count=0; count<MAX_IT; count++ )
  {
- //if ( (count%10)==0 )
- // cout << "Iteration # " << count << endl;
+ if (verbose > 0)
+   cout << "Iteration # " << count << endl;
  Array1D<double> s,q;
  s=S; q=Q;
  for ( int i = 0; i< (int) index.XSize(); i++){
  s(index(i)) = alpha(i)*S(index(i))/(alpha(i)-S(index(i)));
  q(index(i)) = alpha(i)*Q(index(i))/(alpha(i)-S(index(i)));
  }
 
+ // if lambda_init is an empty array, autopopulate it
  if (lambda_init.XSize() == 0)
  {
  double suma=0.0;
  for ( int i = 0; i< (int) alpha.XSize(); i++) suma += 1.0/alpha(i);
  double lambda_ = 2*( index.XSize() - 1 ) / suma;
  lambda_init.Resize(m,lambda_);
  }
-
  Array1D<double> A(m,0.0), B(m,0.0), C(m,0.0);
  Array1D<double> theta(m,0.0), discriminant(m,0.0), nextAlphas(m,0.0) , theta_lambda(m,0.0);
  for ( int i=0; i<m; i++ ){
@@ -181,6 +181,7 @@ void WBCS(Array2D<double> &PHI, Array1D<double> &y, double &sigma2,
  discriminant(i) = pow(B(i),2) - 4.0*A(i)*C(i);
  nextAlphas(i) = (-B(i) - sqrt(discriminant(i)) ) / (2.0*A(i));
  theta_lambda(i)=theta(i)-lambda_init(i);
+
  }
 
  // Choose the next alpha that maximizes marginal likelihood
@@ -191,7 +192,6 @@ void WBCS(Array2D<double> &PHI, Array1D<double> &y, double &sigma2,
  if (ig0.XSize()==0)
  ig0.PushBack(0.0);
 
-
  // Indices for reestimation
  Array1D<int> ire,foo,which ;
  intersect(ig0,index,ire,foo,which);
@@ -286,7 +286,7 @@ void WBCS(Array2D<double> &PHI, Array1D<double> &y, double &sigma2,
  Array1D<double> tmp1, comm;
 
  prodAlphaMatVec (phi,Sigi,1.0, tmp1);
- prodAlphaMatTVec(PHI,tmp1,1.0/sigma2,comm);
+ prodAlphaMatTVec(PHI,tmp1,1.0/sigma2(0),comm);
  addVecAlphaVecPow(S,ki, comm,2);
  addVecAlphaVecPow(Q,ki*mui,comm,1);
  alpha(which(0)) = Alpha;
@@ -305,7 +305,7 @@ void WBCS(Array2D<double> &PHI, Array1D<double> &y, double &sigma2,
  double mui = Sigii*Q(idx);
  Array1D<double> comm1,tmp1;
 
- prodAlphaMatTVec(phi,phii,1.0/sigma2,tmp1);
+ prodAlphaMatTVec(phi,phii,1.0/sigma2(0),tmp1);
  prodAlphaMatVec (Sig,tmp1,1.0, comm1);
  Array1D<double> ei;
  ei=phii;
@@ -326,7 +326,7 @@ void WBCS(Array2D<double> &PHI, Array1D<double> &y, double &sigma2,
  mu.PushBack(mui);
 
  Array1D<double> comm2(m,0.0);
- prodAlphaMatTVec(PHI,ei,1.0/sigma2,comm2);
+ prodAlphaMatTVec(PHI,ei,1.0/sigma2(0),comm2);
 
  addVecAlphaVecPow(S,-Sigii,comm2,2);
  addVecAlphaVecPow(Q,-mui,comm2,1);
@@ -363,7 +363,7 @@ void WBCS(Array2D<double> &PHI, Array1D<double> &y, double &sigma2,
  delCol(mu,which(0));
 
  Array1D<double> comm,tmp1 ;
- prodAlphaMatVec(phi,Sigi,1.0/sigma2,tmp1);
+ prodAlphaMatVec(phi,Sigi,1.0/sigma2(0),tmp1);
  prodAlphaMatTVec(PHI,tmp1,1.0,comm);
  addVecAlphaVecPow(S,1.0/Sigii,comm,2);
  addVecAlphaVecPow(Q,mui/Sigii,comm,1);
@@ -404,11 +404,10 @@ void WBCS(Array2D<double> &PHI, Array1D<double> &y, double &sigma2,
  }
 
 
- sigma2 = sum1/((double) (n-index.XSize())+sum3) ;
+ sigma2(0) = sum1/((double) (n-index.XSize())+sum3) ;
  errbars.Resize(Sig.XSize()) ;
  for ( int i = 0; i<(int) errbars.XSize(); i++) errbars(i) = sqrt(Sig(i,i));
 
-
  // Generate a basis for adaptive CS
  if ( adaptive == 1 ){
 
@@ -437,7 +436,7 @@ void WBCS(Array2D<double> &PHI, Array1D<double> &y, double &sigma2,
 
  else{
  Array2D<double> tmp1;
- prodAlphaMatTMat(phi,phi,1.0/sigma2,tmp1);
+ prodAlphaMatTMat(phi,phi,1.0/sigma2(0),tmp1);
  double tmp1m = 0.0 ;
  for ( int i = 0; i < (int) tmp1.XSize() ; i++ )
  tmp1m += tmp1(i,i);
@@ -468,7 +467,8 @@ void WBCS(Array2D<double> &PHI, Array1D<double> &y, double &sigma2,
  }
 
  }
- //printf("BCS algorithm converged, # iterations : %d \n",count);
+ if (verbose > 0)
+ printf("BCS algorithm converged, # iterations : %d \n",count);
  return ;
 
 }
@@ -518,10 +518,13 @@ void BCS(Array2D<double> &PHI, Array1D<double> &y, double &sigma2,
  Array1D<double> &errbars, Array1D<double> &basis,
  Array1D<double> &alpha, double &lambda)
 {
+ printf("This function is obsolete. Please rework code to use WBCS.");
  int m = (int) PHI.YSize() ;
  Array2D<double> Sig;
- WBCS(PHI, y, sigma2, eta, lambda_init, adaptive, optimal, scale, verbose,
- weights, used, errbars, basis, alpha, Sig);
+ Array1D<double> sigma2_array(1, sigma2);
+ WBCS(PHI, y, sigma2_array, eta, lambda_init, adaptive, optimal, scale, verbose,
+ weights, used, errbars, basis, alpha, Sig);
+
 }
 
 // void BCS(Array2D<double> &PHI, Array1D<double> &y, Array1D<double> &sigma2,

cpp/lib/bcs/bcs.h

@@ -61,7 +61,7 @@
 /// \param[out] basis: : if adaptive==1, then this is the next projection vector, see \cite Ji:2008
 /// \param[out] alpha: : estimated sparse hyperparameters (1/gamma), see \cite Babacan:2010
 /// \param[out] Sig : covariance matrix of the weights
-void WBCS(Array2D<double> &PHI, Array1D<double> &y, double &sigma2,
+void WBCS(Array2D<double> &PHI, Array1D<double> &y, Array1D<double> &sigma2,
  double eta, Array1D<double> &lambda_init,
  int adaptive, int optimal, double scale, int verbose,
  Array1D<double> &weights, Array1D<int> &used,