diff --git a/fastlib/trunk/contrib/dongryel/gp_regression/block_coord_descent.h b/fastlib/trunk/contrib/dongryel/gp_regression/block_coord_descent.h
index 806c724498..1b10e50d30 100644
--- a/fastlib/trunk/contrib/dongryel/gp_regression/block_coord_descent.h
+++ b/fastlib/trunk/contrib/dongryel/gp_regression/block_coord_descent.h
@@ -69,7 +69,7 @@ class BlockCoordDescentResult: boost::noncopyable {
 
     void Init(bool initialize_result,
               const fl::data::MonolithicPoint<double> &right_hand) {
-      if(initialize_result) {
+      if (initialize_result) {
         gradient_.Init(right_hand.length());
         solution_.Init(right_hand.length());
       }
@@ -156,7 +156,7 @@ class BlockCoordDescentDelta: boost::noncopyable {
     }
 
     template < typename Matrix, typename KernelType,
-             typename double >
+    typename double >
     void Update(const Matrix &table, const KernelType &kernel,
                 int new_index, double self_kernel_value,
                 const fl::data::MonolithicPoint<double> &gradient) {
@@ -166,7 +166,7 @@ class BlockCoordDescentDelta: boost::noncopyable {
       beta.Init(active_set_.size());
       beta.SetZero();
 
-      for(int i = 0; i < active_set_.size(); i++) {
+      for (int i = 0; i < active_set_.size(); i++) {
 
         // Compute the kernel value with the current active point.
         int active_point_id = active_set_[i];
@@ -182,7 +182,7 @@ class BlockCoordDescentDelta: boost::noncopyable {
 
         // Compute the beta by multiplying the old inverse by
         // the the column we are adding into the kernel matrix.
-        for(int j = 0; j < active_set_.size(); j++) {
+        for (int j = 0; j < active_set_.size(); j++) {
           beta[j] += kernel_matrix_.get(i, active_set_.size()) *
                      inverse_.get(j, i);
         }
@@ -193,23 +193,23 @@ class BlockCoordDescentDelta: boost::noncopyable {
                          self_kernel_value);
 
       // Compute eta.
-      for(int i = 0; i < active_set_.size(); i++) {
+      for (int i = 0; i < active_set_.size(); i++) {
         eta += kernel_matrix_.get(active_set_.size(), i) * beta[i];
       }
       eta = 1.0 / (self_kernel_value - eta);
 
       // Update the inverse.
-      for(int j = 0; j < active_set_.size() + 1; j++) {
-        for(int i = 0; i < active_set_.size() + 1; i++) {
+      for (int j = 0; j < active_set_.size() + 1; j++) {
+        for (int i = 0; i < active_set_.size() + 1; i++) {
           double increment = eta;
 
-          if(i < active_set_.size()) {
+          if (i < active_set_.size()) {
             increment *= beta[i];
           }
           else {
             increment *= (-1);
           }
-          if(j < active_set_.size()) {
+          if (j < active_set_.size()) {
             increment *= beta[j];
           }
           else {
@@ -223,12 +223,12 @@ class BlockCoordDescentDelta: boost::noncopyable {
       // Update the delta solution by computing the dot product
       // between beta and the current gradient.
       double factor = 0;
-      for(int i = 0; i < active_set_.size(); i++) {
+      for (int i = 0; i < active_set_.size(); i++) {
         factor += beta[i] * gradient[active_set_[i]];
       }
       factor -= gradient[new_index];
       factor *= eta;
-      for(int i = 0; i < active_set_.size(); i++) {
+      for (int i = 0; i < active_set_.size(); i++) {
         delta_solution_[active_set_[i]] -= factor * beta[i];
       }
       delta_solution_[new_index] += factor;
@@ -252,7 +252,7 @@ class BlockCoordDescentDelta: boost::noncopyable {
       temp_gradient_info_.SetZero();
       inactive_set_.resize(num_points_);
       active_set_.resize(0);
-      for(int i = 0; i < num_points_; i++) {
+      for (int i = 0; i < num_points_; i++) {
         inactive_set_[i] = i;
       }
       inverse_.SetZero();
@@ -292,9 +292,9 @@ class SolveSubProblem {
       const std::vector<int> &active_set = delta.active_set();
       fl::dense::Matrix<double, false> &kernel_matrix = delta.kernel_matrix();
 
-      for(int j = 0; j < active_set.size(); j++) {
+      for (int j = 0; j < active_set.size(); j++) {
         int column_index = active_set[j];
-        for(int i = 0; i < active_set.size(); i++) {
+        for (int i = 0; i < active_set.size(); i++) {
           int row_index = active_set[i];
           kernel_matrix.set(
             i, j, fl::ml::KernelValue::Compute(
@@ -310,10 +310,10 @@ class SolveSubProblem {
       fl::data::MonolithicPoint<double> residual;
       residual.Alias(delta.temp_gradient_info().ptr(), kernel_matrix.n_rows());
 
-      for(int i = 0; i < active_set.size(); i++) {
+      for (int i = 0; i < active_set.size(); i++) {
         int row_index = active_set[i];
         residual[i] = right_hand_side[row_index];
-        for(int j = 0; j < active_set.size(); j++) {
+        for (int j = 0; j < active_set.size(); j++) {
           int col_index = active_set[j];
           residual[i] -= kernel_matrix.get(i, j) * current_solution[col_index];
         }
@@ -326,7 +326,7 @@ class SolveSubProblem {
       fl::data::MonolithicPoint<double> &delta_solution =
         delta.delta_solution();
       delta_solution.SetZero();
-      for(int i = 0; i < kernel_matrix.n_rows(); i++) {
+      for (int i = 0; i < kernel_matrix.n_rows(); i++) {
         int row_index = active_set[i];
         delta_solution[row_index] = -tmp_vector_[i];
       }
@@ -365,8 +365,8 @@ class SelectActiveSetTrait<fl::ml::GpRegressionComputation::GREEDY_BC> {
           max_random_set_size,
           (int)inactive_set.size());
 
-      for(int i = ((int) inactive_set.size()) - 1;
-          i >= ((int) inactive_set.size()) - (*random_set_size); i--) {
+      for (int i = ((int) inactive_set.size()) - 1;
+           i >= ((int) inactive_set.size()) - (*random_set_size); i--) {
 
         int random_index = fl::math::Random(0, i + 1);
         std::swap(inactive_set[random_index], inactive_set[i]);
@@ -426,22 +426,22 @@ class SelectActiveSetTrait<fl::ml::GpRegressionComputation::GREEDY_BC> {
         double minimum_value =
           std::numeric_limits<double>::max();
         double cached_self_kernel_value = -1;
-        for(int i = inactive_set.size() - random_set_size;
-            i < inactive_set.size(); i++) {
+        for (int i = inactive_set.size() - random_set_size;
+             i < inactive_set.size(); i++) {
           double kernel_value =
             fl::ml::KernelValue::Compute(
               table, kernel, inactive_set[i], inactive_set[i], true);
           double current_value =
             - fl::math::Sqr(temp_gradient_info[inactive_set[i]]) /
             (2.0 * kernel_value);
-          if(current_value <= minimum_value) {
+          if (current_value <= minimum_value) {
             minimum_value = current_value;
             selected_index = i;
             cached_self_kernel_value = kernel_value;
           }
         }
 
-        if(delta.current_active_set_size() == 0) {
+        if (delta.current_active_set_size() == 0) {
           inverse.set(0, 0, 1.0 / cached_self_kernel_value);
           kernel_matrix.set(0, 0, cached_self_kernel_value);
           delta_solution[inactive_set[selected_index]] =
@@ -466,15 +466,15 @@ class SelectActiveSetTrait<fl::ml::GpRegressionComputation::GREEDY_BC> {
         ChooseSubset(
           result, inactive_set, max_random_set_size, &random_set_size);
 
-        for(int i = inactive_set.size() - random_set_size;
-            i < inactive_set.size(); i++) {
+        for (int i = inactive_set.size() - random_set_size;
+             i < inactive_set.size(); i++) {
 
           // The index of the randomly chosen point.
           int random_point_index = inactive_set[i];
           double dot_product = 0;
 
           // Loop over each point in the current active set.
-          for(int j = 0; j < active_set.size(); j++) {
+          for (int j = 0; j < active_set.size(); j++) {
             int active_point_index = active_set[j];
             double kernel_value =
               fl::ml::KernelValue::Compute(
@@ -487,8 +487,8 @@ class SelectActiveSetTrait<fl::ml::GpRegressionComputation::GREEDY_BC> {
         }
 
       }
-      while(delta.current_active_set_size() < active_set_size &&
-            inactive_set.size() > 0);
+      while (delta.current_active_set_size() < active_set_size &&
+             inactive_set.size() > 0);
     }
 };
 
@@ -521,12 +521,12 @@ class SelectActiveSetTrait<fl::ml::GpRegressionComputation::CYCLIC_BC> {
       std::vector<int> &inactive_set = delta.inactive_set();
       active_set.resize(0);
       inactive_set.resize(0);
-      for(int i = 0; i < std::min(active_set_size, table.n_entries()); i++) {
+      for (int i = 0; i < std::min(active_set_size, table.n_entries()); i++) {
         int active_index = (i + starting_index_) % table.n_entries();
         active_set.push_back(active_index);
       }
-      for(int i = std::min(active_set_size, table.n_entries());
-          i < table.n_entries(); i++) {
+      for (int i = std::min(active_set_size, table.n_entries());
+           i < table.n_entries(); i++) {
         int inactive_index = (i + starting_index_) % table.n_entries();
         inactive_set.push_back(inactive_index);
       }
@@ -591,7 +591,7 @@ class SelectActiveSetTrait<fl::ml::GpRegressionComputation::GRADIENT_BC> {
 
       // Sort the current gradient component by absolute magnitude and
       // select the indices with the greatest ones.
-      for(int i = 0; i < sorted_indices_.size(); i++) {
+      for (int i = 0; i < sorted_indices_.size(); i++) {
         sorted_indices_[i] = i;
       }
 
@@ -603,11 +603,11 @@ class SelectActiveSetTrait<fl::ml::GpRegressionComputation::GRADIENT_BC> {
       std::vector<int> &inactive_set = delta.inactive_set();
       active_set.resize(0);
       inactive_set.resize(0);
-      for(int i = 0; i < std::min(active_set_size, table.n_entries()); i++) {
+      for (int i = 0; i < std::min(active_set_size, table.n_entries()); i++) {
         active_set.push_back(sorted_indices_[i]);
       }
-      for(int i = std::min(active_set_size, table.n_entries());
-          i < table.n_entries(); i++) {
+      for (int i = std::min(active_set_size, table.n_entries());
+           i < table.n_entries(); i++) {
         inactive_set.push_back(sorted_indices_[i]);
       }
 
@@ -639,11 +639,11 @@ class BlockCoordDescent: private boost::noncopyable {
       fl::data::MonolithicPoint<double> &solution = result.solution();
       fl::data::MonolithicPoint<double> &gradient = result.gradient();
 
-      for(int i = 0; i < active_set.size(); i++) {
+      for (int i = 0; i < active_set.size(); i++) {
         solution[active_set[i]] -= delta_solution[active_set[i]];
 
         double dot_product = 0;
-        for(int j = 0; j < active_set.size(); j++) {
+        for (int j = 0; j < active_set.size(); j++) {
           dot_product += kernel_matrix.get(i, j) *
                          delta_solution[active_set[j]];
         }
@@ -674,7 +674,7 @@ class BlockCoordDescent: private boost::noncopyable {
         active_set_size, table.n_entries());
 
       // The main loop.
-      while(!result->IsConverged()) {
+      while (!result->IsConverged()) {
 
         active_set_selection_trait.Select(
           active_set_size, max_random_set_size, table, kernel, right_hand,
diff --git a/fastlib/trunk/contrib/dongryel/gp_regression/dense_matrix_inverse.h b/fastlib/trunk/contrib/dongryel/gp_regression/dense_matrix_inverse.h
index 4b13fdac0a..91ab604f61 100644
--- a/fastlib/trunk/contrib/dongryel/gp_regression/dense_matrix_inverse.h
+++ b/fastlib/trunk/contrib/dongryel/gp_regression/dense_matrix_inverse.h
@@ -24,8 +24,8 @@ class DenseMatrixInverse {
       new_matrix_inverse->Init(previous_inverse.n_rows() + 1,
                                previous_inverse.n_cols() + 1);
 
-      for(int j = 0; j < previous_inverse.n_cols(); j++) {
-        for(int i = 0; i < previous_inverse.n_rows(); i++) {
+      for (int j = 0; j < previous_inverse.n_cols(); j++) {
+        for (int i = 0; i < previous_inverse.n_rows(); i++) {
           new_matrix_inverse->set(
             i, j, previous_inverse.get(i, j) +
             inverse_times_new_column[i] *
@@ -33,7 +33,7 @@ class DenseMatrixInverse {
         }
       }
 
-      for(int j = 0; j < previous_inverse.n_cols(); j++) {
+      for (int j = 0; j < previous_inverse.n_cols(); j++) {
         new_matrix_inverse->set(
           j,
           previous_inverse.n_cols(), - inverse_times_new_column[j] /
diff --git a/fastlib/trunk/contrib/dongryel/gp_regression/dictionary_dev.h b/fastlib/trunk/contrib/dongryel/gp_regression/dictionary_dev.h
index c003bbeff9..4784499c0c 100644
--- a/fastlib/trunk/contrib/dongryel/gp_regression/dictionary_dev.h
+++ b/fastlib/trunk/contrib/dongryel/gp_regression/dictionary_dev.h
@@ -28,10 +28,10 @@ Dictionary::Dictionary() {
 }
 
 Dictionary::~Dictionary() {
-  if(current_kernel_matrix_ != NULL) {
+  if (current_kernel_matrix_ != NULL) {
     delete current_kernel_matrix_;
   }
-  if(current_kernel_matrix_inverse_ != NULL) {
+  if (current_kernel_matrix_inverse_ != NULL) {
     delete current_kernel_matrix_inverse_;
   }
 }
@@ -57,8 +57,8 @@ void Dictionary::inactive_indices(
 
   // Scan the in_dictionary list and build the inactive index set.
   inactive_indices_out->resize(0);
-  for(int i = 0; i < in_dictionary_.size(); i++) {
-    if(in_dictionary_[i] == false) {
+  for (int i = 0; i < in_dictionary_.size(); i++) {
+    if (in_dictionary_[i] == false) {
       inactive_indices_out->push_back(i);
     }
   }
@@ -83,12 +83,12 @@ void Dictionary::UpdateDictionary_(
     current_kernel_matrix_->n_rows() + 1,
     current_kernel_matrix_->n_cols() + 1);
 
-  for(int j = 0; j < current_kernel_matrix_->n_cols(); j++) {
-    for(int i = 0; i < current_kernel_matrix_->n_rows(); i++) {
+  for (int j = 0; j < current_kernel_matrix_->n_cols(); j++) {
+    for (int i = 0; i < current_kernel_matrix_->n_rows(); i++) {
       new_kernel_matrix->set(i, j, current_kernel_matrix_->get(i, j));
     }
   }
-  for(int j = 0; j < current_kernel_matrix_->n_cols(); j++) {
+  for (int j = 0; j < current_kernel_matrix_->n_cols(); j++) {
     new_kernel_matrix->set(
       j, current_kernel_matrix_->n_cols(), new_column_vector[j]);
     new_kernel_matrix->set(
@@ -138,11 +138,11 @@ void Dictionary::AddBasis(
   const std::vector<double> &new_column_vector_in,
   double self_value) {
 
-  if(new_column_vector_in.size() > 0) {
+  if (new_column_vector_in.size() > 0) {
 
     Vector new_column_vector;
     new_column_vector.Init(new_column_vector_in.size());
-    for(int i = 0; i < new_column_vector.length(); i++) {
+    for (int i = 0; i < new_column_vector.length(); i++) {
       new_column_vector[i] = new_column_vector_in[i];
     }
 
@@ -160,7 +160,7 @@ void Dictionary::AddBasis(
 
     // If the projection error is above the threshold, add it to the
     // dictionary.
-    if(projection_error > adding_threshold_) {
+    if (projection_error > adding_threshold_) {
       UpdateDictionary_(
         new_point_index,
         new_column_vector,
@@ -191,7 +191,7 @@ void Dictionary::Init(const Matrix *table_in) {
   // Generate a random permutation and initialize the inital
   // dictionary which consists of the first random point.
   random_permutation_.resize(table_in->n_cols());
-  for(int i = 0; i < table_in->n_cols(); i++) {
+  for (int i = 0; i < table_in->n_cols(); i++) {
     random_permutation_[i] = i;
     in_dictionary_[i] = false;
     training_index_to_dictionary_position_[i] = -1;
diff --git a/fastlib/trunk/contrib/dongryel/gp_regression/kernel_linear_operator.h b/fastlib/trunk/contrib/dongryel/gp_regression/kernel_linear_operator.h
index a795376516..793aa4e777 100644
--- a/fastlib/trunk/contrib/dongryel/gp_regression/kernel_linear_operator.h
+++ b/fastlib/trunk/contrib/dongryel/gp_regression/kernel_linear_operator.h
@@ -45,7 +45,7 @@ class DotProductTrait<false> {
 };
 
 template < typename TableType, typename KernelType, bool do_centering,
-         bool dotproduct_selfcase_special = false >
+bool dotproduct_selfcase_special = false >
 class KernelLinearOperator: public LinearOperator {
 
   private:
@@ -74,23 +74,23 @@ class KernelLinearOperator: public LinearOperator {
       comm_ = &comm_in;
       map_ = &map_in;
 
-      if(do_centering) {
+      if (do_centering) {
         average_row_.Init(table_in.n_entries());
       }
       average_ = 0;
 
-      if(do_centering) {
+      if (do_centering) {
 
         // Precompute the average. This is a naive way of computing it.
-        for(int i = 0; i < table_in.n_entries(); i++) {
+        for (int i = 0; i < table_in.n_entries(); i++) {
           double average_for_i_th_point = 0;
-          for(int j = 0; j < table_in.n_entries(); j++) {
+          for (int j = 0; j < table_in.n_entries(); j++) {
             average_for_i_th_point += kernel_value(i, j);
           }
           average_for_i_th_point /= ((double) table_in.n_entries());
           average_row_[i] = average_for_i_th_point;
         }
-        for(int i = 0; i < table_in.n_entries(); i++) {
+        for (int i = 0; i < table_in.n_entries(); i++) {
           average_ += average_row_[i];
         }
         average_ /= ((double) table_in.n_entries());
@@ -125,11 +125,11 @@ class KernelLinearOperator: public LinearOperator {
 
       prods.PutScalar(0);
 
-      for(int j = 0; j < table_->n_entries(); j++) {
-        for(int i = 0; i < table_->n_entries(); i++) {
+      for (int j = 0; j < table_->n_entries(); j++) {
+        for (int i = 0; i < table_->n_entries(); i++) {
           double pair_kernel_value = (do_centering) ?
                                      centered_kernel_value(i, j) : kernel_value(i, j);
-          for(int k = 0; k < vecs.NumVectors(); k++) {
+          for (int k = 0; k < vecs.NumVectors(); k++) {
             prods.Pointers()[k][i] += pair_kernel_value * vecs.Pointers()[k][j];
           }
         }
@@ -152,7 +152,7 @@ class KernelLinearOperator: public LinearOperator {
 
 template<typename TableType, typename KernelType, bool do_centering>
 class OperatorTraits < double, Epetra_MultiVector,
-    KernelLinearOperator<TableType, KernelType, do_centering> > {
+      KernelLinearOperator<TableType, KernelType, do_centering> > {
   public:
 
     static void Apply(const Epetra_Operator& Op,
diff --git a/fastlib/trunk/contrib/dongryel/gp_regression/lbfgs_dev.h b/fastlib/trunk/contrib/dongryel/gp_regression/lbfgs_dev.h
index e976ddfce8..d62f7ef476 100644
--- a/fastlib/trunk/contrib/dongryel/gp_regression/lbfgs_dev.h
+++ b/fastlib/trunk/contrib/dongryel/gp_regression/lbfgs_dev.h
@@ -60,7 +60,7 @@ double Lbfgs<FunctionType>::ChooseScalingFactor_(
   const Vector &gradient) {
 
   double scaling_factor = 1.0;
-  if(iteration_num > 0) {
+  if (iteration_num > 0) {
     int previous_pos = (iteration_num - 1) % num_basis_;
     Vector s_basis;
     Vector y_basis;
@@ -113,7 +113,7 @@ bool Lbfgs<FunctionType>::LineSearch_(
     la::Dot(gradient, search_direction);
 
   // If it is not a descent direction, just report failure.
-  if(initial_search_direction_dot_gradient > 0.0) {
+  if (initial_search_direction_dot_gradient > 0.0) {
     return false;
   }
 
@@ -131,7 +131,7 @@ bool Lbfgs<FunctionType>::LineSearch_(
   const double inc = 2.1;
   const double dec = 0.5;
   double width = 0;
-  for(; ;) {
+  for (; ;) {
 
     // Perform a step and evaluate the gradient and the function
     // values at that point.
@@ -141,7 +141,7 @@ bool Lbfgs<FunctionType>::LineSearch_(
     function_->Gradient(new_iterate_tmp_, &gradient);
     num_iterations++;
 
-    if(function_value > initial_function_value + step_size *
+    if (function_value > initial_function_value + step_size *
         linear_approx_function_value_decrease) {
       width = dec;
     }
@@ -151,12 +151,12 @@ bool Lbfgs<FunctionType>::LineSearch_(
       double search_direction_dot_gradient =
         la::Dot(gradient, search_direction);
 
-      if(search_direction_dot_gradient < param_.wolfe() *
+      if (search_direction_dot_gradient < param_.wolfe() *
           initial_search_direction_dot_gradient) {
         width = inc;
       }
       else {
-        if(search_direction_dot_gradient > -param_.wolfe() *
+        if (search_direction_dot_gradient > -param_.wolfe() *
             initial_search_direction_dot_gradient) {
           width = dec;
         }
@@ -168,13 +168,13 @@ bool Lbfgs<FunctionType>::LineSearch_(
 
     // Terminate when the step size gets too small or too big or it
     // exceeds the max number of iterations.
-    if(step_size < param_.min_step()) {
+    if (step_size < param_.min_step()) {
       return false;
     }
-    if(step_size > param_.max_step()) {
+    if (step_size > param_.max_step()) {
       return false;
     }
-    if(num_iterations >= param_.max_line_search()) {
+    if (num_iterations >= param_.max_line_search()) {
       return false;
     }
 
@@ -204,7 +204,7 @@ void Lbfgs<FunctionType>::SearchDirection_(
 
 
   int limit = std::max(iteration_num - num_basis_, 0);
-  for(int i = iteration_num - 1; i >= limit; i--) {
+  for (int i = iteration_num - 1; i >= limit; i--) {
     int translated_position = i % num_basis_;
     Vector y_basis, s_basis;
     s_lbfgs_.MakeColumnVector(translated_position, &s_basis);
@@ -214,7 +214,7 @@ void Lbfgs<FunctionType>::SearchDirection_(
                                      la::Dot(s_basis, q);
   }
   la::ScaleOverwrite(scaling_factor, q, search_direction);
-  for(int i = limit; i <= iteration_num - 1; i++) {
+  for (int i = limit; i <= iteration_num - 1; i++) {
     int translated_position = i % num_basis_;
     Vector y_basis, s_basis;
     s_lbfgs_.MakeColumnVector(translated_position, &s_basis);
@@ -254,7 +254,7 @@ double Lbfgs<FunctionType>::Evaluate_(
   // value encountered during the optimization.
   double function_value = function_->Evaluate(iterate);
 
-  if(function_value < min_point_iterate_.second) {
+  if (function_value < min_point_iterate_.second) {
     min_point_iterate_.first.CopyValues(iterate);
     min_point_iterate_.second = function_value;
   }
@@ -311,11 +311,11 @@ bool Lbfgs<FunctionType>::Optimize(int num_iterations,
 
   // The main optimization loop.
   int it_num;
-  for(it_num = 0; optimize_until_convergence ||
-      it_num < num_iterations; it_num++) {
+  for (it_num = 0; optimize_until_convergence ||
+       it_num < num_iterations; it_num++) {
 
     // Break when the norm of the gradient becomes too small.
-    if(GradientNormTooSmall_(gradient)) {
+    if (GradientNormTooSmall_(gradient)) {
       break;
     }
 
@@ -336,7 +336,7 @@ bool Lbfgs<FunctionType>::Optimize(int num_iterations,
       LineSearch_(function_value, *iterate, gradient, search_direction,
                   step_size);
 
-    if(search_is_success == false) {
+    if (search_is_success == false) {
       break;
     }
 
diff --git a/fastlib/trunk/contrib/dongryel/gp_regression/linear_operator.h b/fastlib/trunk/contrib/dongryel/gp_regression/linear_operator.h
index 42cd22be67..b4c792062c 100644
--- a/fastlib/trunk/contrib/dongryel/gp_regression/linear_operator.h
+++ b/fastlib/trunk/contrib/dongryel/gp_regression/linear_operator.h
@@ -110,8 +110,8 @@ class LinearOperator: public virtual Epetra_Operator {
 
     void PrintDebug(const char *name = "", FILE *stream = stderr) const {
       fprintf(stream, "----- MATRIX ------: %s\n", name);
-      for(int r = 0; r < this->n_rows(); r++) {
-        for(int c = 0; c < this->n_cols(); c++) {
+      for (int r = 0; r < this->n_rows(); r++) {
+        for (int c = 0; c < this->n_cols(); c++) {
           fprintf(stream, "%+3.3f ", this->get(r, c));
         }
         fprintf(stream, "\n");
diff --git a/fastlib/trunk/contrib/dongryel/gp_regression/sg_gp_regression.h b/fastlib/trunk/contrib/dongryel/gp_regression/sg_gp_regression.h
index f9f11bf668..0c98012d9e 100644
--- a/fastlib/trunk/contrib/dongryel/gp_regression/sg_gp_regression.h
+++ b/fastlib/trunk/contrib/dongryel/gp_regression/sg_gp_regression.h
@@ -38,14 +38,27 @@ class SparseGreedyGprModel {
       const Vector &right_hand_side_in,
       Vector *solution_out) const;
 
+    void ExtractWeightedTargetSubset_(
+      const ml::Dictioanay &dictionary_in,
+      const std::vector<double> &additional_kernel_values_in,
+      Vector *target_subset_out) const;
+
     void ExtractTargetSubset_(
       const ml::Dictioanay &dictionary_in,
       Vector *target_subset_out) const;
 
     double QuadraticObjective_(
       const ml::Dictionary &dictionary_in,
+      const std::vector<double> &kernel_values_in,
       bool for_coeffs) const;
 
+    template<typename CovarianceType>
+    void ComputeKernelValues_(
+      const CovarianceType &covariance_in,
+      const std::vector<int> &point_indices_in_dictionary,
+      const Vector &point,
+      Vector *kernel_values_out) const;
+
     template<typename CovarianceType>
     void ComputeKernelValues_(
       const CovarianceType &covariance_in,
@@ -84,12 +97,12 @@ class SparseGreedyGprModel {
     void FinalizeModel();
 
     template<typename CovarianceType>
-    void PredictMean(
+    double PredictMean(
       const CovarianceType &covariance,
       const Vector &point) const;
 
     template<typename CovarianceType>
-    void PredictVariance(
+    double PredictVariance(
       const CovarianceType &covariance,
       const Vector &point) const;
 };
diff --git a/fastlib/trunk/contrib/dongryel/gp_regression/sg_gp_regression_dev.h b/fastlib/trunk/contrib/dongryel/gp_regression/sg_gp_regression_dev.h
index 43079fc7a7..6147a2ea90 100644
--- a/fastlib/trunk/contrib/dongryel/gp_regression/sg_gp_regression_dev.h
+++ b/fastlib/trunk/contrib/dongryel/gp_regression/sg_gp_regression_dev.h
@@ -15,27 +15,45 @@ namespace ml {
 namespace gp_regression {
 
 template<typename CovarianceType>
-void SparseGreedyGprModel::PredictMean(
+double SparseGreedyGprModel::PredictMean(
   const CovarianceType &covariance,
   const Vector &point) const {
 
   // Get the indices of the points in the dictionary.
   const std::vector<int> &point_indices_in_dictionary =
     dictionary_for_error_.point_indices_in_dictionary();
+
+  // Compute the kernel values with the basis points and take the dot
+  // product with the coefficients.
+  Vector kernel_values;
+  ComputeKernelValues_(
+    covariance, point_indices_in_dictionary, point, &kernel_values);
+
+  return la::Dot(kernel_values, coefficients_);
 }
 
 template<typename CovarianceType>
-void SparseGreedyGprModel::PredictVariance(
+double SparseGreedyGprModel::PredictVariance(
   const CovarianceType &covariance,
   const Vector &point) const {
 
   // Get the indices of the points in the dictionary.
   const std::vector<int> &point_indices_in_dictionary =
     dictionary_for_error_.point_indices_in_dictionary();
+
+  // The variance is self-correlation minus the variance explained by
+  // the sparse GPR model.
+  Vector kernel_values;
+  ComputeKernelValues_(
+    covariance, point_indices_in_dictionary, point, &kernel_values);
+
 }
 
 void SparseGreedyGprModel::FinalizeModel() {
-  SolveSystem_(dictionary_, &coefficients_);
+  Vector weighted_target_subset;
+  ExtractWeightedTargetSubset_(
+    dictionary_in, kernel_values_in, &weighted_target_subset);
+  SolveSystem_(dictionary_, weighted_target_subset, &coefficients);
 }
 
 double SparseGreedyGprModel::frobenius_norm_targets() const {
@@ -53,6 +71,30 @@ void SparseGreedyGprModel::SolveSystem_(
               solution_out);
 }
 
+void SparseGreedyGprModel::ExtractWeightedTargetSubset_(
+  const ml::Dictioanay &dictionary_in,
+  const std::vector<double> &additional_kernel_values_in,
+  Vector *target_subset_out) const {
+
+  // Get the indices of the points in the dictionary.
+  const std::vector<int> &point_indices_in_dictionary =
+    dictionary_in.point_indices_in_dictionary();
+
+  target_subset_out->Init(point_indices_in_dictionary.size());
+  for (int i = 0; i < target_subset_out->length(); i++) {
+    int basis_index = point_indices_in_dictionary[i];
+    const std::vector<double> &kernel_values =
+      (i == target_subset_out->length() - 1) ?
+      additional_kernel_values_in : kernel_matrix_columns_[basis_index];
+
+    double dot_product = 0;
+    for (int j = 0; j < dataset_->n_cols(); j++) {
+      dot_product += (*targets_)[ j ] * kernel_values[j];
+    }
+    (*target_subset_out)[i] = dot_product;
+  }
+}
+
 void SparseGreedyGprModel::ExtractTargetSubset_(
   const ml::Dictioanay &dictionary_in,
   Vector *target_subset_out) const {
@@ -62,28 +104,31 @@ void SparseGreedyGprModel::ExtractTargetSubset_(
     dictionary_in.point_indices_in_dictionary();
 
   target_subset_out->Init(point_indices_in_dictionary.size());
-  for(int i = 0; i < target_subset_out->length(); i++) {
+  for (int i = 0; i < target_subset_out->length(); i++) {
     (*target_subset_out)[i] = (*targets_)[ point_indices_in_dictionary[i] ];
   }
 }
 
 double SparseGreedyGprModel::QuadraticObjective_(
   const ml::Dictionary &dictionary_in,
+  const std::vector<double> &kernel_values_in,
   bool for_coeffs) const {
 
   // Compute the objective, -0.5 y^T K^{1} y.
   Vector product;
-  if(for_coeffs) {
+  if (for_coeffs) {
     Vector weighted_target_subset;
+    ExtractWeightedTargetSubset_(
+      dictionary_in, kernel_values_in, &weighted_target_subset);
     SolveSystem_(dictionary_in, weighted_target_subset, &product);
+    return -0.5 * la::Dot(product, weighted_target_subset);
   }
   else {
     Vector target_subset;
     ExtractTargetSubset_(dictionary_in, &target_subset);
     SolveSystem_(dictionary_in, target_subset, &product);
+    return -0.5 * la::Dot(product, target_subset);
   }
-
-  return - 0.5 * la::Dot(product, target_subset);
 }
 
 void SparseGreedyGprModel::FillSquaredKernelMatrix_(
@@ -97,15 +142,15 @@ void SparseGreedyGprModel::FillSquaredKernelMatrix_(
   const std::vector<int> &point_indices_in_dictionary =
     dictionary_.point_indices_in_dictionary();
 
-  for(int i = 0; i < point_indices_in_dictionary.size(); i++) {
+  for (int i = 0; i < point_indices_in_dictionary.size(); i++) {
     int basis_index = point_indices_in_dictionary[i];
     const std::vector<double> &cached_kernel_values =
       kernel_matrix_columns_[basis_index];
 
     // Take the dot product.
     double dot_product = 0;
-    for(int j = 0; j < kernel_values.length(); j++) {
-      if(j == candidate_index) {
+    for (int j = 0; j < kernel_values.length(); j++) {
+      if (j == candidate_index) {
         dot_product += cached_kernel_values[j] *
                        (kernel_values[j] + noise_level_in);
       }
@@ -132,13 +177,29 @@ void SparseGreedyGprModel::FillKernelMatrix_(
 
   // Simply the necessary kernel values. For the self value, add the
   // noise.
-  for(int i = 0; i < point_indices_in_dictionary.size(); i++) {
+  for (int i = 0; i < point_indices_in_dictionary.size(); i++) {
     (*new_column_vector_out)[i] =
       kernel_values[ point_indices_in_dictionary[i] ];
   }
   *new_self_value_out = kernel_values[candidate_index] + noise_level_in;
 }
 
+template<typename CovarianceType>
+void SparseGreedyGprModel::ComputeKernelValues_(
+  const CovarianceType &covariance_in,
+  const std::vector<int> &point_indices_in_dictionary,
+  const Vector &point,
+  Vector *kernel_values_out) const {
+
+  kernel_values_out->Init(point_indices_in_dictionary.size());
+  for (int i = 0; i < kernel_values_out->length(); i++) {
+    int basis_index = point_indices_in_dictionary[i];
+    Vector basis_point;
+    dataset_->MakeColumnVector(i, &basis_point);
+    (*kernel_values_out)[i] = covariance_in.Dot(point, basis_point, false);
+  }
+}
+
 template<typename CovarianceType>
 void SparseGreedyGprModel::ComputeKernelValues_(
   const CovarianceType &covariance_in,
@@ -149,7 +210,7 @@ void SparseGreedyGprModel::ComputeKernelValues_(
   dataset_->MakeColumnVector(candidate_index, &candidate_point);
 
   // Fill out the kernel values sequentially.
-  for(int i = 0; i < dataset_->n_cols(); i++) {
+  for (int i = 0; i < dataset_->n_cols(); i++) {
     Vector point;
     dataset_->MakeColumnVector(i, &point);
     (*kernel_values_out)[i] = covariance_in.Dot(
@@ -172,11 +233,11 @@ double SparseGreedyGprModel::AddOptimalPoint(
   double optimum_value = std::numeric_limits<double>::max();
 
   // Loop over candidates and decide to add the optimal.
-  for(int i = 0; i < candidate_indices.size(); i++) {
+  for (int i = 0; i < candidate_indices.size(); i++) {
 
     // Make a copy of the dictionaries.
     Dictionary dictionary_copy;
-    if(for_coeffs) {
+    if (for_coeffs) {
       dictionary_copy = dictionary_;
     }
     else {
@@ -193,7 +254,7 @@ double SparseGreedyGprModel::AddOptimalPoint(
 
     // Compute the additional quantities to be appended to grow the
     // matrix and update the dictionary.
-    if(for_coeffs) {
+    if (for_coeffs) {
       FillSquaredKernelMatrix_(
         candidate_index, kernel_values, &new_column_vector, &new_self_value);
     }
@@ -205,9 +266,10 @@ double SparseGreedyGprModel::AddOptimalPoint(
       candidate_index, new_column_vector, new_self_value);
 
     // Compute the objective function value for the coefficients.
-    double objective_value = QuadraticObjective_(dictionary_copy, for_coeffs);
+    double objective_value = QuadraticObjective_(
+                               dictionary_copy, kernel_values, for_coeffs);
 
-    if(objective_value < optimum_value) {
+    if (objective_value < optimum_value) {
       optimal_point_index = i;
       optimum_value = objective_value;
     }
@@ -217,7 +279,7 @@ double SparseGreedyGprModel::AddOptimalPoint(
   int final_candidate_index = candidate_indices[optimal_point_index];
   std::vector<double> final_column_vector;
   double final_self_value;
-  if(for_coeffs) {
+  if (for_coeffs) {
 
     // Grow the kernel cache in this case.
     kernel_matrix_columns_[final_candidate_index].resize(dataset_->n_cols());
@@ -265,13 +327,13 @@ void SparseGreedyGpr::ChooseRandomSubset_(
 
   // Copy
   subset_out->resize(inactive_set.size());
-  for(int i = 0; i < inactive_set.size(); i++) {
+  for (int i = 0; i < inactive_set.size(); i++) {
     (*subset_out)[i] = inactive_set[i];
   }
 
   // Then shuffle, and truncate.
-  if(inactive_set.size() > subset_size) {
-    for(int i = subset_out->size() - 1; i >= 1; i--) {
+  if (inactive_set.size() > subset_size) {
+    for (int i = subset_out->size() - 1; i >= 1; i--) {
 
       // Pick a random index between 0 and i, inclusive.
       int random_index = math::RandInt(0, i);
@@ -321,7 +383,7 @@ void SparseGreedyGpr::Compute(
   // sets).
   std::vector<int> inactive_indices(dataset_->n_cols());
   std::vector<int> inactive_indices_for_error(dataset_-> n_cols());
-  for(int i = 0; i < dataset_->n_cols(); i++) {
+  for (int i = 0; i < dataset_->n_cols(); i++) {
     inactive_indices[i] = i;
     inactive_indices_for_error[i] = i;
   }
@@ -351,9 +413,9 @@ void SparseGreedyGpr::Compute(
     dictionary_.inactive_indices(&inactive_indices);
     dictionary_for_error_.inactive_indices(&inactive_indices_for_error_);
   }
-  while(Done_(
-          model_out->frobenius_norm_targets(), noise_level_in, precision_in,
-          optimum_value, optimum_value_for_error));
+  while (Done_(
+           model_out->frobenius_norm_targets(), noise_level_in, precision_in,
+           optimum_value, optimum_value_for_error));
 
   // Using the final model, compute the coefficients.
   model_out->FinalizeModel();