Mercurial > hg > kwantix

changeset 60:0288895c1071

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix code formatting and most opening curly brackets to their own line
author Jordi Gutiérrez Hermoso <jordigh@gmail.com>
date Sat, 29 May 2010 15:34:17 -0500
parents f8bd0fd20cc8
children b3bf4ac981ec
files src/linalg.cpp src/rbf.cpp src/utils.cpp
diffstat 3 files changed, 351 insertions(+), 187 deletions(-) [+]
line wrap: on
line diff
--- a/src/linalg.cpp
+++ b/src/linalg.cpp
@@ -18,26 +18,30 @@
 
 // *************** Matrix allocation stuff ***************************
 
-matrix::matrix(){
+matrix::matrix()
+{
   A = gsl_matrix_calloc(1,1); //Allocate 1x1 matrix zero matrix.    
   LUfactored = false;
   SVDfactored = false;
 }
 
-matrix::matrix(size_t m, size_t n, double fillvalue){
+matrix::matrix(size_t m, size_t n, double fillvalue)
+{
   A = gsl_matrix_alloc(m,n);
   gsl_matrix_set_all(A,fillvalue);
   LUfactored = false;
   SVDfactored = false;
 }
 
-matrix::matrix(gsl_matrix *M){
+matrix::matrix(gsl_matrix *M)
+{
   A = M;
   LUfactored = false;
   SVDfactored = false;
 }
 
-matrix::matrix(const matrix& M){
+matrix::matrix(const matrix& M)
+{
   A = gsl_matrix_alloc(M.A->size1, M.A->size2);
   gsl_matrix_memcpy(A,M.A);
 
@@ -48,14 +52,17 @@
   condition_number = M.condition_number;
 
   SVDfactored = M.SVDfactored;
-  if(SVDfactored){
+  if(SVDfactored)
+  {
     SVD = gsl_vector_alloc(M.SVD->size);
     gsl_vector_memcpy(SVD,M.SVD);
   } 
 }
 
-matrix matrix::operator=(const matrix& M){
-  if(this != &M){
+matrix matrix::operator=(const matrix& M)
+{
+  if(this != &M)
+  {
     gsl_matrix_free(A);
     A = gsl_matrix_alloc(M.A->size1, M.A->size2);
     gsl_matrix_memcpy(A,M.A);
@@ -64,7 +71,8 @@
       delete LUptr;
 
     LUfactored = M.LUfactored;
-    if(LUfactored){
+    if(LUfactored)
+    {
       LUptr = new LUmatrix(*M.LUptr);
     }
       
@@ -72,7 +80,8 @@
       gsl_vector_free(SVD);
 
     SVDfactored = M.SVDfactored;
-    if(SVDfactored){
+    if(SVDfactored)
+    {
       condition_number = M.condition_number;
       SVD = gsl_vector_alloc(M.SVD->size);
       gsl_vector_memcpy(SVD,M.SVD);
@@ -81,7 +90,8 @@
   return *this;
 }
 
-matrix::~matrix(){
+matrix::~matrix()
+{
   if(A != 0) //Has subclass already deleted this matrix?
     gsl_matrix_free(A);
   if(LUfactored)
@@ -98,7 +108,8 @@
   return precsn;
 }
 
-void matrix::set_precision(size_t p){
+void matrix::set_precision(size_t p)
+{
   precsn = p;
 }
 
@@ -153,18 +164,21 @@
 }
 
 // ******************* Matrix arithmetic stuff ******************
-matrix matrix::operator*(const double  a) const{
+matrix matrix::operator*(const double  a) const
+{
   matrix Z = *this;
   gsl_matrix_scale(Z.A,a);
   return Z;
 }
 
-matrix matrix::operator+(const matrix& N) const{
+matrix matrix::operator+(const matrix& N) const
+{
   matrix Z = *this;
   try{
     gsl_matrix_add(Z.A,N.A);
   }
-  catch(inconformantSizes& exc){
+  catch(inconformantSizes& exc)
+  {
     exc.n_A = A->size1;
     exc.m_A = A->size2;
     exc.n_B = N.A -> size1;
@@ -174,12 +188,15 @@
   return Z;
 }
 
-matrix matrix::operator*(const matrix& N)const{
+matrix matrix::operator*(const matrix& N)const
+{
   matrix Z(A->size1,N.A->size2);
-  try{
+  try
+  {
     gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1,  A, N.A, 0, Z.A);
   }
-  catch(inconformantSizes& exc){
+  catch(inconformantSizes& exc)
+  {
     exc.n_A = A->size1;
     exc.m_A = A->size2;
     exc.n_B = N.A -> size1;
@@ -189,12 +206,15 @@
   return Z;
 }
 
-matrix matrix::operator-(const matrix& N) const{
+matrix matrix::operator-(const matrix& N) const
+{
   matrix Z = *this;
-  try{
+  try
+  {
     gsl_matrix_sub(Z.A,N.A);
   }
-  catch(inconformantSizes& exc){
+  catch(inconformantSizes& exc)
+  {
     exc.n_A = A->size1;
     exc.m_A = A->size2;
     exc.n_B = N.A -> size1;
@@ -204,12 +224,15 @@
   return Z;
 }
 
-vector matrix::operator*(const vector& v) const{
+vector matrix::operator*(const vector& v) const
+{
   vector w(A -> size1);
-  try{
+  try
+  {
     gsl_blas_dgemv(CblasNoTrans, 1, A, v.x, 0, w.x);
   }
-  catch(inconformantSizes& exc){
+  catch(inconformantSizes& exc)
+  {
     exc.n_A = A->size1;
     exc.m_A = A->size2;
     exc.n_B = v.x -> size;
@@ -222,8 +245,10 @@
 // ***************** Other arithmetic functions ******************
 
 
-matrix::LUmatrix* matrix::LU() const{
-  if(A -> size1 != A -> size2){
+matrix::LUmatrix* matrix::LU() const
+{
+  if(A -> size1 != A -> size2)
+  {
     matrixNotSquare exc;
     exc.reason = "Cannot LU-factorise non-square matrices.";
     exc.file = __FILE__;
@@ -231,7 +256,8 @@
     throw exc;
   }
 
-  if(!LUfactored){
+  if(!LUfactored)
+  {
     LUptr = new LUmatrix(A);
     gsl_linalg_LU_decomp(LUptr->matrix_ptr(), LUptr->perm_ptr(), 
                          LUptr->sgn_ptr());
@@ -240,8 +266,10 @@
   return LUptr;    
 }
 
-matrix matrix::inv() const{
-  if(A -> size1 != A -> size2){
+matrix matrix::inv() const
+{
+  if(A -> size1 != A -> size2)
+  {
     matrixNotSquare exc;
     exc.reason = "Cannot invert non-square matrices.";
     exc.file = __FILE__;
@@ -249,7 +277,8 @@
     throw exc;
   }
     
-  if(!LUfactored){
+  if(!LUfactored)
+  {
     LUptr =  LU();    
     LUfactored = true;
   }
@@ -260,7 +289,8 @@
   return Z;
 }
 
-matrix matrix::T()  const{
+matrix matrix::T()  const
+{
   matrix Z(A->size2, A->size1);
   for(size_t i = 1; i <= A->size1; i++)
     for(size_t j = 1; j <= A->size2; j++)
@@ -268,8 +298,10 @@
   return Z;
 }
   
-double matrix::tr() const{
-  if(A -> size1 != A -> size2){
+double matrix::tr() const
+{
+  if(A -> size1 != A -> size2)
+  {
     matrixNotSquare exc;
     exc.reason = "Cannot find trace of non-square matrix.";
     exc.file = __FILE__;
@@ -283,8 +315,10 @@
   return result;
 }
 
-double matrix::det() const{
-  if(A -> size1 != A -> size2){
+double matrix::det() const
+{
+  if(A -> size1 != A -> size2)
+  {
     matrixNotSquare exc;
     exc.reason = "Cannot find determinant of non-square matrices.";
     exc.file = __FILE__;
@@ -292,7 +326,8 @@
     throw exc;
   }
     
-  if( !LUfactored ){
+  if( !LUfactored )
+  {
     LUptr = LU();
     LUfactored = true;
   }
@@ -301,8 +336,10 @@
   return out;
 }
 
-vector matrix::inv(const vector& w) const{
-  if(A -> size1 != A -> size2){
+vector matrix::inv(const vector& w) const
+{
+  if(A -> size1 != A -> size2)
+  {
     matrixNotSquare exc;
     exc.reason = "Cannot invert non-square matrices.";
     exc.file = __FILE__;
@@ -311,7 +348,8 @@
   }
 
   vector v(w.size());
-  if( !LUfactored){
+  if( !LUfactored)
+  {
     LUptr = LU();
     LUfactored = true;
   }
@@ -320,7 +358,8 @@
   return v;
 }
  
-double matrix::cond() const{
+double matrix::cond() const
+{
   if(SVDfactored)
     return condition_number;
 
@@ -330,16 +369,19 @@
   return condition_number;
 }
 
-void matrix::SVDfactor() const{
+void matrix::SVDfactor() const
+{
   if(SVDfactored)
     return;
 
   gsl_matrix * M;
-  if(A -> size1 < A -> size2){ //Transpose so that m>=n
+  if(A -> size1 < A -> size2)
+  { //Transpose so that m>=n
     M = gsl_matrix_alloc(A -> size2, A -> size1);
     gsl_matrix_transpose_memcpy(M,A);
   }
-  else{
+  else
+  {
     M = gsl_matrix_alloc(A -> size1, A -> size2);
     gsl_matrix_memcpy(M,A);
   }
@@ -353,20 +395,21 @@
   gsl_matrix_free(V);
   gsl_matrix_free(M);
 }
-}
+
 
 
 
 //***************** LUmatrix stuff **********************************
-namespace kwantix{
-matrix::LUmatrix::LUmatrix(gsl_matrix* M){
+matrix::LUmatrix::LUmatrix(gsl_matrix* M)
+{
   A = gsl_matrix_alloc(M -> size1, M -> size2);
   gsl_matrix_memcpy(A,M);
   p = gsl_permutation_alloc(A->size1);
   signum = 0;
 }
   
-matrix::LUmatrix::LUmatrix(const LUmatrix& M){
+matrix::LUmatrix::LUmatrix(const LUmatrix& M)
+{
   p = gsl_permutation_alloc(M.p->size);
   gsl_permutation_memcpy(p,M.p);
   A = gsl_matrix_alloc(M.A -> size1, M.A -> size2);
@@ -374,8 +417,10 @@
   signum = M.signum;
 }
 
-matrix::LUmatrix matrix::LUmatrix::operator=(const LUmatrix& M){
-  if(this != &M){
+matrix::LUmatrix matrix::LUmatrix::operator=(const LUmatrix& M)
+{
+  if(this != &M)
+  {
     gsl_permutation_free(p);
     p = gsl_permutation_alloc(M.p->size);
     gsl_permutation_memcpy(p,M.p);
@@ -389,67 +434,77 @@
   return *this;
 }
 
-gsl_matrix* matrix::LUmatrix::matrix_ptr(){
+gsl_matrix* matrix::LUmatrix::matrix_ptr()
+{
   return A;
 }
-gsl_permutation* matrix::LUmatrix::perm_ptr(){
+gsl_permutation* matrix::LUmatrix::perm_ptr()
+{
   return p;
 }
 
-int* matrix::LUmatrix::sgn_ptr(){
+int* matrix::LUmatrix::sgn_ptr()
+{
   return &signum;
 }
   
-int matrix::LUmatrix::sgn(){
+int matrix::LUmatrix::sgn()
+{
   return signum;
 }
 
-matrix::LUmatrix::LUmatrix(){
+matrix::LUmatrix::LUmatrix()
+{
   A = 0;
   p = 0;
   signum = 0;
 }
 
-matrix::LUmatrix::~LUmatrix(){
+matrix::LUmatrix::~LUmatrix()
+{
   gsl_permutation_free(p);
   p = 0;
   gsl_matrix_free(A);
   A = 0;
 }
-}
+
 
 
 //Vector stuff
-
-namespace kwantix{
-
 // **************** Vector allocation stuff ********************* 
   
-vector::vector(){
+vector::vector()
+{
   x = gsl_vector_calloc(1); //Allocate zero vector of size one.
 }
 
-vector::vector(const size_t n, const double fillvalue){
+vector::vector(const size_t n, const double fillvalue)
+{
   x = gsl_vector_alloc(n);
   gsl_vector_set_all(x,fillvalue);
 }
 
-vector::vector(gsl_vector *y){
+vector::vector(gsl_vector *y)
+{
   x = y;
 }
 
-vector::vector(const gsl_vector *y){
+vector::vector(const gsl_vector *y)
+{
   x = gsl_vector_alloc(y -> size);
   gsl_vector_memcpy(x,y);
 }
 
-vector::vector(const vector &y){
+vector::vector(const vector &y)
+{
   x = gsl_vector_alloc(y.x->size);
   gsl_vector_memcpy(x,y.x);
 }
  
-vector& vector::operator=(const vector &y){
-  if(this != &y){
+vector& vector::operator=(const vector &y)
+{
+  if(this != &y)
+  {
     gsl_vector_free(x);
     x = gsl_vector_alloc(y.x -> size);
     gsl_vector_memcpy(x,y.x);
@@ -457,7 +512,8 @@
   return *this;
 }
 
-vector::~vector(){
+vector::~vector()
+{
   if(x != 0) //Has subclass vector_view already deleted this vector?
     gsl_vector_free(x);
 }
@@ -469,7 +525,8 @@
   return precsn;
 }
 
-void vector::set_precision(size_t p){
+void vector::set_precision(size_t p)
+{
   precsn = p;
 }
 
@@ -477,7 +534,8 @@
   return x->size;
 }
   
-vector_view vector::operator()(const slice &s) {
+vector_view vector::operator()(const slice &s) 
+{
   vector_view x_sub(x,s);
   return x_sub;
 }
@@ -502,7 +560,8 @@
 }
 
 vector vector::operator+(const vector& w) const{
-  if(x -> size != w.x->size){
+  if(x -> size != w.x->size)
+  {
     inconformantSizes exc;
     exc.reason = "Cannot add vectors of different sizes.";
     exc.file = __FILE__;
@@ -517,7 +576,8 @@
 }
 
 vector vector::operator-(const vector& w) const{
-  if(x -> size != w.x->size){
+  if(x -> size != w.x->size)
+  {
     inconformantSizes exc;
     exc.reason = "Cannot subtract vectors of different sizes.";
     exc.file = __FILE__;
@@ -537,7 +597,8 @@
   try{
     gsl_vector_mul(u.x, w.x);
   }
-  catch(inconformantSizes& exc){
+  catch(inconformantSizes& exc)
+  {
     exc.reason = "Can't multiply vectors of different sizes.";
     exc.file = __FILE__;
     exc.line = __LINE__;
@@ -554,7 +615,8 @@
   try{
     gsl_vector_div(u.x, w.x);
   }
-  catch(inconformantSizes& exc){
+  catch(inconformantSizes& exc)
+  {
     exc.reason = "Can't divide vectors of different sizes.";
     exc.file = __FILE__;
     exc.line = __LINE__;
@@ -581,7 +643,8 @@
   try{
     gsl_blas_ddot(x, w.x, &a);
   }
-  catch(inconformantSizes& exc){
+  catch(inconformantSizes& exc)
+  {
     exc.reason = "Can't dot product vectors of different sizes.";
     exc.file = __FILE__;
     exc.line = __LINE__;
@@ -594,7 +657,8 @@
 
 //Comparison
 bool vector::operator==(const vector& w) const{
-  if(x -> size != w.x -> size){
+  if(x -> size != w.x -> size)
+  {
     badArgument exc;
     exc.reason = "Cannot compare vectors of different sizes.";
     exc.file = __FILE__;
@@ -613,7 +677,8 @@
   if(x -> size < w.x -> size) //Smaller vectors go first in this order.
     return true;
 
-  for(size_t i = 0; i < x -> size; i++){
+  for(size_t i = 0; i < x -> size; i++)
+  {
     double L = gsl_vector_get(x,i);
     double R = gsl_vector_get(w.x,i);
     if(L < R )
@@ -624,15 +689,18 @@
    
   return false; //Then vectors are equal.
 }
-}
+
 
 //************************* Vector view stuff *******************************
-namespace kwantix{
-vector_view::vector_view(){
+vector_view::vector_view()
+{
   x = 0;
 }
-vector_view::vector_view(gsl_vector* y, const slice &s):vector(new gsl_vector){
-  if(s.end() > y->size - 1){
+vector_view::vector_view(gsl_vector* y, const slice &s)
+  : vector(new gsl_vector)
+{
+  if(s.end() > y->size - 1)
+  {
     indexOutOfRange exc;
     exc.file = __FILE__;
     exc.line = __LINE__;
@@ -648,8 +716,10 @@
   x -> owner = 0;
 }
 
-vector_view::vector_view(gsl_matrix* A, const slice& a, const size_t j ):vector(new gsl_vector){
-  if(a.end() > A -> size1 - 1 or j-1 > A -> size2-1){
+vector_view::vector_view(gsl_matrix* A, const slice& a, const size_t j ):vector(new gsl_vector)
+{
+  if(a.end() > A -> size1 - 1 or j-1 > A -> size2-1)
+  {
     indexOutOfRange exc;
     exc.file = __FILE__;
     exc.line = __LINE__;
@@ -668,8 +738,10 @@
   x -> owner = 0;
 }
 
-vector_view::vector_view(gsl_matrix* A, const size_t i, const slice& b ):vector(new gsl_vector){
-  if(b.end() > A -> size2 - 1 or i-1 > A-> size1-1){
+vector_view::vector_view(gsl_matrix* A, const size_t i, const slice& b ):vector(new gsl_vector)
+{
+  if(b.end() > A -> size2 - 1 or i-1 > A-> size1-1)
+  {
     indexOutOfRange exc;
     exc.file = __FILE__;
     exc.line = __LINE__;
@@ -687,8 +759,10 @@
   x -> owner = 0;
 }
 
-vector_view& vector_view::operator=(const vector& w){
-  if(x->size != w.size()){
+vector_view& vector_view::operator=(const vector& w)
+{
+  if(x->size != w.size())
+  {
     badArgument exc;
     exc.reason =  "Cannot assign to vector view: incomformant sizes.";
     exc.file = __FILE__;
@@ -699,8 +773,10 @@
   return *this;
 }
 
-vector_view& vector_view::operator=(const vector_view& w){
-  if(x->size != w.size()){
+vector_view& vector_view::operator=(const vector_view& w)
+{
+  if(x->size != w.size())
+  {
     badArgument exc;
     exc.reason =  "Cannot assign to vector view: incomformant sizes.";
     exc.file = __FILE__;
@@ -712,26 +788,29 @@
     
 }
 
-vector_view::~vector_view(){
+vector_view::~vector_view()
+{
   delete x;
   x = 0;
 }
-}
+
 
 
 //Slice stuff
-namespace kwantix{
-
-slice::slice(size_t a, size_t b, size_t k) {
+slice::slice(size_t a, size_t b, size_t k) 
+{
   set(a,b,k);
 } 
 
-slice slice::operator()(size_t a, size_t b, size_t k) {
+slice slice::operator()(size_t a, size_t b, size_t k) 
+{
   return set(a,b,k);
 }
 
-slice slice::set(size_t a, size_t b, size_t k) {
-  if(b < a){
+slice slice::set(size_t a, size_t b, size_t k) 
+{
+  if(b < a)
+  {
     badArgument exc;
     exc.reason = 
       "Invalid arguments to slice::set. "
@@ -739,7 +818,8 @@
     exc.line = __LINE__; exc.file = __FILE__;
     throw exc;
   }
-  else if(k==0){
+  else if(k==0)
+  {
     badArgument exc;
     exc.reason = 
       "Invalid arguments to slice::set. "
@@ -747,7 +827,8 @@
     exc.line = __LINE__; exc.file = __FILE__;
     throw exc;
   }
-  else if( k > b-a){
+  else if( k > b-a)
+  {
     badArgument exc;
     exc.reason = "Invalid arguments to slice::set. " 
       "Step size cannot be greater than range.";
@@ -760,29 +841,29 @@
   return *this;
 }
 
-}
 
 //Non-member functions
-namespace kwantix{
-  
 // *************      I/O       ************************
-  
-std::ostream& operator<<(std::ostream& os, const vector &v){
+std::ostream& operator<<(std::ostream& os, const vector &v)
+{
   os.setf(std::ios::scientific);
   os << std::setprecision(int(v.precision()) );
-  for(size_t i = 1; i <= v.size(); i++){
+  for(size_t i = 1; i <= v.size(); i++)
+  {
     os << " " << std::setw(int(v.precision()+6) ) << v(i) ;
   }
   return os;
 }
 
-vector operator>>(std::istream& is, vector& v){
+vector operator>>(std::istream& is, vector& v)
+{
   using namespace std;
   string s;
   list<double> data;
   bool colvector = true;
   bool shouldbedone = false;
-  while(getline(is, s)){
+  while(getline(is, s))
+  {
     s = kwantix::trim(s);
     if(s[0] == '#' or s.size() == 0) //Blank line or comment character
       continue;
@@ -792,8 +873,10 @@
 
     double x;
     size_t i = 0;
-    while(ss >> x){
-      if( (i > 0 and colvector == false) or (shouldbedone == true)){
+    while(ss >> x)
+    {
+      if( (i > 0 and colvector == false) or (shouldbedone == true))
+      {
         badArgument exc;
         exc.reason = "Cannot read vector: bad format in input";
         exc.file = __FILE__;
@@ -803,13 +886,15 @@
       data.push_back(x);
       i++;
     }
-    if(i > 1){
+    if(i > 1)
+    {
       colvector = false; //So it must be a row vector instead
       shouldbedone = true;
     }
   }
 
-  if(data.size() == 0){
+  if(data.size() == 0)
+  {
     endOfFile exc;
     exc.reason = "Cannot read empty vector from input.";
     exc.file = __FILE__;
@@ -817,10 +902,11 @@
     throw exc;
   }
 
-  vector w(data.size());
+.  vector w(data.size());
   typedef list<double>::iterator LI;
   size_t k = 1;
-  for(LI i = data.begin(); i != data.end(); i++){
+  for(LI i = data.begin(); i != data.end(); i++)
+  {
     w(k) = *i;
     k++;
   }
@@ -828,10 +914,12 @@
   return v;
 }
 
-std::ostream& operator<<(std::ostream& os, const matrix& A){
+std::ostream& operator<<(std::ostream& os, const matrix& A)
+{
   os.setf(std::ios::scientific);
   os << std::setprecision(int(A.precision()) );
-  for(size_t i = 1; i <= A.rows(); i++){
+  for(size_t i = 1; i <= A.rows(); i++)
+  {
     for(size_t j = 1; j <= A.cols(); j++)
       os << " " << std::setw(int(A.precision()+6) ) << A(i,j) << " ";
     os << std::endl;
@@ -839,7 +927,8 @@
   return os;
 }
 
-matrix operator>>(std::istream& is, matrix& A){
+matrix operator>>(std::istream& is, matrix& A)
+{
   using namespace std;
   string line, token;
   bool rowset = false;
@@ -848,7 +937,8 @@
   size_t rowsize = 0;
   size_t rows = 0;
   size_t cols = 0;
-  while(getline(is, line)){
+  while(getline(is, line))
+  {
     line = kwantix::trim(line);
     //Blank row or comment character.
     if(line[0] == '#' or line.length() == 0) 
@@ -857,8 +947,10 @@
     stringstream ss_line;      
     cols = 0;
     ss_line << line;
-    while(ss_line >> token){
-      if(token[0] == '#'){
+    while(ss_line >> token)
+    {
+      if(token[0] == '#')
+      {
         break; //Rest of line is comment.
       }
 
@@ -866,22 +958,26 @@
       //obey IEEE arithmetic (IEC 559). We could check for those,
       //but do we really want to compile Octave on C++
       //implementations that don't follow IEEE arithmetic?
-      else if(token == "NaN"){ 
+      else if(token == "NaN")
+      { 
         double x = std::numeric_limits<double>::quiet_NaN();
         data.push_back(x);
         cols++;
       }
-      else if(token == "Inf"){
+      else if(token == "Inf")
+      {
         double x = std::numeric_limits<double>::infinity();
         data.push_back(x);
         cols++;
       }
-      else if(token == "-Inf"){
+      else if(token == "-Inf")
+      {
         double x = -std::numeric_limits<double>::infinity();
         data.push_back(x);
         cols++;
       }
-      else if(token == ","){
+      else if(token == ",")
+      {
         ss_line >> token;
       }
 
@@ -901,12 +997,14 @@
        
     //First row gives the number of columns, and all successive rows
     //must have the same number of elements. 
-    if(!rowset){
+    if(!rowset)
+    {
       rowset = true;
       rowsize = cols;
     }
        
-    if (cols != rowsize){
+    if (cols != rowsize)
+    {
       badArgument exc;
       exc.reason = "Cannot read matrix: bad format in input";
       exc.file = __FILE__;
@@ -915,7 +1013,8 @@
     }
     rows++;
   }
-  if(rows == 0){
+  if(rows == 0)
+  {
     endOfFile exc;
     exc.reason = "Cannot read empty matrix from input.";
     exc.file = __FILE__;
@@ -927,8 +1026,10 @@
   typedef list<double>::iterator LI;
     
   LI k = data.begin();
-  for(size_t i = 1; i <= rows; i++){
-    for(size_t j = 1; j <= cols; j++){
+  for(size_t i = 1; i <= rows; i++)
+  {
+    for(size_t j = 1; j <= cols; j++)
+    {
       M(i,j) = *k;
       k++;
     }
@@ -944,6 +1045,7 @@
 {
   return v*a;
 }
+
 double norm(const vector& v)
 {
   return v.norm();
@@ -953,18 +1055,22 @@
 {
   return M*a;
 }
+
 matrix inv(const matrix& A)
 {
   return A.inv();
 }
+
 matrix T(const matrix& A)
 {
   return A.T();
 }
+
 double tr(const matrix& A)
 {
   return A.tr();
 }
+
 double det(matrix& A)
 {
   return A.det();
@@ -975,4 +1081,4 @@
   return A.cond();
 }
 
-}
+}//namespace kwantix
--- a/src/rbf.cpp
+++ b/src/rbf.cpp
@@ -17,37 +17,44 @@
 size_t piecewise_polynomial::n = 0;
 size_t thin_plate_spline::n = 0;
 
-radial_basis_function::radial_basis_function(){
-  if(dimension != 0){
+radial_basis_function::radial_basis_function()
+{
+  if(dimension != 0)
+  {
     point zero(dimension); //Centred by default at the origin.
     centre = zero;
   }
 }
 
-radial_basis_function::radial_basis_function(const point& c){
+radial_basis_function::radial_basis_function(const point& c)
+{
   if(c.size() != dimension)
     bad_dimension(__FILE__, __LINE__, c.size());
   else
     centre = c;
 }
   
-radial_basis_function::~radial_basis_function(){
+radial_basis_function::~radial_basis_function()
+{
   //Nothing to destroy!
 }
   
-void radial_basis_function::set_centre(const point& c){
+void radial_basis_function::set_centre(const point& c)
+{
   if(c.size () != dimension)
     bad_dimension(__FILE__, __LINE__, c.size());
   else
     centre = c;
 }
 
-void radial_basis_function::set_dimension(size_t dim){
+void radial_basis_function::set_dimension(size_t dim)
+{
   dimension = dim;
 }
 
 void radial_basis_function::bad_dimension(string file, 
-                                          int line, size_t dim) const{
+                                          int line, size_t dim) const
+{
   kwantix::badDimension exc;
   if(dimension == 0)
     exc.reason = 
@@ -66,20 +73,24 @@
   throw exc;
 }
 
-double radial_basis_function::at(const point& x) const{
+double radial_basis_function::at(const point& x) const
+{
   return operator()(x);
 }
   
-double radial_basis_function::operator()(const point& x) const{
+double radial_basis_function::operator()(const point& x) const
+{
   if(x.size() != dimension)
     bad_dimension(__FILE__, __LINE__, x.size());
   return operator()( norm(x - centre) );
 }
   
-double radial_basis_function::d(const point& x, size_t k) const{
+double radial_basis_function::d(const point& x, size_t k) const
+{
   if(x.size() != dimension)
     bad_dimension(__FILE__, __LINE__, x.size());
-  else if(k < 1 or k > dimension){
+  else if(k < 1 or k > dimension)
+  {
     kwantix::badArgument exc;
     exc.reason = "Cannot differentiate wrt given index: out of bounds.";
     exc.line = __LINE__;
@@ -96,10 +107,12 @@
 }
 
 double radial_basis_function::d2(const point& x, 
-                                 size_t k1, size_t k2) const{
+                                 size_t k1, size_t k2) const
+{
   if(x.size() != dimension)
     bad_dimension(__FILE__, __LINE__, x.size());
-  else if(k1 < 1 or k1 > dimension or k2 < 1 or k2 > dimension){
+  else if(k1 < 1 or k1 > dimension or k2 < 1 or k2 > dimension)
+  {
     kwantix::badArgument exc;
     exc.reason = "Cannot differentiate wrt given indices: out of bounds.";
     exc.line = __LINE__;
@@ -110,7 +123,8 @@
     return d2(0);
   else if(x == centre and k1 != k2)
     return 0;
-  else if(k1 != k2){
+  else if(k1 != k2)
+  {
     double r = norm(x-centre);
     double r2 = r*r;
     double top = (x(k1) - centre(k1)) * (x(k2)-centre(k2));
@@ -147,24 +161,30 @@
 
 namespace kwantix{
 
-piecewise_smooth_rbf::piecewise_smooth_rbf(){
+piecewise_smooth_rbf::piecewise_smooth_rbf()
+{
   //Nothing to create!
 }
   
-piecewise_smooth_rbf::~piecewise_smooth_rbf(){
+piecewise_smooth_rbf::~piecewise_smooth_rbf()
+{
   //Nothing to destroy!
 }
 
-c_infty_rbf::c_infty_rbf(){
+c_infty_rbf::c_infty_rbf()
+{
   //Nothing to create!
 }
 
-c_infty_rbf::~c_infty_rbf(){
+c_infty_rbf::~c_infty_rbf()
+{
   //Nothing to destroy!
 }
 
-void piecewise_polynomial::set_n(size_t n_new){
-  if(n_new % 2 != 1){ 
+void piecewise_polynomial::set_n(size_t n_new)
+{
+  if(n_new % 2 != 1)
+  { 
     badArgument exc;
     exc.reason = "Cannot assign an even n to a piecewise polynomial RBF.";
     exc.line = __LINE__;
@@ -174,8 +194,10 @@
   piecewise_polynomial::n = n_new;    
 }
 
-void thin_plate_spline::set_n(size_t n_new){
-  if(n_new % 2 != 0){ 
+void thin_plate_spline::set_n(size_t n_new)
+{
+  if(n_new % 2 != 0)
+  { 
     badArgument exc;
     exc.reason = "Cannot assign an odd n to a thin-plate spline RBF.";
     exc.line = __LINE__;
@@ -185,8 +207,10 @@
   thin_plate_spline::n = n_new;    
 }  
 
-void c_infty_rbf::set_epsilon(double e){
-  if(e <= 0){
+void c_infty_rbf::set_epsilon(double e)
+{
+  if(e <= 0)
+  {
     badArgument exc;
     exc.reason = "C-infinity RBFs must have a positive epsilon.";
     exc.line = __LINE__;
@@ -201,8 +225,10 @@
 
 //Piecewise polynomial
 namespace kwantix{
-double piecewise_polynomial::operator()(double r) const{
-  if(n == 0){
+double piecewise_polynomial::operator()(double r) const
+{
+  if(n == 0)
+  {
     badArgument exc;
     exc.reason = 
       "Parameter n not set for piecewise_polynomial. \n"
@@ -226,8 +252,10 @@
   return pow(r,double(n));
 }
   
-double piecewise_polynomial::d(double r) const {
-  if(n == 0){
+double piecewise_polynomial::d(double r) const
+{
+  if(n == 0)
+  {
     badArgument exc;
     exc.reason = 
       "Parameter n not set for piecewise_polynomial. \n"
@@ -250,8 +278,10 @@
   return double(n)*pow(r,double(n-1));
 }
 
-double piecewise_polynomial::d2(double r) const {
-  if(n == 0){
+double piecewise_polynomial::d2(double r) const
+{
+  if(n == 0)
+  {
     badArgument exc;
     exc.reason = 
       "Parameter n not set for piecewise_polynomial. \n"
@@ -279,8 +309,10 @@
 
 //Thin-plate spline
 namespace kwantix{
-double thin_plate_spline::operator()(double r) const {
-  if(n == 0){
+double thin_plate_spline::operator()(double r) const
+{
+  if(n == 0)
+  {
     badArgument exc;
     exc.reason = 
       "Parameter n not set for thin_plate_spline. \n"
@@ -305,8 +337,10 @@
 }
   
  
-double thin_plate_spline::d(double r) const {
-  if(n == 0){
+double thin_plate_spline::d(double r) const
+{
+  if(n == 0)
+  {
     badArgument exc;
     exc.reason = 
       "Parameter n not set for thin_plate_spline. \n"
@@ -330,8 +364,10 @@
   return pow(r,double(n-1))*(double(n)*log(r) + 1);
 }
 
-double thin_plate_spline::d2(double r) const {
-  if(n == 0){
+double thin_plate_spline::d2(double r) const
+{
+  if(n == 0)
+  {
     badArgument exc;
     exc.reason = 
       "Parameter n not set for thin_plate_spline. \n"
@@ -340,7 +376,8 @@
     exc.file = __FILE__;
     throw exc;
   }
-  if (n == 2 and r == 0){
+  if (n == 2 and r == 0)
+  {
     badDomain exc;
     exc.reason = 
       "For n = 2, thin-plate splines do no have a derivative at zero.";
@@ -366,15 +403,18 @@
 
 //Multiquadric
 namespace kwantix{
-double multiquadric::operator()(double r) const {
+double multiquadric::operator()(double r) const
+{
   return sqrt(1 + gsl_pow_2(eps*r));
 }
   
-double multiquadric::d(double r)const{
+double multiquadric::d(double r)const
+{
   return eps*eps*fabs(r)/sqrt(1 + gsl_pow_2(eps*r));
 }
 
-double multiquadric::d2(double r)const{
+double multiquadric::d2(double r)const
+{
   double temp = sqrt(1 + gsl_pow_2(eps*r));
   return 
     eps*eps/temp - gsl_pow_4(eps)*r*r/gsl_pow_3(temp);
@@ -383,15 +423,18 @@
 
 //Inverse multiquadric
 namespace kwantix{
-double inverse_multiquadric::operator()(double r) const {
+double inverse_multiquadric::operator()(double r) const
+{
   return 1/sqrt(1 + gsl_pow_2(eps*r));
 }
   
-double inverse_multiquadric::d(double r) const {
+double inverse_multiquadric::d(double r) const
+{
   return -eps*eps*fabs(r)/gsl_pow_3(sqrt(1 + gsl_pow_2(eps*r)));
 }
 
-double inverse_multiquadric::d2(double r) const {
+double inverse_multiquadric::d2(double r) const
+{
   double temp = sqrt(1 + gsl_pow_2(eps*r));
   return
     - eps*eps/gsl_pow_3(temp) + 3*gsl_pow_4(eps)*r*r/gsl_pow_5(temp);
@@ -400,15 +443,18 @@
 
 //Inverse quadratic
 namespace kwantix{
-double inverse_quadratic::operator()(double r) const {
+double inverse_quadratic::operator()(double r) const
+{
   return 1/(1 + gsl_pow_2(eps*r));
 }
   
-double inverse_quadratic::d(double r)const{
+double inverse_quadratic::d(double r)const
+{
   return -2*eps*eps*fabs(r)/gsl_pow_2(1 + gsl_pow_2(eps*r));
 }
 
-double inverse_quadratic::d2(double r)const{
+double inverse_quadratic::d2(double r)const
+{
   double temp = 1 + gsl_pow_2(eps*r);
   double temp2 = temp*temp, temp3 = temp2*temp;
   return 
@@ -419,15 +465,18 @@
 
 //Gaussian
 namespace kwantix{
-double gaussian::operator()(double r) const {
+double gaussian::operator()(double r) const
+{
   return exp(-gsl_pow_2(eps*r));
 }
   
-double gaussian::d(double r) const {
+double gaussian::d(double r) const
+{
   return -2*eps*eps*fabs(r) * exp(-gsl_pow_2(eps*r));
 }
 
-double gaussian::d2(double r) const {
+double gaussian::d2(double r) const
+{
   double temp = exp(-gsl_pow_2(eps*r));
   return eps*eps*(4*gsl_pow_2(eps*r) - 2) * temp;
       
--- a/src/utils.cpp
+++ b/src/utils.cpp
@@ -8,7 +8,8 @@
 #include "include/error.hpp"
 
 namespace kwantix{
-std::string trim(const std::string& s){
+std::string trim(const std::string& s)
+{
   if(s.length() == 0)
     return s;
   std::size_t beg = s.find_first_not_of(" \a\b\f\n\r\t\v");
@@ -19,24 +20,28 @@
 }
 
 template<typename K, typename V> 
-bool contains(const std::map<K,V>& m, K thing){
+bool contains(const std::map<K,V>& m, K thing)
+{
   return m.find(thing) != m.end();
 }
 
 template<typename E>
-bool contains(const std::set<E>& s, E thing){
+bool contains(const std::set<E>& s, E thing)
+{
   return s.find(thing) != s.end();
 }
   
 template<typename E>
-bool includes(const std::set<E>& s1, const std::set<E>& s2){
+bool includes(const std::set<E>& s1, const std::set<E>& s2)
+{
   return std::includes(s2.begin(), s2.end(), s1.begin(), s1.end());
 }
   
 using kwantix::vector;
 using kwantix::matrix;
 
-matrix read_matrix(std::string filename){
+matrix read_matrix(std::string filename)
+{
   std::ifstream ifs(filename.c_str());
   if(!ifs){
     kwantix::badArgument exc;
@@ -51,7 +56,8 @@
   return v;
 }
 
-vector read_vector(std::string filename){
+vector read_vector(std::string filename)
+{
   std::ifstream ifs(filename.c_str());
   if(!ifs){
     kwantix::badArgument exc;
@@ -66,7 +72,8 @@
   return v;
 }
 
-std::map<kwantix::point, double> read_pd_map(std::string filename){
+std::map<kwantix::point, double> read_pd_map(std::string filename)
+{
   std::ifstream ifs(filename.c_str());
   if(!ifs){
     kwantix::badArgument exc;
@@ -97,17 +104,19 @@
   return result;
 }
 
-void show_exception(kwantix::error exc){
+void show_exception(kwantix::error exc)
+{
   using namespace std;
 
   cerr << "Caught an exception!" << endl;
   cerr << exc.file << ":" << exc.line << " error: "<< exc.reason << endl;
 }
 
-}
+}//namespace kwantix
 
 #include <boost/shared_ptr.hpp>
 #include "include/ddm.hpp"
+
 //Instantiations
 namespace kwantix{
 using boost::shared_ptr;
@@ -124,4 +133,4 @@
                        <const kwantix::overlapping_domain> >&, 
                        shared_ptr<const kwantix::overlapping_domain> E); 
 
-}
+}//namespace kwantix