playing with subsetting in Eigen

sjb-subsettest.cpp

#include <iostream>
#include <iomanip>
#include <fstream>
#include "mkl.h"
#include "math.h"
#include <vector>
#include <cmath>
#include <string>



#include <boost/random/linear_congruential.hpp>
#include <boost/random/uniform_int.hpp>
#include <boost/random/uniform_real.hpp>
#include <boost/random/variate_generator.hpp>


#define EIGEN_USE_MKL_BLAS
#define EIGEN_USE_MKL_LAPACKE

#ifndef PI
#define PI 3.141592653589793
#endif


#include <Eigen/Dense>

using namespace Eigen;
using namespace std;

// This is a typedef for a random number generator.
// Try boost::mt19937 or boost::ecuyer1988 instead of boost::minstd_rand
typedef boost::minstd_rand base_generator_type;
////////////////////////////////////////////////


// THIS IS A TYPEDEF FOR A ROWMAJOR MATRIX
typedef Eigen::Matrix<double,Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> MatrixRMXd;

MatrixRMXd negSubset(MatrixXd xmat) {
  
  vector<double> nn;
  int N = xmat.rows();
  int T = xmat.cols();

  for(int n = 0; n < N; n++) {
    for(int t = 0; t < T; t++) {
      if( xmat(n,t) < 0) {
	nn.push_back(xmat(n,t));
      }
    }
  } // end double for
    
  MatrixRMXd mysubset = MatrixRMXd(nn.size(), 1);
  for(int i = 0; i < nn.size(); i++) {
    mysubset(i,0) = nn[i];
  }
  
  return(mysubset);
}


// TEST Policy function getter
int main() {


  int N, T;
  N = 15;
  T = 8;
  MatrixXd xmat(N, T);

  for(int n = 0; n < N; n++) {
    for(int t = 0; t < T; t++) {
      int tmp = n * t + t;
      if ( (tmp % 2 ) == 0) {
	tmp = -1 * tmp;
      }
      xmat(n,t) = double(tmp);
    }
  }

  cout << xmat << endl;

  double meanneg = (xmat.array() < 0).select(xmat, 0).sum() / (xmat.array() < 0).count();
  cout << "MEAN NEG: " << meanneg << endl;

  cout << "////////////////////////////////////////" << endl;
  MatrixRMXd myneg = negSubset(xmat);
  cout << myneg.mean() << endl;
}