Use more iter_rows and select_rows (#107)

* use iter_rows on GaussianProcess::ker_mat

* use iter_rows on SVM::ker_mat

* use iter_rows on SVM::train

* use select_rows in DBSCAN::expand_cluster

* use iter_rows in NaiveBayes::update_param

* use iter_rows_mut in GaussianMixtureModel::update_param

* use iter_rows in NaiveBayes::get_class

* remove unecessary &[..]

* use select_rows in SVM::train

* use matrix scalar product in tests

Author: tafia

src/learning/dbscan.rs

@@ -178,11 +178,8 @@ impl DBSCAN {
             let visited = self._visited[*data_point_idx];
             if !visited {
                 self._visited[*data_point_idx] = true;
-                let sub_neighbours =
-                    self.region_query(&inputs.data()[data_point_idx * inputs.cols()..(data_point_idx +
-                                                                                   1) *
-                                                                                  inputs.cols()],
-                                      inputs);
+                let sub_neighbours = self.region_query(inputs.select_rows(&[*data_point_idx]).data(),
+                                                       inputs);
 
                 if sub_neighbours.len() >= self.min_points {
                     self.expand_cluster(inputs, *data_point_idx, sub_neighbours, cluster);

src/learning/gmm.rs

@@ -279,9 +279,9 @@ impl GaussianMixtureModel {
 
         let mut new_means = membership_weights.transpose() * inputs;
 
-        for (idx, mean) in new_means.mut_data().chunks_mut(d).enumerate() {
-            for m in mean {
-                *m /= sum_weights[idx];
+        for (mean, w) in new_means.iter_rows_mut().zip(sum_weights.data().iter()) {
+            for m in mean.iter_mut() {
+                *m /= *w;
             }
         }
 

src/learning/gp.rs

@@ -122,19 +122,14 @@ impl<T: Kernel, U: MeanFunc> GaussianProcess<T, U> {
     /// Construct a kernel matrix
     fn ker_mat(&self, m1: &Matrix<f64>, m2: &Matrix<f64>) -> Matrix<f64> {
         assert_eq!(m1.cols(), m2.cols());
-        let cols = m1.cols();
 
         let dim1 = m1.rows();
         let dim2 = m2.rows();
 
         let mut ker_data = Vec::with_capacity(dim1 * dim2);
-
-        for i in 0..dim1 {
-            for j in 0..dim2 {
-                ker_data.push(self.ker.kernel(&m1.data()[i * cols..(i + 1) * cols],
-                                              &m2.data()[j * cols..(j + 1) * cols]));
-            }
-        }
+        ker_data.extend(
+            m1.iter_rows().flat_map(|row1| m2.iter_rows()
+                                    .map(move |row2| self.ker.kernel(row1, row2))));
 
         Matrix::new(dim1, dim2, ker_data)
     }

src/learning/naive_bayes.rs

@@ -151,7 +151,7 @@ impl<T: Distribution> NaiveBayes<T> {
         self.class_counts = vec![0; class_count];
         let mut class_data = vec![Vec::new(); class_count];
 
-        for (idx, row) in targets.data().chunks(class_count).enumerate() {
+        for (idx, row) in targets.iter_rows().enumerate() {
             // Find the class of this input
             let class = NaiveBayes::<T>::find_class(row);
 
@@ -195,15 +195,13 @@ impl<T: Distribution> NaiveBayes<T> {
     }
 
     fn get_classes(log_probs: Matrix<f64>) -> Vec<usize> {
-        let class_count = log_probs.cols();
         let mut data_classes = Vec::with_capacity(log_probs.rows());
 
-        // Argmax each class log-probability per input
-        for row in log_probs.data().chunks(class_count) {
+        data_classes.extend(log_probs.iter_rows().map(|row| {
+            // Argmax each class log-probability per input
             let (class, _) = utils::argmax(row);
-
-            data_classes.push(class);
-        }
+            class
+        }));
 
         data_classes
     }

src/learning/optim/fmincg.rs

@@ -115,7 +115,7 @@ impl<M: Optimizable> OptimAlgorithm<M> for ConjugateGD {
 
             x = x + &s * z1;
 
-            let cost = model.compute_grad(&x.data()[..], inputs, targets);
+            let cost = model.compute_grad(x.data(), inputs, targets);
             f2 = cost.0;
             df2 = Vector::new(cost.1);
 
@@ -167,7 +167,7 @@ impl<M: Optimizable> OptimAlgorithm<M> for ConjugateGD {
 
                     z1 += z2;
                     x = x + &s * z2;
-                    let cost_grad = model.compute_grad(&x.data()[..], inputs, targets);
+                    let cost_grad = model.compute_grad(x.data(), inputs, targets);
                     f2 = cost_grad.0;
                     df2 = Vector::new(cost_grad.1);
 
@@ -215,7 +215,7 @@ impl<M: Optimizable> OptimAlgorithm<M> for ConjugateGD {
                 z1 += z2;
                 x = x + &s * z2;
 
-                let cost_grad = model.compute_grad(&x.data()[..], inputs, targets);
+                let cost_grad = model.compute_grad(x.data(), inputs, targets);
                 f2 = cost_grad.0;
                 df2 = Vector::new(cost_grad.1);
 

src/learning/optim/grad_desc.rs

@@ -80,7 +80,7 @@ impl<M: Optimizable> OptimAlgorithm<M> for GradientDesc {
 
         for _ in 0..self.iters {
             // Compute the cost and gradient for the current parameters
-            let (cost, grad) = model.compute_grad(&optimizing_val.data()[..], inputs, targets);
+            let (cost, grad) = model.compute_grad(optimizing_val.data(), inputs, targets);
 
             // Early stopping
             if (start_iter_cost - cost).abs() < LEARNING_EPS {

src/learning/svm.rs

@@ -103,19 +103,15 @@ impl<K: Kernel> SVM<K> {
     /// Construct a kernel matrix
     fn ker_mat(&self, m1: &Matrix<f64>, m2: &Matrix<f64>) -> Matrix<f64> {
         assert_eq!(m1.cols(), m2.cols());
-        let cols = m1.cols();
 
         let dim1 = m1.rows();
         let dim2 = m2.rows();
 
         let mut ker_data = Vec::with_capacity(dim1 * dim2);
 
-        for i in 0..dim1 {
-            for j in 0..dim2 {
-                ker_data.push(self.ker.kernel(&m1.data()[i * cols..(i + 1) * cols],
-                                              &m2.data()[j * cols..(j + 1) * cols]));
-            }
-        }
+        ker_data.extend(
+            m1.iter_rows().flat_map(|row1| m2.iter_rows()
+                                    .map(move |row2| self.ker.kernel(row1, row2))));
 
         Matrix::new(dim1, dim2, ker_data)
     }
@@ -150,17 +146,13 @@ impl<K: Kernel> SupModel<Matrix<f64>, Vector<f64>> for SVM<K> {
 
         let ones = Matrix::<f64>::ones(inputs.rows(), 1);
         let full_inputs = ones.hcat(inputs);
-        let m = full_inputs.cols();
 
         for t in 0..self.optim_iters {
             let i = rng.gen_range(0, n);
-            let mut sum = 0f64;
-            for j in 0..n {
-                sum += alpha[j] * targets[j] *
-                       self.ker.kernel(&full_inputs.data()[i * m..(i + 1) * m],
-                                       &full_inputs.data()[j * m..(j + 1) * m]);
-            }
-            sum *= targets[i] / (self.lambda * (t as f64));
+            let row_i = full_inputs.select_rows(&[i]);
+            let sum =  full_inputs.iter_rows()
+                .fold(0f64, |sum, row| sum + self.ker.kernel(row_i.data(), row)) * 
+                targets[i] / (self.lambda * (t as f64));
 
             if sum < 1f64 {
                 alpha[i] += 1f64;

src/learning/toolkit/regularization.rs

@@ -141,11 +141,8 @@ mod tests {
 
         assert!((a - (input_mat.norm() / 12f64)) < 1e-18);
 
-        let true_grad = input_mat.data()
-            .iter()
-            .map(|x| x / 6f64)
-            .collect::<Vec<_>>();
-        for eps in (b - Matrix::new(3, 4, true_grad)).into_vec() {
+        let true_grad = &input_mat / 6f64;
+        for eps in (b - true_grad).into_vec() {
             assert!(eps < 1e-18);
         }
     }
@@ -162,16 +159,14 @@ mod tests {
 
         assert!(a - ((input_mat.norm() / 24f64) + (42f64 / 12f64)) < 1e-18);
 
-        let l1_true_grad = vec![-1., -1., -1., 1., 1., 1., 1., 1., 1., 1., 1., 1.]
+        let l1_true_grad = Matrix::new(3, 4,
+            vec![-1., -1., -1., 1., 1., 1., 1., 1., 1., 1., 1., 1.]
             .into_iter()
             .map(|x| x / 12f64)
-            .collect::<Vec<_>>();
-        let l2_true_grad = input_mat.data()
-            .iter()
-            .map(|x| x / 12f64)
-            .collect::<Vec<_>>();
+            .collect::<Vec<_>>());
+        let l2_true_grad = &input_mat / 12f64;
 
-        for eps in (b - Matrix::new(3, 4, l1_true_grad) - Matrix::new(3, 4, l2_true_grad))
+        for eps in (b - l1_true_grad - l2_true_grad)
             .into_vec() {
             // Slightly lower boundary than others - more numerical error as more ops.
             assert!(eps < 1e-12);