diff --git a/examples/k-means_generating_cluster.rs b/examples/k-means_generating_cluster.rs
index 7b8c1368..0c05f1d5 100644
--- a/examples/k-means_generating_cluster.rs
+++ b/examples/k-means_generating_cluster.rs
@@ -58,16 +58,17 @@ fn main() {
     // Create a new model with 2 clusters
     let mut model = KMeansClassifier::new(2);
 
-    println!("Training the model...");
     // Train the model
-    model.train(&samples);
+    println!("Training the model...");
+    // Our train function returns a Result<(), E>
+    model.train(&samples).unwrap();
 
     let centroids = model.centroids().as_ref().unwrap();
     println!("Model Centroids:\n{:.3}", centroids);
 
     // Predict the classes and partition into
     println!("Classifying the samples...");
-    let classes = model.predict(&samples);
+    let classes = model.predict(&samples).unwrap();
     let (first, second): (Vec<usize>, Vec<usize>) = classes.data().iter().partition(|&x| *x == 0);
 
     println!("Samples closest to first centroid: {}", first.len());
diff --git a/examples/nnet-and_gate.rs b/examples/nnet-and_gate.rs
index bf319bf1..b4c6b156 100644
--- a/examples/nnet-and_gate.rs
+++ b/examples/nnet-and_gate.rs
@@ -44,7 +44,8 @@ fn main() {
     let mut model = NeuralNet::new(layers, criterion, StochasticGD::default());
 
     println!("Training...");
-    model.train(&inputs, &targets);
+    // Our train function returns a Result<(), E>
+    model.train(&inputs, &targets).unwrap();
 
     let test_cases = vec![
         0.0, 0.0,
@@ -59,7 +60,7 @@ fn main() {
         0.0,
         ];
     let test_inputs = Matrix::new(test_cases.len() / 2, 2, test_cases);
-    let res = model.predict(&test_inputs);
+    let res = model.predict(&test_inputs).unwrap();
 
     println!("Evaluation...");
     let mut hits = 0;
diff --git a/examples/svm-sign_learner.rs b/examples/svm-sign_learner.rs
index f1a27203..b1a16907 100644
--- a/examples/svm-sign_learner.rs
+++ b/examples/svm-sign_learner.rs
@@ -30,7 +30,8 @@ fn main() {
 
     // Trainee
     let mut svm_mod = SVM::new(HyperTan::new(100., 0.), 0.3);
-    svm_mod.train(&inputs, &targets);
+    // Our train function returns a Result<(), E>
+    svm_mod.train(&inputs, &targets).unwrap();
 
     println!("Evaluation...");
     let mut hits = 0;
@@ -41,7 +42,7 @@ fn main() {
     for n in (-1000..1000).filter(|&x| x % 100 == 0) {
         let nf = n as f64;
         let input = Matrix::new(1, 1, vec![nf]);
-        let out = svm_mod.predict(&input);
+        let out = svm_mod.predict(&input).unwrap();
         let res = if out[0] * nf > 0. {
             hits += 1;
             true
diff --git a/src/learning/dbscan.rs b/src/learning/dbscan.rs
index 1d82ba0c..4dbfefb6 100644
--- a/src/learning/dbscan.rs
+++ b/src/learning/dbscan.rs
@@ -31,12 +31,13 @@
 //!                                     -2.2, 3.1]);
 //!
 //! let mut model = DBSCAN::new(0.5, 2);
-//! model.train(&inputs);
+//! model.train(&inputs).unwrap();
 //!
 //! let clustering = model.clusters().unwrap();
 //! ```
 
-use learning::UnSupModel;
+use learning::{LearningResult, UnSupModel};
+use learning::error::{Error, ErrorKind};
 
 use linalg::{Matrix, Vector};
 use rulinalg::utils;
@@ -75,7 +76,7 @@ impl Default for DBSCAN {
 
 impl UnSupModel<Matrix<f64>, Vector<Option<usize>>> for DBSCAN {
     /// Train the classifier using input data.
-    fn train(&mut self, inputs: &Matrix<f64>) {
+    fn train(&mut self, inputs: &Matrix<f64>) -> LearningResult<()> {
         self.init_params(inputs.rows());
         let mut cluster = 0;
 
@@ -95,11 +96,13 @@ impl UnSupModel<Matrix<f64>, Vector<Option<usize>>> for DBSCAN {
         }
 
         if self.predictive {
-            self._cluster_data = Some(inputs.clone())
+            self._cluster_data = Some(inputs.clone());
         }
+
+        Ok(())
     }
 
-    fn predict(&self, inputs: &Matrix<f64>) -> Vector<Option<usize>> {
+    fn predict(&self, inputs: &Matrix<f64>) -> LearningResult<Vector<Option<usize>>> {
         if self.predictive {
             if let (&Some(ref cluster_data), &Some(ref clusters)) = (&self._cluster_data,
                                                                      &self.clusters) {
@@ -122,12 +125,13 @@ impl UnSupModel<Matrix<f64>, Vector<Option<usize>>> for DBSCAN {
                     }
                 }
 
-                Vector::new(classes)
+                Ok(Vector::new(classes))
             } else {
-                panic!("The model has not been trained.");
+                Err(Error::new_untrained())
             }
         } else {
-            panic!("Model must be set to predictive. Use `self.set_predictive(true)`.");
+            Err(Error::new(ErrorKind::InvalidState,
+                           "Model must be set to predictive. Use `self.set_predictive(true)`."))
         }
     }
 }
diff --git a/src/learning/error.rs b/src/learning/error.rs
index 135c23a0..955ef684 100644
--- a/src/learning/error.rs
+++ b/src/learning/error.rs
@@ -25,6 +25,8 @@ pub enum ErrorKind {
     InvalidData,
     /// The action could not be carried out as the model was in an invalid state.
     InvalidState,
+    /// The model has not been trained
+    UntrainedModel
 }
 
 impl Error {
@@ -38,6 +40,13 @@ impl Error {
         }
     }
 
+    /// Returns a new error for an untrained model
+    ///
+    /// This function is unstable and may be removed with changes to the API.
+    pub fn new_untrained() -> Error {
+        Error::new(ErrorKind::UntrainedModel, "The model has not been trained.")
+    }
+
     /// Get the kind of this `Error`.
     pub fn kind(&self) -> &ErrorKind {
         &self.kind
diff --git a/src/learning/glm.rs b/src/learning/glm.rs
index 874985ec..5321fd18 100644
--- a/src/learning/glm.rs
+++ b/src/learning/glm.rs
@@ -25,11 +25,11 @@
 //! let mut log_mod = GenLinearModel::new(Bernoulli);
 //!
 //! // Train the model
-//! log_mod.train(&inputs, &targets);
+//! log_mod.train(&inputs, &targets).unwrap();
 //!
 //! // Now we'll predict a new point
 //! let new_point = Matrix::new(1,1,vec![10.]);
-//! let output = log_mod.predict(&new_point);
+//! let output = log_mod.predict(&new_point).unwrap();
 //!
 //! // Hopefully we classified our new point correctly!
 //! assert!(output[0] > 0.5, "Our classifier isn't very good!");
@@ -38,7 +38,8 @@
 use linalg::Vector;
 use linalg::Matrix;
 
-use learning::SupModel;
+use learning::{LearningResult, SupModel};
+use learning::error::{Error, ErrorKind};
 
 /// The Generalized Linear Model
 ///
@@ -78,22 +79,24 @@ impl<C: Criterion> GenLinearModel<C> {
 /// The model is trained using Iteratively Re-weighted Least Squares.
 impl<C: Criterion> SupModel<Matrix<f64>, Vector<f64>> for GenLinearModel<C> {
     /// Predict output from inputs.
-    fn predict(&self, inputs: &Matrix<f64>) -> Vector<f64> {
+    fn predict(&self, inputs: &Matrix<f64>) -> LearningResult<Vector<f64>> {
         if let Some(ref v) = self.parameters {
             let ones = Matrix::<f64>::ones(inputs.rows(), 1);
             let full_inputs = ones.hcat(inputs);
-            self.criterion.apply_link_inv(full_inputs * v)
+            Ok(self.criterion.apply_link_inv(full_inputs * v))
         } else {
-            panic!("The model has not been trained.");
+            Err(Error::new_untrained())
         }
     }
 
     /// Train the model using inputs and targets.
-    fn train(&mut self, inputs: &Matrix<f64>, targets: &Vector<f64>) {
+    fn train(&mut self, inputs: &Matrix<f64>, targets: &Vector<f64>) -> LearningResult<()> {
         let n = inputs.rows();
 
-        assert!(n == targets.size(),
-                "Training data do not have the same dimensions.");
+        if n != targets.size() {
+            return Err(Error::new(ErrorKind::InvalidData,
+                                  "Training data do not have the same dimensions"));
+        }
 
         // Construct initial estimate for mu
         let mut mu = Vector::new(self.criterion.initialize_mu(targets.data()));
@@ -132,6 +135,7 @@ impl<C: Criterion> SupModel<Matrix<f64>, Vector<f64>> for GenLinearModel<C> {
         }
 
         self.parameters = Some(beta);
+        Ok(())
     }
 }
 
diff --git a/src/learning/gmm.rs b/src/learning/gmm.rs
index b9906477..3e1eae76 100644
--- a/src/learning/gmm.rs
+++ b/src/learning/gmm.rs
@@ -18,14 +18,14 @@
 //! model.cov_option = CovOption::Diagonal;
 //!
 //! // Where inputs is a Matrix with features in columns.
-//! model.train(&inputs);
+//! model.train(&inputs).unwrap();
 //!
 //! // Print the means and covariances of the GMM
 //! println!("{:?}", model.means());
 //! println!("{:?}", model.covariances());
 //!
 //! // Where test_inputs is a Matrix with features in columns.
-//! let post_probs = model.predict(&test_inputs);
+//! let post_probs = model.predict(&test_inputs).unwrap();
 //!
 //! // Probabilities that each point comes from each Gaussian.
 //! println!("{:?}", post_probs.data());
@@ -34,8 +34,9 @@
 use linalg::{Matrix, MatrixSlice, Vector};
 use rulinalg::utils;
 
-use learning::UnSupModel;
+use learning::{LearningResult, UnSupModel};
 use learning::toolkit::rand_utils;
+use learning::error::{Error, ErrorKind};
 
 /// Covariance options for GMMs.
 ///
@@ -68,7 +69,7 @@ pub struct GaussianMixtureModel {
 
 impl UnSupModel<Matrix<f64>, Matrix<f64>> for GaussianMixtureModel {
     /// Train the model using inputs.
-    fn train(&mut self, inputs: &Matrix<f64>) {
+    fn train(&mut self, inputs: &Matrix<f64>) -> LearningResult<()> {
         // Initialization:
         let k = self.comp_count;
 
@@ -98,14 +99,16 @@ impl UnSupModel<Matrix<f64>, Matrix<f64>> for GaussianMixtureModel {
 
             self.update_params(inputs, weights);
         }
+
+        Ok(())
     }
 
     /// Predict output from inputs.
-    fn predict(&self, inputs: &Matrix<f64>) -> Matrix<f64> {
+    fn predict(&self, inputs: &Matrix<f64>) -> LearningResult<Matrix<f64>> {
         if let (&Some(_), &Some(_)) = (&self.model_means, &self.model_covars) {
-            self.membership_weights(inputs).0
+            Ok(self.membership_weights(inputs).0)
         } else {
-            panic!("Model has not been trained.");
+            Err(Error::new_untrained())
         }
 
     }
@@ -148,32 +151,33 @@ impl GaussianMixtureModel {
     ///
     /// let mix_weights = Vector::new(vec![0.25, 0.25, 0.5]);
     ///
-    /// let _ = GaussianMixtureModel::with_weights(3, mix_weights);
+    /// let gmm = GaussianMixtureModel::with_weights(3, mix_weights).unwrap();
     /// ```
     ///
-    /// # Panics
+    /// # Failures
     ///
-    /// Panics if either of the following conditions are met:
+    /// Fails if either of the following conditions are met:
     ///
     /// - Mixture weights do not have length k.
     /// - Mixture weights have a negative entry.
-    pub fn with_weights(k: usize, mixture_weights: Vector<f64>) -> GaussianMixtureModel {
-        assert!(mixture_weights.size() == k,
-                "Mixture weights must have length k.");
-        assert!(!mixture_weights.data().iter().any(|&x| x < 0f64),
-                "Mixture weights must have only non-negative entries.");
-
-        let sum = mixture_weights.sum();
-        let normalized_weights = mixture_weights / sum;
-
-        GaussianMixtureModel {
-            comp_count: k,
-            mix_weights: normalized_weights,
-            model_means: None,
-            model_covars: None,
-            log_lik: 0f64,
-            max_iters: 100,
-            cov_option: CovOption::Full,
+    pub fn with_weights(k: usize, mixture_weights: Vector<f64>) -> LearningResult<GaussianMixtureModel> {
+        if mixture_weights.size() != k {
+            Err(Error::new(ErrorKind::InvalidParameters, "Mixture weights must have length k."))
+        } else if mixture_weights.data().iter().any(|&x| x < 0f64) {
+            Err(Error::new(ErrorKind::InvalidParameters, "Mixture weights must have only non-negative entries.")) 
+        } else {
+            let sum = mixture_weights.sum();
+            let normalized_weights = mixture_weights / sum;
+
+            Ok(GaussianMixtureModel {
+                comp_count: k,
+                mix_weights: normalized_weights,
+                model_means: None,
+                model_covars: None,
+                log_lik: 0f64,
+                max_iters: 100,
+                cov_option: CovOption::Full,
+            })
         }
     }
 
@@ -292,7 +296,7 @@ impl GaussianMixtureModel {
 
             for i in 0..n {
                 let inputs_i = MatrixSlice::from_matrix(inputs, [i, 0], 1, d);
-                let diff = inputs_i - new_means_k; 
+                let diff = inputs_i - new_means_k;
                 cov_mat += self.compute_cov(diff, membership_weights[[i, k]]);
             }
             new_covs.push(cov_mat / sum_weights[k]);
@@ -332,16 +336,16 @@ mod tests {
     }
 
     #[test]
-    #[should_panic]
     fn test_negative_mixtures() {
         let mix_weights = Vector::new(vec![-0.25, 0.75, 0.5]);
-        let _ = GaussianMixtureModel::with_weights(3, mix_weights);
+        let gmm_res = GaussianMixtureModel::with_weights(3, mix_weights);
+        assert!(gmm_res.is_err());
     }
 
     #[test]
-    #[should_panic]
     fn test_wrong_length_mixtures() {
         let mix_weights = Vector::new(vec![0.1, 0.25, 0.75, 0.5]);
-        let _ = GaussianMixtureModel::with_weights(3, mix_weights);
+        let gmm_res = GaussianMixtureModel::with_weights(3, mix_weights);
+        assert!(gmm_res.is_err());
     }
 }
diff --git a/src/learning/gp.rs b/src/learning/gp.rs
index 12ec0b33..b1f0eefb 100644
--- a/src/learning/gp.rs
+++ b/src/learning/gp.rs
@@ -16,21 +16,24 @@
 //! let train_data = Matrix::new(10,1,vec![0.,1.,2.,3.,4.,5.,6.,7.,8.,9.]);
 //! let target = Vector::new(vec![0.,1.,2.,3.,4.,4.,3.,2.,1.,0.]);
 //!
-//! gaussp.train(&train_data, &target);
+//! gaussp.train(&train_data, &target).unwrap();
 //!
 //! let test_data = Matrix::new(5,1,vec![2.3,4.4,5.1,6.2,7.1]);
 //!
-//! let outputs = gaussp.predict(&test_data);
+//! let outputs = gaussp.predict(&test_data).unwrap();
 //! ```
 //! Alternatively one could use `gaussp.get_posterior()` which would return both
 //! the predictive mean and covariance. However, this is likely to change in
 //! a future release.
 
 use learning::toolkit::kernel::{Kernel, SquaredExp};
-use learning::SupModel;
+use learning::{LearningResult, SupModel};
+use learning::error::{Error, ErrorKind};
+
 use linalg::Matrix;
 use linalg::Vector;
 
+
 /// Trait for GP mean functions.
 pub trait MeanFunc {
     /// Compute the mean function applied elementwise to a matrix.
@@ -120,46 +123,49 @@ 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 dim1 = m1.rows();
-        let dim2 = m2.rows();
-
-        let mut ker_data = Vec::with_capacity(dim1 * dim2);
-        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)
+    fn ker_mat(&self, m1: &Matrix<f64>, m2: &Matrix<f64>) -> LearningResult<Matrix<f64>> {
+        if m1.cols() != m2.cols() {
+            Err(Error::new(ErrorKind::InvalidState,
+                           "Inputs to kernel matrices have different column counts."))
+        } else {
+            let dim1 = m1.rows();
+            let dim2 = m2.rows();
+
+            let mut ker_data = Vec::with_capacity(dim1 * dim2);
+            ker_data.extend(m1.iter_rows().flat_map(|row1| {
+                m2.iter_rows()
+                    .map(move |row2| self.ker.kernel(row1, row2))
+            }));
+
+            Ok(Matrix::new(dim1, dim2, ker_data))
+        }
     }
 }
 
 impl<T: Kernel, U: MeanFunc> SupModel<Matrix<f64>, Vector<f64>> for GaussianProcess<T, U> {
     /// Predict output from inputs.
-    fn predict(&self, inputs: &Matrix<f64>) -> Vector<f64> {
+    fn predict(&self, inputs: &Matrix<f64>) -> LearningResult<Vector<f64>> {
 
         // Messy referencing for succint syntax
         if let (&Some(ref alpha), &Some(ref t_data)) = (&self.alpha, &self.train_data) {
             let mean = self.mean.func(inputs.clone());
-
-            let post_mean = self.ker_mat(inputs, t_data) * alpha;
-
-            return mean + post_mean;
-
+            let post_mean = try!(self.ker_mat(inputs, t_data)) * alpha;
+            Ok(mean + post_mean)
+        } else {
+            Err(Error::new(ErrorKind::UntrainedModel, "The model has not been trained."))
         }
-
-        panic!("The model has not been trained.");
     }
 
     /// Train the model using data and outputs.
-    fn train(&mut self, inputs: &Matrix<f64>, targets: &Vector<f64>) {
+    fn train(&mut self, inputs: &Matrix<f64>, targets: &Vector<f64>) -> LearningResult<()> {
         let noise_mat = Matrix::identity(inputs.rows()) * self.noise;
 
-        let ker_mat = self.ker_mat(inputs, inputs);
+        let ker_mat = self.ker_mat(inputs, inputs).unwrap();
 
-        let train_mat =
-            (ker_mat + noise_mat).cholesky().expect("Could not compute Cholesky decomposition.");
+        let train_mat = try!((ker_mat + noise_mat).cholesky().map_err(|_| {
+            Error::new(ErrorKind::InvalidState,
+                       "Could not compute Cholesky decomposition.")
+        }));
 
         let x = train_mat.solve_l_triangular(targets - self.mean.func(inputs.clone())).unwrap();
         let alpha = train_mat.transpose().solve_u_triangular(x).unwrap();
@@ -167,6 +173,8 @@ impl<T: Kernel, U: MeanFunc> SupModel<Matrix<f64>, Vector<f64>> for GaussianProc
         self.train_mat = Some(train_mat);
         self.train_data = Some(inputs.clone());
         self.alpha = Some(alpha);
+
+        Ok(())
     }
 }
 
@@ -176,15 +184,17 @@ impl<T: Kernel, U: MeanFunc> GaussianProcess<T, U> {
     /// Requires the model to be trained first.
     ///
     /// Outputs the posterior mean and covariance matrix.
-    pub fn get_posterior(&self, inputs: &Matrix<f64>) -> (Vector<f64>, Matrix<f64>) {
+    pub fn get_posterior(&self,
+                         inputs: &Matrix<f64>)
+                         -> LearningResult<(Vector<f64>, Matrix<f64>)> {
         if let (&Some(ref t_mat), &Some(ref alpha), &Some(ref t_data)) = (&self.train_mat,
                                                                           &self.alpha,
                                                                           &self.train_data) {
             let mean = self.mean.func(inputs.clone());
 
-            let post_mean = mean + self.ker_mat(inputs, t_data) * alpha;
+            let post_mean = mean + try!(self.ker_mat(inputs, t_data)) * alpha;
 
-            let test_mat = self.ker_mat(inputs, t_data);
+            let test_mat = try!(self.ker_mat(inputs, t_data));
             let mut var_data = Vec::with_capacity(inputs.rows() * inputs.cols());
             for row in test_mat.iter_rows() {
                 let test_point = Vector::new(row.to_vec());
@@ -193,11 +203,11 @@ impl<T: Kernel, U: MeanFunc> GaussianProcess<T, U> {
 
             let v_mat = Matrix::new(test_mat.rows(), test_mat.cols(), var_data);
 
-            let post_var = self.ker_mat(inputs, inputs) - &v_mat * v_mat.transpose();
+            let post_var = try!(self.ker_mat(inputs, inputs)) - &v_mat * v_mat.transpose();
 
-            return (post_mean, post_var);
+            Ok((post_mean, post_var))
+        } else {
+            Err(Error::new_untrained())
         }
-
-        panic!("The model has not been trained.");
     }
 }
diff --git a/src/learning/k_means.rs b/src/learning/k_means.rs
index a1bf38e2..6c8ad3cc 100644
--- a/src/learning/k_means.rs
+++ b/src/learning/k_means.rs
@@ -16,10 +16,10 @@
 //! let mut model = KMeansClassifier::new(2);
 //!
 //! // Where inputs is a Matrix with features in columns.
-//! model.train(&inputs);
+//! model.train(&inputs).unwrap();
 //!
 //! // Where test_inputs is a Matrix with features in columns.
-//! let a = model.predict(&test_inputs);
+//! let a = model.predict(&test_inputs).unwrap();
 //! ```
 //!
 //! Additionally you can control the initialization
@@ -45,7 +45,7 @@
 use linalg::BaseSlice;
 use linalg::{Matrix, MatrixSlice, Axes};
 use linalg::Vector;
-use learning::UnSupModel;
+use learning::{LearningResult, UnSupModel};
 use learning::error::{Error, ErrorKind};
 
 use rand::{Rng, thread_rng};
@@ -82,24 +82,22 @@ impl<InitAlg: Initializer> UnSupModel<Matrix<f64>, Vector<usize>> for KMeansClas
     /// Predict classes from data.
     ///
     /// Model must be trained.
-    fn predict(&self, inputs: &Matrix<f64>) -> Vector<usize> {
+    fn predict(&self, inputs: &Matrix<f64>) -> LearningResult<Vector<usize>> {
         if let Some(ref centroids) = self.centroids {
-            return KMeansClassifier::<InitAlg>::find_closest_centroids(centroids.as_slice(),
-                                                                       inputs)
-                .0;
+            Ok(KMeansClassifier::<InitAlg>::find_closest_centroids(centroids.as_slice(), inputs).0)
         } else {
-            panic!("Model has not been trained.");
+            Err(Error::new_untrained())
         }
     }
 
     /// Train the classifier using input data.
-    fn train(&mut self, inputs: &Matrix<f64>) {
-        self.init_centroids(inputs).expect("Could not initialize centroids.");
+    fn train(&mut self, inputs: &Matrix<f64>) -> LearningResult<()> {
+        try!(self.init_centroids(inputs));
         let mut cost = 0.0;
         let eps = 1e-14;
 
         for _i in 0..self.iters {
-            let (idx, distances) = self.get_closest_centroids(inputs);
+            let (idx, distances) = try!(self.get_closest_centroids(inputs));
             self.update_centroids(inputs, idx);
 
             let cost_i = distances.sum();
@@ -109,6 +107,8 @@ impl<InitAlg: Initializer> UnSupModel<Matrix<f64>, Vector<usize>> for KMeansClas
 
             cost = cost_i;
         }
+
+        Ok(())
     }
 }
 
@@ -185,7 +185,7 @@ impl<InitAlg: Initializer> KMeansClassifier<InitAlg> {
     /// Initialize the centroids.
     ///
     /// Used internally within model.
-    fn init_centroids(&mut self, inputs: &Matrix<f64>) -> Result<(), Error> {
+    fn init_centroids(&mut self, inputs: &Matrix<f64>) -> LearningResult<()> {
         if self.k > inputs.rows() {
             Err(Error::new(ErrorKind::InvalidData,
                            format!("Number of clusters ({0}) exceeds number of data points \
@@ -194,12 +194,17 @@ impl<InitAlg: Initializer> KMeansClassifier<InitAlg> {
                                    inputs.rows())))
         } else {
             let centroids = try!(self.init_algorithm.init_centroids(self.k, inputs));
-            assert!(centroids.rows() == self.k,
-                    "Initial centroids must have exactly k rows.");
-            assert!(centroids.cols() == inputs.cols(),
-                    "Initial centroids must have the same column count as inputs.");
-            self.centroids = Some(centroids);
-            Ok(())
+
+            if centroids.rows() != self.k {
+                Err(Error::new(ErrorKind::InvalidState,
+                                    "Initial centroids must have exactly k rows."))
+            } else if centroids.cols() != inputs.cols() {
+                Err(Error::new(ErrorKind::InvalidState,
+                                    "Initial centroids must have the same column count as inputs."))
+            } else {
+                self.centroids = Some(centroids);
+                Ok(())
+            }
         }
 
     }
@@ -223,11 +228,14 @@ impl<InitAlg: Initializer> KMeansClassifier<InitAlg> {
         self.centroids = Some(Matrix::new(self.k, inputs.cols(), new_centroids));
     }
 
-    fn get_closest_centroids(&self, inputs: &Matrix<f64>) -> (Vector<usize>, Vector<f64>) {
+    fn get_closest_centroids(&self,
+                             inputs: &Matrix<f64>)
+                             -> LearningResult<(Vector<usize>, Vector<f64>)> {
         if let Some(ref c) = self.centroids {
-            return KMeansClassifier::<InitAlg>::find_closest_centroids(c.as_slice(), inputs);
+            Ok(KMeansClassifier::<InitAlg>::find_closest_centroids(c.as_slice(), inputs))
         } else {
-            panic!("Centroids not correctly initialized.");
+            Err(Error::new(ErrorKind::InvalidState,
+                           "Centroids not correctly initialized."))
         }
     }
 
@@ -261,7 +269,7 @@ pub trait Initializer: Debug {
     /// Initialize the centroids for the initial state of the K-Means model.
     ///
     /// The `Matrix` returned must have `k` rows and the same column count as `inputs`.
-    fn init_centroids(&self, k: usize, inputs: &Matrix<f64>) -> Result<Matrix<f64>, Error>;
+    fn init_centroids(&self, k: usize, inputs: &Matrix<f64>) -> LearningResult<Matrix<f64>>;
 }
 
 /// The Forgy initialization scheme.
@@ -269,7 +277,7 @@ pub trait Initializer: Debug {
 pub struct Forgy;
 
 impl Initializer for Forgy {
-    fn init_centroids(&self, k: usize, inputs: &Matrix<f64>) -> Result<Matrix<f64>, Error> {
+    fn init_centroids(&self, k: usize, inputs: &Matrix<f64>) -> LearningResult<Matrix<f64>> {
         let mut random_choices = Vec::with_capacity(k);
         let mut rng = thread_rng();
         while random_choices.len() < k {
@@ -289,14 +297,16 @@ impl Initializer for Forgy {
 pub struct RandomPartition;
 
 impl Initializer for RandomPartition {
-    fn init_centroids(&self, k: usize, inputs: &Matrix<f64>) -> Result<Matrix<f64>, Error> {
+    fn init_centroids(&self, k: usize, inputs: &Matrix<f64>) -> LearningResult<Matrix<f64>> {
 
         // Populate so we have something in each class.
         let mut random_assignments = (0..k).map(|i| vec![i]).collect::<Vec<Vec<usize>>>();
         let mut rng = thread_rng();
         for i in k..inputs.rows() {
             let idx = rng.gen_range(0, k);
-            unsafe { random_assignments.get_unchecked_mut(idx).push(i); }
+            unsafe {
+                random_assignments.get_unchecked_mut(idx).push(i);
+            }
         }
 
         let mut init_centroids = Vec::with_capacity(k * inputs.cols());
@@ -315,14 +325,14 @@ impl Initializer for RandomPartition {
 pub struct KPlusPlus;
 
 impl Initializer for KPlusPlus {
-    fn init_centroids(&self, k: usize, inputs: &Matrix<f64>) -> Result<Matrix<f64>, Error> {
+    fn init_centroids(&self, k: usize, inputs: &Matrix<f64>) -> LearningResult<Matrix<f64>> {
         let mut rng = thread_rng();
 
         let mut init_centroids = Vec::with_capacity(k * inputs.cols());
         let first_cen = rng.gen_range(0usize, inputs.rows());
 
-        unsafe { 
-            init_centroids.extend_from_slice(inputs.get_row_unchecked(first_cen)); 
+        unsafe {
+            init_centroids.extend_from_slice(inputs.get_row_unchecked(first_cen));
         }
 
         for i in 1..k {
@@ -342,7 +352,7 @@ impl Initializer for KPlusPlus {
                 }
 
                 let next_cen = sample_discretely(dist);
-                init_centroids.extend_from_slice(inputs.get_row_unchecked(next_cen)); 
+                init_centroids.extend_from_slice(inputs.get_row_unchecked(next_cen));
             }
         }
 
diff --git a/src/learning/lin_reg.rs b/src/learning/lin_reg.rs
index a6de93d9..8577d2e0 100644
--- a/src/learning/lin_reg.rs
+++ b/src/learning/lin_reg.rs
@@ -20,24 +20,25 @@
 //! let mut lin_mod = LinRegressor::default();
 //!
 //! // Train the model
-//! lin_mod.train(&inputs, &targets);
+//! lin_mod.train(&inputs, &targets).unwrap();
 //!
 //! // Now we'll predict a new point
 //! let new_point = Matrix::new(1,1,vec![10.]);
-//! let output = lin_mod.predict(&new_point);
+//! let output = lin_mod.predict(&new_point).unwrap();
 //!
 //! // Hopefully we classified our new point correctly!
 //! assert!(output[0] > 17f64, "Our regressor isn't very good!");
 //! ```
 
-use learning::SupModel;
-use linalg::Matrix;
-use linalg::Vector;
+use learning::{LearningResult, SupModel};
 use learning::toolkit::cost_fn::CostFunc;
 use learning::toolkit::cost_fn::MeanSqError;
 use learning::optim::grad_desc::GradientDesc;
-use learning::optim::OptimAlgorithm;
-use learning::optim::Optimizable;
+use learning::optim::{OptimAlgorithm, Optimizable};
+use learning::error::Error;
+
+use linalg::Matrix;
+use linalg::Vector;
 
 /// Linear Regression Model.
 ///
@@ -80,9 +81,9 @@ impl SupModel<Matrix<f64>, Vector<f64>> for LinRegressor {
     /// let inputs = Matrix::new(3,1, vec![2.0, 3.0, 4.0]);
     /// let targets = Vector::new(vec![5.0, 6.0, 7.0]);
     ///
-    /// lin_mod.train(&inputs, &targets);
+    /// lin_mod.train(&inputs, &targets).unwrap();
     /// ```
-    fn train(&mut self, inputs: &Matrix<f64>, targets: &Vector<f64>) {
+    fn train(&mut self, inputs: &Matrix<f64>, targets: &Vector<f64>) -> LearningResult<()> {
         let ones = Matrix::<f64>::ones(inputs.rows(), 1);
         let full_inputs = ones.hcat(inputs);
 
@@ -91,18 +92,20 @@ impl SupModel<Matrix<f64>, Vector<f64>> for LinRegressor {
         self.parameters =
             Some(((&xt * full_inputs).inverse().expect("Could not compute (X_T X) inverse.") *
                   &xt) * targets);
+
+        Ok(())
     }
 
     /// Predict output value from input data.
     ///
     /// Model must be trained before prediction can be made.
-    fn predict(&self, inputs: &Matrix<f64>) -> Vector<f64> {
+    fn predict(&self, inputs: &Matrix<f64>) -> LearningResult<Vector<f64>> {
         if let Some(ref v) = self.parameters {
             let ones = Matrix::<f64>::ones(inputs.rows(), 1);
             let full_inputs = ones.hcat(inputs);
-            full_inputs * v
+            Ok(full_inputs * v)
         } else {
-            panic!("Model has not been trained.");
+            Err(Error::new_untrained())
         }
     }
 }
@@ -148,7 +151,7 @@ impl LinRegressor {
     ///
     /// // Now we'll predict a new point
     /// let new_point = Matrix::new(1,1,vec![10.]);
-    /// let _ = lin_mod.predict(&new_point);
+    /// let _ = lin_mod.predict(&new_point).unwrap();
     /// ```
     pub fn train_with_optimization(&mut self, inputs: &Matrix<f64>, targets: &Vector<f64>) {
         let ones = Matrix::<f64>::ones(inputs.rows(), 1);
diff --git a/src/learning/logistic_reg.rs b/src/learning/logistic_reg.rs
index ed85b4c6..8e149fb8 100644
--- a/src/learning/logistic_reg.rs
+++ b/src/learning/logistic_reg.rs
@@ -20,11 +20,11 @@
 //! let mut log_mod = LogisticRegressor::default();
 //!
 //! // Train the model
-//! log_mod.train(&inputs, &targets);
+//! log_mod.train(&inputs, &targets).unwrap();
 //!
 //! // Now we'll predict a new point
 //! let new_point = Matrix::new(1,1,vec![10.]);
-//! let output = log_mod.predict(&new_point);
+//! let output = log_mod.predict(&new_point).unwrap();
 //!
 //! // Hopefully we classified our new point correctly!
 //! assert!(output[0] > 0.5, "Our classifier isn't very good!");
@@ -34,13 +34,15 @@
 //! by using the `new` constructor instead. This allows us to provide
 //! a `GradientDesc` object with custom parameters.
 
-use learning::SupModel;
-use linalg::Matrix;
-use linalg::Vector;
+use learning::{LearningResult, SupModel};
 use learning::toolkit::activ_fn::{ActivationFunc, Sigmoid};
 use learning::toolkit::cost_fn::{CostFunc, CrossEntropyError};
 use learning::optim::grad_desc::GradientDesc;
 use learning::optim::{OptimAlgorithm, Optimizable};
+use learning::error::Error;
+
+use linalg::Matrix;
+use linalg::Vector;
 
 /// Logistic Regression Model.
 ///
@@ -110,9 +112,9 @@ impl<A> SupModel<Matrix<f64>, Vector<f64>> for LogisticRegressor<A>
     /// let inputs = Matrix::new(3,2, vec![1.0, 2.0, 1.0, 3.0, 1.0, 4.0]);
     /// let targets = Vector::new(vec![5.0, 6.0, 7.0]);
     ///
-    /// logistic_mod.train(&inputs, &targets);
+    /// logistic_mod.train(&inputs, &targets).unwrap();
     /// ```
-    fn train(&mut self, inputs: &Matrix<f64>, targets: &Vector<f64>) {
+    fn train(&mut self, inputs: &Matrix<f64>, targets: &Vector<f64>) -> LearningResult<()> {
         let ones = Matrix::<f64>::ones(inputs.rows(), 1);
         let full_inputs = ones.hcat(inputs);
 
@@ -120,18 +122,19 @@ impl<A> SupModel<Matrix<f64>, Vector<f64>> for LogisticRegressor<A>
 
         let optimal_w = self.alg.optimize(&self.base, &initial_params[..], &full_inputs, targets);
         self.base.set_parameters(Vector::new(optimal_w));
+        Ok(())
     }
 
     /// Predict output value from input data.
     ///
     /// Model must be trained before prediction can be made.
-    fn predict(&self, inputs: &Matrix<f64>) -> Vector<f64> {
+    fn predict(&self, inputs: &Matrix<f64>) -> LearningResult<Vector<f64>> {
         if let Some(v) = self.base.parameters() {
             let ones = Matrix::<f64>::ones(inputs.rows(), 1);
             let full_inputs = ones.hcat(inputs);
-            (full_inputs * v).apply(&Sigmoid::func)
+            Ok((full_inputs * v).apply(&Sigmoid::func))
         } else {
-            panic!("Model has not been trained.");
+            Err(Error::new_untrained())
         }
     }
 }
diff --git a/src/learning/naive_bayes.rs b/src/learning/naive_bayes.rs
index 3d578aea..20d3619e 100644
--- a/src/learning/naive_bayes.rs
+++ b/src/learning/naive_bayes.rs
@@ -31,18 +31,20 @@
 //! let mut model = NaiveBayes::<Gaussian>::new();
 //!
 //! // Train the model.
-//! model.train(&inputs, &targets);
+//! model.train(&inputs, &targets).unwrap();
 //!
 //! // Predict the classes on the input data
-//! let outputs = model.predict(&inputs);
+//! let outputs = model.predict(&inputs).unwrap();
 //!
 //! // Will output the target classes - otherwise our classifier is bad!
 //! println!("Final outputs --\n{}", outputs);
 //! ```
 
+use learning::{LearningResult, SupModel};
+use learning::error::{Error, ErrorKind};
+
 use linalg::{Matrix, Axes};
 use rulinalg::utils;
-use learning::SupModel;
 
 use std::f64::consts::PI;
 
@@ -105,14 +107,14 @@ impl<T: Distribution> NaiveBayes<T> {
 /// the input class. e.g. [[1,0,0],[0,0,1]] shows class 1 first, then class 3.
 impl<T: Distribution> SupModel<Matrix<f64>, Matrix<f64>> for NaiveBayes<T> {
     /// Train the model using inputs and targets.
-    fn train(&mut self, inputs: &Matrix<f64>, targets: &Matrix<f64>) {
+    fn train(&mut self, inputs: &Matrix<f64>, targets: &Matrix<f64>) -> LearningResult<()> {
         self.distr = Some(T::from_model_params(targets.cols(), inputs.cols()));
-        self.update_params(inputs, targets);
+        self.update_params(inputs, targets)
     }
 
     /// Predict output from inputs.
-    fn predict(&self, inputs: &Matrix<f64>) -> Matrix<f64> {
-        let log_probs = self.get_log_probs(inputs);
+    fn predict(&self, inputs: &Matrix<f64>) -> LearningResult<Matrix<f64>> {
+        let log_probs = try!(self.get_log_probs(inputs));
         let input_classes = NaiveBayes::<T>::get_classes(log_probs);
 
         if let Some(cluster_count) = self.cluster_count {
@@ -125,26 +127,26 @@ impl<T: Distribution> SupModel<Matrix<f64>, Matrix<f64>> for NaiveBayes<T> {
                 class_data.append(&mut row);
             }
 
-            Matrix::new(inputs.rows(), cluster_count, class_data)
+            Ok(Matrix::new(inputs.rows(), cluster_count, class_data))
         } else {
-            panic!("The model has not been trained.");
+            Err(Error::new(ErrorKind::UntrainedModel, "The model has not been trained."))
         }
     }
 }
 
 impl<T: Distribution> NaiveBayes<T> {
     /// Get the log-probabilities per class for each input.
-    pub fn get_log_probs(&self, inputs: &Matrix<f64>) -> Matrix<f64> {
+    pub fn get_log_probs(&self, inputs: &Matrix<f64>) -> LearningResult<Matrix<f64>> {
 
         if let (&Some(ref distr), &Some(ref prior)) = (&self.distr, &self.class_prior) {
             // Get the joint log likelihood from the distribution
-            distr.joint_log_lik(inputs, prior)
+            Ok(distr.joint_log_lik(inputs, prior))
         } else {
-            panic!("Model has not been trained.");
+            Err(Error::new_untrained())
         }
     }
 
-    fn update_params(&mut self, inputs: &Matrix<f64>, targets: &Matrix<f64>) {
+    fn update_params(&mut self, inputs: &Matrix<f64>, targets: &Matrix<f64>) -> LearningResult<()> {
         let class_count = targets.cols();
         let total_data = inputs.rows();
 
@@ -153,7 +155,7 @@ impl<T: Distribution> NaiveBayes<T> {
 
         for (idx, row) in targets.iter_rows().enumerate() {
             // Find the class of this input
-            let class = NaiveBayes::<T>::find_class(row);
+            let class = try!(NaiveBayes::<T>::find_class(row));
 
             // Note the class of the input
             class_data[class].push(idx);
@@ -180,16 +182,18 @@ impl<T: Distribution> NaiveBayes<T> {
 
         self.class_prior = Some(class_prior);
         self.cluster_count = Some(class_count);
+        Ok(())
     }
 
-    fn find_class(row: &[f64]) -> usize {
+    fn find_class(row: &[f64]) -> LearningResult<usize> {
         // Find the `1` entry in the row
         for (idx, r) in row.into_iter().enumerate() {
             if *r == 1f64 {
-                return idx;
+                return Ok(idx);
             }
         }
-        panic!("No class found for entry in targets");
+        
+        Err(Error::new(ErrorKind::InvalidState, "No class found for entry in targets"))
     }
 
     fn get_classes(log_probs: Matrix<f64>) -> Vec<usize> {
@@ -441,10 +445,9 @@ mod tests {
                                        0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0]);
 
         let mut model = NaiveBayes::<Gaussian>::new();
-        model.train(&inputs, &targets);
-
-        let outputs = model.predict(&inputs);
+        model.train(&inputs, &targets).unwrap();
 
+        let outputs = model.predict(&inputs).unwrap();
         assert_eq!(outputs.into_vec(), targets.into_vec());
     }
 
@@ -457,11 +460,9 @@ mod tests {
         let targets = Matrix::new(4, 2, vec![1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0]);
 
         let mut model = NaiveBayes::<Bernoulli>::new();
-        model.train(&inputs, &targets);
-
-        let outputs = model.predict(&inputs);
-        println!("{}", outputs);
+        model.train(&inputs, &targets).unwrap();
 
+        let outputs = model.predict(&inputs).unwrap();
         assert_eq!(outputs.into_vec(), targets.into_vec());
     }
 
@@ -475,11 +476,9 @@ mod tests {
         let targets = Matrix::new(4, 2, vec![1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0]);
 
         let mut model = NaiveBayes::<Multinomial>::new();
-        model.train(&inputs, &targets);
-
-        let outputs = model.predict(&inputs);
-        println!("{}", outputs);
+        model.train(&inputs, &targets).unwrap();
 
+        let outputs = model.predict(&inputs).unwrap();
         assert_eq!(outputs.into_vec(), targets.into_vec());
     }
 }
diff --git a/src/learning/nnet.rs b/src/learning/nnet.rs
index 06a4a731..5b8a6a53 100644
--- a/src/learning/nnet.rs
+++ b/src/learning/nnet.rs
@@ -26,12 +26,12 @@
 //! let mut model = NeuralNet::new(layers, criterion, StochasticGD::default());
 //!
 //! // Train the model!
-//! model.train(&inputs, &targets);
+//! model.train(&inputs, &targets).unwrap();
 //!
 //! let test_inputs = Matrix::new(2,3, vec![1.5,1.5,1.5,5.1,5.1,5.1]);
 //!
 //! // And predict new output from the test inputs
-//! model.predict(&test_inputs);
+//! let outputs = model.predict(&test_inputs).unwrap();
 //! ```
 //!
 //! The neural networks are specified via a criterion - similar to
@@ -44,7 +44,8 @@
 use linalg::{Matrix, MatrixSlice};
 use linalg::BaseSlice;
 
-use learning::SupModel;
+use learning::{LearningResult, SupModel};
+use learning::error::{Error, ErrorKind};
 use learning::toolkit::activ_fn;
 use learning::toolkit::activ_fn::ActivationFunc;
 use learning::toolkit::cost_fn;
@@ -77,14 +78,15 @@ impl<'a, T, A> SupModel<Matrix<f64>, Matrix<f64>> for NeuralNet<'a, T, A>
           A: OptimAlgorithm<BaseNeuralNet<'a, T>>
 {
     /// Predict neural network output using forward propagation.
-    fn predict(&self, inputs: &Matrix<f64>) -> Matrix<f64> {
+    fn predict(&self, inputs: &Matrix<f64>) -> LearningResult<Matrix<f64>> {
         self.base.forward_prop(inputs)
     }
 
     /// Train the model using gradient optimization and back propagation.
-    fn train(&mut self, inputs: &Matrix<f64>, targets: &Matrix<f64>) {
+    fn train(&mut self, inputs: &Matrix<f64>, targets: &Matrix<f64>) -> LearningResult<()> {
         let optimal_w = self.alg.optimize(&self.base, &self.base.weights, inputs, targets);
         self.base.weights = optimal_w;
+        Ok(())
     }
 }
 
@@ -198,8 +200,7 @@ impl<'a, T: Criterion> BaseNeuralNet<'a, T> {
     fn create_weights(layer_sizes: &[usize]) -> Vec<f64> {
         let mut between = range::Range::new(0f64, 1f64);
         let mut rng = thread_rng();
-        layer_sizes
-            .windows(2)
+        layer_sizes.windows(2)
             .flat_map(|w| {
                 let l_in = w[0] + 1;
                 let l_out = w[1];
@@ -207,7 +208,8 @@ impl<'a, T: Criterion> BaseNeuralNet<'a, T> {
                 (0..l_in * l_out)
                     .map(|_i| (between.sample(&mut rng) * 2f64 * eps_init) - eps_init)
                     .collect::<Vec<_>>()
-            }).collect()
+            })
+            .collect()
     }
 
     /// Gets matrix of weights between specified layer and forward layer for the weights.
@@ -342,22 +344,25 @@ impl<'a, T: Criterion> BaseNeuralNet<'a, T> {
     }
 
     /// Forward propagation of the model weights to get the outputs.
-    fn forward_prop(&self, inputs: &Matrix<f64>) -> Matrix<f64> {
-        assert_eq!(inputs.cols(), self.layer_sizes[0]);
+    fn forward_prop(&self, inputs: &Matrix<f64>) -> LearningResult<Matrix<f64>> {
+        if inputs.cols() != self.layer_sizes[0] {
+            Err(Error::new(ErrorKind::InvalidData,
+                           "The input data dimensions must match the first layer."))
+        } else {
+            let net_data = Matrix::ones(inputs.rows(), 1).hcat(inputs);
 
-        let net_data = Matrix::ones(inputs.rows(), 1).hcat(inputs);
+            let mut z = net_data * self.get_net_weights(0);
+            let mut a = self.criterion.activate(z.clone());
 
-        let mut z = net_data * self.get_net_weights(0);
-        let mut a = self.criterion.activate(z.clone());
+            for l in 1..self.layer_sizes.len() - 1 {
+                let ones = Matrix::ones(a.rows(), 1);
+                a = ones.hcat(&a);
+                z = a * self.get_net_weights(l);
+                a = self.criterion.activate(z.clone());
+            }
 
-        for l in 1..self.layer_sizes.len() - 1 {
-            let ones = Matrix::ones(a.rows(), 1);
-            a = ones.hcat(&a);
-            z = a * self.get_net_weights(l);
-            a = self.criterion.activate(z.clone());
+            Ok(a)
         }
-
-        a
     }
 }
 
diff --git a/src/learning/svm.rs b/src/learning/svm.rs
index 4e106035..abd02e82 100644
--- a/src/learning/svm.rs
+++ b/src/learning/svm.rs
@@ -22,11 +22,11 @@
 //! let mut svm_mod = SVM::default();
 //!
 //! // Train the model
-//! svm_mod.train(&inputs, &targets);
+//! svm_mod.train(&inputs, &targets).unwrap();
 //!
 //! // Now we'll predict a new point
 //! let new_point = Matrix::new(1,1,vec![10.]);
-//! let output = svm_mod.predict(&new_point);
+//! let output = svm_mod.predict(&new_point).unwrap();
 //!
 //! // Hopefully we classified our new point correctly!
 //! assert!(output[0] == 1f64, "Our classifier isn't very good!");
@@ -37,7 +37,8 @@ use linalg::Matrix;
 use linalg::Vector;
 
 use learning::toolkit::kernel::{Kernel, SquaredExp};
-use learning::SupModel;
+use learning::{LearningResult, SupModel};
+use learning::error::{Error, ErrorKind};
 
 use rand;
 use rand::Rng;
@@ -101,43 +102,46 @@ impl<K: Kernel> SVM<K> {
 
 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 dim1 = m1.rows();
-        let dim2 = m2.rows();
-
-        let mut ker_data = Vec::with_capacity(dim1 * dim2);
+    fn ker_mat(&self, m1: &Matrix<f64>, m2: &Matrix<f64>) -> LearningResult<Matrix<f64>> {
+        if m1.cols() != m2.cols() {
+            Err(Error::new(ErrorKind::InvalidState,
+                           "Inputs to kernel matrices have different column counts."))
+        } else {
+            let dim1 = m1.rows();
+            let dim2 = m2.rows();
 
-        ker_data.extend(
-            m1.iter_rows().flat_map(|row1| m2.iter_rows()
-                                    .map(move |row2| self.ker.kernel(row1, row2))));
+            let mut ker_data = Vec::with_capacity(dim1 * dim2);
+            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)
+            Ok(Matrix::new(dim1, dim2, ker_data))
+        }
     }
 }
 
 /// Train the model using the Pegasos algorithm and
 /// predict the model output from new data.
 impl<K: Kernel> SupModel<Matrix<f64>, Vector<f64>> for SVM<K> {
-    fn predict(&self, inputs: &Matrix<f64>) -> Vector<f64> {
+    fn predict(&self, inputs: &Matrix<f64>) -> LearningResult<Vector<f64>> {
         let ones = Matrix::<f64>::ones(inputs.rows(), 1);
         let full_inputs = ones.hcat(inputs);
 
         if let (&Some(ref alpha), &Some(ref train_inputs), &Some(ref train_targets)) =
                (&self.alpha, &self.train_inputs, &self.train_targets) {
-            let ker_mat = self.ker_mat(&full_inputs, train_inputs);
+            let ker_mat = try!(self.ker_mat(&full_inputs, train_inputs));
             let weight_vec = alpha.elemul(train_targets) / self.lambda;
 
             let plane_dist = ker_mat * weight_vec;
 
-            plane_dist.apply(&|d| d.signum())
+            Ok(plane_dist.apply(&|d| d.signum()))
         } else {
-            panic!("Model has not been trained.");
+            Err(Error::new_untrained())
         }
     }
 
-    fn train(&mut self, inputs: &Matrix<f64>, targets: &Vector<f64>) {
+    fn train(&mut self, inputs: &Matrix<f64>, targets: &Vector<f64>) -> LearningResult<()> {
         let n = inputs.rows();
 
         let mut rng = rand::thread_rng();
@@ -150,9 +154,9 @@ impl<K: Kernel> SupModel<Matrix<f64>, Vector<f64>> for SVM<K> {
         for t in 0..self.optim_iters {
             let i = rng.gen_range(0, n);
             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));
+            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;
@@ -162,5 +166,7 @@ impl<K: Kernel> SupModel<Matrix<f64>, Vector<f64>> for SVM<K> {
         self.alpha = Some(Vector::new(alpha) / (self.optim_iters as f64));
         self.train_inputs = Some(full_inputs);
         self.train_targets = Some(targets.clone());
+
+        Ok(())
     }
 }
diff --git a/src/lib.rs b/src/lib.rs
index 5323c2e2..82f23115 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -83,10 +83,10 @@
 //! // Now we can train and predict from the model.
 //!
 //! // Train the model!
-//! gp.train(&inputs, &targets);
+//! gp.train(&inputs, &targets).unwrap();
 //!
 //! // Predict output from test datae]
-//! let outputs = gp.predict(&test_inputs);
+//! let outputs = gp.predict(&test_inputs).unwrap();
 //! ```
 //!
 //! This code could have been a lot simpler if we had simply adopted
@@ -109,7 +109,7 @@ extern crate rulinalg;
 extern crate num as libnum;
 extern crate rand;
 
-pub mod prelude;
+pub mod prelude; 
 
 /// The linear algebra module
 ///
@@ -141,22 +141,25 @@ pub mod learning {
 
     pub mod error;
 
+    /// A new type which provides clean access to the learning errors
+    pub type LearningResult<T> = Result<T, error::Error>;
+
     /// Trait for supervised model.
     pub trait SupModel<T, U> {
         /// Predict output from inputs.
-        fn predict(&self, inputs: &T) -> U;
+        fn predict(&self, inputs: &T) -> LearningResult<U>;
 
         /// Train the model using inputs and targets.
-        fn train(&mut self, inputs: &T, targets: &U);
+        fn train(&mut self, inputs: &T, targets: &U) -> LearningResult<()>;
     }
 
     /// Trait for unsupervised model.
     pub trait UnSupModel<T, U> {
         /// Predict output from inputs.
-        fn predict(&self, inputs: &T) -> U;
+        fn predict(&self, inputs: &T) -> LearningResult<U>;
 
         /// Train the model using inputs.
-        fn train(&mut self, inputs: &T);
+        fn train(&mut self, inputs: &T) -> LearningResult<()>;
     }
 
     /// Module for optimization in machine learning setting.
diff --git a/tests/learning/dbscan.rs b/tests/learning/dbscan.rs
index 05338cb9..5761e193 100644
--- a/tests/learning/dbscan.rs
+++ b/tests/learning/dbscan.rs
@@ -13,7 +13,7 @@ fn test_basic_clusters() {
                                         -2.2, 3.1]);
 
     let mut model = DBSCAN::new(0.5, 2);
-    model.train(&inputs);
+    model.train(&inputs).unwrap();
 
     let clustering = model.clusters().unwrap();
 
@@ -33,11 +33,11 @@ fn test_basic_prediction() {
 
     let mut model = DBSCAN::new(0.5, 2);
     model.set_predictive(true);
-    model.train(&inputs);
+    model.train(&inputs).unwrap();
 
     let new_points = Matrix::new(2,2, vec![1.0, 2.0, 4.0, 4.0]);
 
-    let classes = model.predict(&new_points);
+    let classes = model.predict(&new_points).unwrap();
     assert!(classes[0] == Some(0));
     assert!(classes[1] == None);
 }
diff --git a/tests/learning/gp.rs b/tests/learning/gp.rs
index 359e23d4..34c92993 100644
--- a/tests/learning/gp.rs
+++ b/tests/learning/gp.rs
@@ -11,9 +11,9 @@ fn test_default_gp() {
 	let inputs = Matrix::new(10,1,vec![0.,1.,2.,3.,4.,5.,6.,7.,8.,9.]);
 	let targets = Vector::new(vec![0.,1.,2.,3.,4.,4.,3.,2.,1.,0.]);
 
-	gp.train(&inputs, &targets);
+	gp.train(&inputs, &targets).unwrap();
 
 	let test_inputs = Matrix::new(5,1,vec![2.3,4.4,5.1,6.2,7.1]);
 
-	let _outputs = gp.predict(&test_inputs);
+	let _outputs = gp.predict(&test_inputs).unwrap();
 }
diff --git a/tests/learning/k_means.rs b/tests/learning/k_means.rs
index f73da019..168f2597 100644
--- a/tests/learning/k_means.rs
+++ b/tests/learning/k_means.rs
@@ -9,9 +9,9 @@ fn test_model_default() {
     let inputs = Matrix::new(3, 2, vec![1.0, 2.0, 1.0, 3.0, 1.0, 4.0]);
     let targets = Matrix::new(3,2, vec![1.0, 2.0, 1.0, 3.0, 1.0, 4.0]);
 
-    model.train(&inputs);
+    model.train(&inputs).unwrap();
 
-    let outputs = model.predict(&targets);
+    let outputs = model.predict(&targets).unwrap();
 
     assert_eq!(outputs.size(), 3);
 }
@@ -23,9 +23,9 @@ fn test_model_iter() {
     let targets = Matrix::new(3,2, vec![1.0, 2.0, 1.0, 3.0, 1.0, 4.0]);
 
     model.set_iters(1000);
-    model.train(&inputs);
+    model.train(&inputs).unwrap();
 
-    let outputs = model.predict(&targets);
+    let outputs = model.predict(&targets).unwrap();
 
     assert_eq!(outputs.size(), 3);
 }
@@ -36,9 +36,9 @@ fn test_model_forgy() {
     let inputs = Matrix::new(3, 2, vec![1.0, 2.0, 1.0, 3.0, 1.0, 4.0]);
     let targets = Matrix::new(3,2, vec![1.0, 2.0, 1.0, 3.0, 1.0, 4.0]);
 
-    model.train(&inputs);
+    model.train(&inputs).unwrap();
 
-    let outputs = model.predict(&targets);
+    let outputs = model.predict(&targets).unwrap();
 
     assert_eq!(outputs.size(), 3);
 }
@@ -49,9 +49,9 @@ fn test_model_ran_partition() {
     let inputs = Matrix::new(3, 2, vec![1.0, 2.0, 1.0, 3.0, 1.0, 4.0]);
     let targets = Matrix::new(3,2, vec![1.0, 2.0, 1.0, 3.0, 1.0, 4.0]);
 
-    model.train(&inputs);
+    model.train(&inputs).unwrap();
 
-    let outputs = model.predict(&targets);
+    let outputs = model.predict(&targets).unwrap();
 
     assert_eq!(outputs.size(), 3);
 }
@@ -62,9 +62,9 @@ fn test_model_kplusplus() {
     let inputs = Matrix::new(3, 2, vec![1.0, 2.0, 1.0, 3.0, 1.0, 4.0]);
     let targets = Matrix::new(3,2, vec![1.0, 2.0, 1.0, 3.0, 1.0, 4.0]);
 
-    model.train(&inputs);
+    model.train(&inputs).unwrap();
 
-    let outputs = model.predict(&targets);
+    let outputs = model.predict(&targets).unwrap();
 
     assert_eq!(outputs.size(), 3);
 }
@@ -75,7 +75,7 @@ fn test_no_train_predict() {
     let model = KMeansClassifier::<KPlusPlus>::new(3);
     let inputs = Matrix::new(3, 2, vec![1.0, 2.0, 1.0, 3.0, 1.0, 4.0]);
 
-    model.predict(&inputs);
+    model.predict(&inputs).unwrap();
 
 }
 
@@ -89,9 +89,9 @@ fn test_two_centroids() {
                                         314.59375, 174.6875,
                                         350.59375, 161.6875]);
 
-    model.train(&inputs);
+    model.train(&inputs).unwrap();
 
-    let classes = model.predict(&inputs);
+    let classes = model.predict(&inputs).unwrap();
     let class_a = classes[0];
 
     let class_b = if class_a == 0 { 1 } else { 0 };
diff --git a/tests/learning/lin_reg.rs b/tests/learning/lin_reg.rs
index 6871a4cd..92f936ad 100644
--- a/tests/learning/lin_reg.rs
+++ b/tests/learning/lin_reg.rs
@@ -21,7 +21,7 @@ fn test_regression() {
     let inputs = Matrix::new(3, 1, vec![2.0, 3.0, 4.0]);
     let targets = Vector::new(vec![5.0, 6.0, 7.0]);
 
-    lin_mod.train(&inputs, &targets);
+    lin_mod.train(&inputs, &targets).unwrap();
 
     let parameters = lin_mod.parameters().unwrap();
 
@@ -46,5 +46,5 @@ fn test_no_train_predict() {
     let lin_mod = LinRegressor::default();
     let inputs = Matrix::new(3, 2, vec![1.0, 2.0, 1.0, 3.0, 1.0, 4.0]);
 
-    let _ = lin_mod.predict(&inputs);
+    let _ = lin_mod.predict(&inputs).unwrap();
 }
\ No newline at end of file