1

My program is short and simple, designed to test my PCA regression class.

#include <iostream>
#include <vector>
#include "data.generator.h"
#include "reg.class.h"

int main(int argc, char* argv[]) {
  std::cout << "Flag A" << std::endl;
  basicGenerator bg;
  std::cout << "Sample size: " << bg.get_sampleSize() << std::endl;
  bg.makeData();
  std::cout << "Flag B" << std::endl;
  std::vector<std::vector<double> > x;
  std::vector<double> y;
  bg.getDataForRegression(x,y);
  unsigned int imax = y.size();
  for (unsigned int i = 0 ; i < imax ; i++) {
    std::cout << i << "\t" << y[i] << "\t" << x[i][0] << "\t" << x[i][1] << std::endl;
  }
  std::cout << "Flag C" << std::endl;
  regressionPCA rpca(x,y);
      return 0;
}

Here is the tail end of the output it produces:

    89      0.416356        0.565407        -0.114605
    90      -0.258883       -0.585597       -0.183372
    91      0.833355        0.914344        0.81171
    92      -0.0338814      -0.00264442     -0.118973
    93      1.13764 0.41599 1.33175
    94      -1.86323        -1.90867        -1.35118
    95      0.907604        1.14917 0.621669
    96      2.1166  1.06194 1.1703
    97      0.159543        0.14446 -0.665135
    98      -0.508617       -0.370597       -0.703225
    99      2.69086 2.75267 1.40633
    Flag C

    Ted@Ted-acer-i7w7 ~/New.System/tests
    $

As you can see, the last output is from the flag I set immediately before creating my instance of regressionPCA on the stack.

Here is the declaration of that class.

#ifndef REG_CLASS_H
#define REG_CLASS_H

#include <iosfwd>
#include <vector>
//#include <boost/smart_ptr/shared_array.hpp>
#include <gsl/gsl_linalg.h>

class regressionPCA {
 private:
  unsigned int nrows, ncols;
  //  boost::shared_array<double> B, Y, U, V, S;
  regressionPCA(void) {}; // makes default construction impossible
 public:
  regressionPCA(const std::vector<std::vector<double> >&, const std::vector<double>&);
  void Reset(const std::vector<std::vector<double> >&, const std::vector<double>&);
  inline void setNrows(unsigned int v) { nrows = v;};
  inline void setNcols(unsigned int v) { ncols = v;};
  inline unsigned int getNrows(void) const { return nrows; };
  inline unsigned int getNcols(void) const { return ncols; };
};

#endif

Please note that instead of just dying without producing output, the program died most ungracefully with a core dump when I included boost's shared_array. I have no idea why.

Here is the implementation of the regressionPCA class:

#include "reg.class.h"

#include <iostream>
#include <algorithm>
#include <gsl/gsl_linalg.h>

regressionPCA::regressionPCA(const std::vector<std::vector<double> >&data, 
                             const std::vector<double>& Ydata) {
  std::cout << "Flag ALPHA" << std::endl;
  Reset(data,Ydata);
}

void regressionPCA::Reset(const std::vector<std::vector<double> >&data, 
                             const std::vector<double>& Ydata) {
  unsigned int r(0),c(0),n(0);
  std::cout << "r: " << r << "\tc: " << c << "\tn: " << n << std::endl;
  c = data.size();
  std::cout << "r: " << r << "\tc: " << c << "\tn: " << n << std::endl;
  r = data[0].size();
  std::cout << "r: " << r << "\tc: " << c << "\tn: " << n << std::endl;
  setNrows(r);
  setNcols(c);
  n = r * c;
  std::cout << "r: " << r << "\tc: " << c << "\tn: " << n << std::endl;
  double* B = new double[c];
  double* Y = new double[r];
  double* U = new double[n];
  double* V = new double[c*c];
  double* S = new double[c];
  /*  boost::shared_array<double> B(new double[c]), 
  boost::shared_array<double> B(new double[c]), 
                              Y(new double[r]), 
                              U(new double[n]),  
                              V(new double[c*c]),  
                              S(new double[c]);
  B = Btmp;
  Y = Ytmp;
  U = Utmp;
  V = Vtmp;
  S = Stmp;*/
  double *bptr = U.get();
  std::cout << "Flag1" << std::endl;
  std::vector<std::vector<double> >::const_iterator it = data.begin(), end = data.end();
  while (it != end) {
    bptr = std::copy(it->begin(), it->end(),bptr);
    ++it;
  }
  std::cout << "Flag2" << std::endl;
  bptr = Y.get();
  std::copy(Ydata.begin(),Ydata.end(),bptr);
  gsl_matrix_view Um = gsl_matrix_view_array(U.get(), getNrows(), getNcols());
  gsl_matrix_view Ym = gsl_vector_view_array(Y.get(), getNrows());
  gsl_matrix_view Bm = gsl_vector_view_array(B.get(), getNcols());
  gsl_matrix_view Sm = gsl_vector_view_array(S.get(), getNcols());
  gsl_matrix_view Vm = gsl_matrix_view_array(U.get(), getNcols(), getNcols());
  std::cout << "Flag3" << std::endl;
  gsl_linalg_SV_decomp_jacobi(Um,Vm,Sm);
  std::cout << "Flag4" << std::endl;
  gsl_linalg_SV_solve(Um,Vm,Sm,Ym,Bm);
  std::cout << std::endl << std::endl << "Sv = " << gsl_vector_get(Sm,0) << "\t" << gsl_vector_get(Sm,1) << std::endl;
  std::cout << std::endl << std::endl << "V = " << std::endl;
  std::cout << "\t" << gsl_matrix_get(Vm,0,0) << "\t" << gsl_matrix_get(Vm,0,1) << std::endl;
  std::cout << "\t" << gsl_matrix_get(Vm,1,0) << "\t" << gsl_matrix_get(Vm,1,1) << std::endl;
  std::cout << std::endl << std::endl << "Beta = " << gsl_vector_get(Bm,0) << "\t" <<      gsl_vector_get(Bm,1) << std::endl;
}

Note, I know the execution path does not get into the constructor becuase the line "std::cout << "Flag ALPHA" << std::endl;" does not produce output.

So, a couple mysteries.

1) Why does merely including boost's shared_array cause a core dump? 2) when I resort to naked array pointers (yes, I know about the memory leak, but that is the least of my troubles right now, and I would have avoided that had boost's shared array worked for me), why does execution not enter my constructor?

Any insights would be appreciated.

BTW: If you notice me doing something silly with GSL, I'd appreciate knowing about it, even though my problems do not appear to be related to GSL.

Mankarse
  • 39,818
  • 11
  • 97
  • 141
  • +1 This is a fantastically presented question. Welcome to StackOverflow :) – Justin ᚅᚔᚈᚄᚒᚔ Mar 24 '12 at 04:35
  • While debugging by printf is a tried and true technique, sometimes you just have to fire up gdb. – ergosys Mar 24 '12 at 04:52
  • Perhaps better to print traces trhough std::cerr as std::cout may be buffered. – J. Daniel Garcia Mar 24 '12 at 09:59
  • Thanks guys. I tried both, but neither told me anything I didn't already know. We get a segfault because of an access violation, but I don't see how that can happen after Flag C is printed but before the constructor is entered. – user1289485 Mar 24 '12 at 19:59
  • In exasperation, I blew away all the binaries in my development tree, and did a build instead of merely running make. That highlighted some silliness with my use of GSL. I suppose some inconsistency between the object file produced for that class and what function main called caused the segfault, and I suppose that make didn't detect that the source file for that class had been changed (for reasons that remain a mystery). Anyway, thanks for your time, and I will charge ahead debugging this thing. – user1289485 Mar 24 '12 at 21:28

0 Answers0