Alec's Web Log

As a naturally follow up to my last post, I have extended my C++ implementation of matlab like array indexing to work for “zipped matrices”. I often store matrices zipped into std::vectors, column-wise. Now I can access them like how I do in matlab. Suppose in matlab I had:

A = [10,11,12;13,14,15];

You could then suck out a re-ordered block of rows and columns of the matrix by calling:

B = A([1],[3,2]);

Revealing:

B =

    12    11

This is especially convenient for chopping up and reordering rows and columns in linear systems.

I can do this same sort of indexing with C++, the caveat being that I have to pass the number of columns as a third argument. For now, I think I would rather do this than have it be a field of the SmartMatrix.

#include <vector>
#include <stdio.h>

template<typename T>
class SmartMatrix: public std::vector<T>{
  public:
    // act like operator[]
    T operator()(
        size_t _row,
        size_t _col,
        size_t number_of_columns){
      return (*this)[_row*number_of_columns+_col];
    }

    // act like matlab operator()
    SmartMatrix<T> operator()(
        std::vector<size_t>& row_positions,
        std::vector<size_t>& column_positions, 
        size_t number_of_columns){
      SmartMatrix<T> sub;
      sub.resize(row_positions.size()*column_positions.size());
      size_t sub_i = 0;
      for(
          std::vector<size_t>::iterator rpit = row_positions.begin();
          rpit != row_positions.end();
          rpit++,sub_i++){
        size_t sub_j = 0;
        for(
          std::vector<size_t>::iterator cpit = column_positions.begin();
          cpit != column_positions.end();
          cpit++,sub_j++){
          sub[sub_i*column_positions.size() + sub_j] = 
            (*this)[(*rpit)*number_of_columns + (*cpit)];
        }
      }
      return sub;
    }
};

int main(){
  SmartMatrix<int> A;
  size_t cols_in_A = 3;
  size_t rows_in_A = 2;
  A.push_back(10);
  A.push_back(11);
  A.push_back(12);
  A.push_back(13);
  A.push_back(14);
  A.push_back(15);
  A.push_back(16);

  printf("A=\n");
  for(int i=0;i<rows_in_A;i++){
    for(int j=0;j<cols_in_A;j++)
      printf("%d ",A[i*cols_in_A+j]);
    printf("\n");
  }
  printf("\nA[%d,%d]=%d\n",0,2,A(0,2,cols_in_A));

  // make some list of indices
  std::vector<size_t> sub_row_indices;
  sub_row_indices.push_back(0);
  std::vector<size_t> sub_column_indices;
  sub_column_indices.push_back(2);
  sub_column_indices.push_back(1);

  // use operator() to extract entries at those indices
  SmartMatrix<int> B = A(sub_row_indices,sub_column_indices,cols_in_A);
  printf("\nB = A([ ");
  for(int i =0;i<sub_row_indices.size();i++)
    printf("%d ",(int)sub_row_indices[i]);
  printf("],[ ");
  for(int i =0;i<sub_column_indices.size();i++)
    printf("%d ",(int)sub_column_indices[i]);
  printf("]);\n\n");

  printf("B=\n");
  for(int i=0;i<sub_row_indices.size();i++){
    for(int j=0;j<sub_column_indices.size();j++)
      printf("%d ",B[i*sub_column_indices.size()+j]);
    printf("\n");
  }
  
  return 0;
}

Which should output:

A=
10 11 12 
13 14 15 

A[0,2]=12

B = A([ 0 ],[ 2 1 ]);

B=
12 11 

I’ve been doing some C++ coding, working with matrices. In a previous project it was convenient to think of our systems in block form. This translates easily into matlab’s indexing. So say you have in matlab:

a = [10,11,12,13,14,15];

Then you can extract elements 1,3,5 (remember matlab is one-indexed) by calling

b = a([1,3,5]);

which reveals

b =

    10    12    14

So in C++, I’ve made my first attempt at useful templating to extend the std::vector class to have this feature. It is no where near implemented the full features of matlab’s operator(), but it seems to do this particular feature correctly.

#include <vector>
#include <stdio.h>

template<typename T>
class SmartVector : public std::vector<T>{
  public:
    // act like operator[]
    T operator()(size_t _Pos){
      return (*this)[_Pos];
    }
    // act like matlab operator()
    SmartVector<T> operator()(std::vector<size_t>& positions){
      SmartVector<T> sub;
      sub.resize(positions.size());
      size_t sub_i = 0;
      for(
          std::vector<size_t>::iterator pit = positions.begin();
          pit != positions.end();
          pit++,sub_i++){
        sub[sub_i] = (*this)[*pit];
      }
      return sub;
    }
};

int main(){
  // make some vector
  SmartVector<int> a;
  a.push_back(10);
  a.push_back(11);
  a.push_back(12);
  a.push_back(13);
  a.push_back(14);
  a.push_back(15);
  for(int i =0;i<a.size();i++)
    printf("a[%d] = %d\n",i,a[i]);

  // make some list of indices
  std::vector<size_t> sub_indices;
  sub_indices.push_back(0);
  sub_indices.push_back(2);
  sub_indices.push_back(4);

  // use operator() to extract entries at those indices
  SmartVector<int> b = a(sub_indices);
  printf("\nb = a([ ");
  for(int i =0;i<sub_indices.size();i++)
    printf("%d ",(int)sub_indices[i]);
  printf("]);\n\n");

  for(int i =0;i<b.size();i++)
    printf("b[%d] = %d\n",i,b[i]);
  
  return 0;
}

Then if you run ./main you should see:

a[0] = 10
a[1] = 11
a[2] = 12
a[3] = 13
a[4] = 14
a[5] = 15

b = a([ 0 2 4 ]);

b[0] = 10
b[1] = 12
b[2] = 14