Added some general sorting functions

miscmaths.cc

@@ -294,10 +294,12 @@ namespace MISCMATHS {
   
   int write_ascii_matrix(const Matrix& mat, ofstream& fs, int precision)
   {
+    fs.setf(ios::floatfield);  // use fixed or scientific notation as appropriate
     if (precision>0)  { 
-      fs.setf(ios::scientific | ios::showpos); 
       fs.precision(precision); 
-    } 
+    } else {
+      fs.precision(10);    // default precision
+    }
 #ifdef PPC64	
     int n=0;
 #endif
@@ -618,6 +620,53 @@ namespace MISCMATHS {
     }
 
 
+  vector<int> get_sortindex(const Matrix& vals, const string& mode, int col)
+  {
+    // mode is either "new2old" or "old2new"
+    // return the mapping of old and new indices in the *ascending* sort of vals (from column=col)
+    int length=vals.Nrows();
+    vector<pair<double, int> > sortlist(length);
+    for (int n=0; n<length; n++) {
+      sortlist[n] = pair<double, int>((double) vals(n+1,col),n+1);
+    }
+    sort(sortlist.begin(),sortlist.end());  // O(N.log(N))
+    vector<int> idx(length);
+    for (int n=0; n<length; n++) {
+      if (mode=="old2new") {
+	//  here idx[n] is the where in the ordered list the old n'th row is mapped to (i.e. idx[n] = rank)
+	idx[sortlist[n].second-1] = n+1;
+      } else if (mode=="new2old") {
+	// here idx[n] is the the old row number of the n'th ordered item (i.e. idx[n] is old row number with rank = n)
+	idx[n] = sortlist[n].second;  
+      } else {
+	cerr << "ERROR:: unknown mode in get_sortidx = " << mode << endl; 
+      }
+    }
+    return idx;
+  }
+
+  Matrix apply_sortindex(const Matrix& vals, vector<int> sidx, const string& mode)
+  {
+    // mode is either "new2old" or "old2new"
+    // apply the index mapping from get_sortindex to the whole matrix (swapping rows)
+    Matrix res(vals);
+    res=0.0;
+    int ncols=vals.Ncols();
+    for (unsigned int n=0; n<sidx.size(); n++) {
+      int row = sidx[n];
+      if (mode=="old2new") {
+	res.SubMatrix(row,row,1,ncols)=vals.SubMatrix(n+1,n+1,1,ncols);
+      } else if (mode=="new2old") {
+	res.SubMatrix(n+1,n+1,1,ncols)=vals.SubMatrix(row,row,1,ncols);
+      } else {
+	cerr << "ERROR:: unknown mode in apply_sortidx = " << mode << endl; 
+      }
+    }
+    return res;
+  }
+
+
+
   //------------------------------------------------------------------------//
 
   // Handy MATLAB-like functions
@@ -1143,9 +1192,8 @@ float interp1(const ColumnVector& x, const ColumnVector& y, float xi)
     if(xi <= x.Minimum()) 
       ans = y(1); 
     else{
-      int ind=1;
-      while(xi >= x(ind))
-	ind++;      
+      int ind=2;
+      while(xi >= x(ind)) { ind++; }
       float xa = x(ind-1), xb = x(ind), ya = y(ind-1), yb = y(ind);
       ans = ya + (xi - xa)/(xb - xa) * (yb - ya);
     }

miscmaths.h

@@ -133,6 +133,10 @@ namespace MISCMATHS {
   int rank(const Matrix& X);
   ReturnMatrix sqrtaff(const Matrix& mat);
 
+  // in the following mode = "new2old" or "old2new" (see .cc for more info)
+  vector<int> get_sortindex(const Matrix& vals, const string& mode, int col=1);
+  Matrix apply_sortindex(const Matrix& vals, vector<int> sidx, const string& mode);
+
   void reshape(Matrix& r, const Matrix& m, int nrows, int ncols);
   ReturnMatrix reshape(const Matrix& m, int nrows, int ncols);
   int addrow(Matrix& m, int ncols);