text: #include "brook.h"
  
  using Brook::Graph;
  
  Graph::Graph()
  : brook_elemwisesub (1, 9)
  {
      brook_elemwisesub(0, 0) = 0.951322;
      brook_elemwisesub(0, 1) = 1.73086;
      brook_elemwisesub(0, 2) = 2.7857;
      brook_elemwisesub(0, 3) = 3.43227;
      brook_elemwisesub(0, 4) = 4.55248;
      brook_elemwisesub(0, 5) = 5.22439;
      brook_elemwisesub(0, 6) = 6.21409;
      brook_elemwisesub(0, 7) = 7.18887;
      brook_elemwisesub(0, 8) = 7.9813;
      
  }
  
  Eigen::MatrixXd Graph::convertToEigenMatrix(std::vector<std::vector<double> > data) {
      Eigen::MatrixXd mat(data.size(), data[0].size());
      for (uint i = 0; i < data.size(); i++){
          mat.row(i) = Eigen::VectorXd::Map(&data[i][0], data[i].size());
      }
  
      return mat;
  }
  
  Eigen::MatrixXd sub(Eigen::MatrixXd a, Eigen::MatrixXd b) {
  
      if ((a.rows() == b.rows()) and (a.cols() == b.cols())) // No need for broadcasting
          return a - b;
  
      // Full broadcast from a
      else if ((a.rows() == 1) and (a.cols() == 1) )
          return b.array() - a(0, 0);
  
      // Full broadcast from b
      else if ((b.rows() == 1) and (b.cols() == 1) )
          return a.array() - b(0, 0);
  
      // Row-wise broadcast from a
      else if ((a.rows() == 1) && (a.cols() == b.cols())) {
          Eigen::Map<Eigen::RowVectorXd> *a_vector = new Eigen::Map<Eigen::RowVectorXd>(a.data(), a.cols());
          return b.rowwise() - *a_vector;
  
      // Row-wise broadcast from b
      } else if ((b.rows() == 1) && (a.cols() == b.cols())) {
          Eigen::Map<Eigen::RowVectorXd>  *b_vector = new Eigen::Map<Eigen::RowVectorXd>(b.data(), b.cols());
          return a.rowwise() - *b_vector;
  
      // Column-wise broadcast from a
      } else if ((a.cols() == 1) && (a.rows() == b.rows())) {
           Eigen::Map<Eigen::VectorXd>  *a_vector = new Eigen::Map<Eigen::VectorXd>(a.data(), a.rows());
           return b.colwise() - *a_vector;
  
      // Column-wise broadcast from b
      } else if ((b.cols() == 1) && (a.rows() == b.rows())) {
           Eigen::Map<Eigen::VectorXd>  *b_vector = new Eigen::Map<Eigen::VectorXd>(b.data(), b.rows());
           return a.colwise() - *b_vector;
  
      } else { // Must never come here normally
          throw "Abnormal matrix dimensions for addition. Exiting...";
      }
  }
  Eigen::MatrixXd Graph::infer(Eigen::MatrixXd brook_inputx){
  
      Eigen::MatrixXd brook_outputx = sub(brook_inputx, brook_elemwisesub);
      return brook_outputx;
      
  }
  