Hi Sven,

 

Thanks for this. 

 

 

Second, in my first post, the problem was due to some unpleasant behaviour of Gnu make I wasn't aware of that resulting in a certain object file not being remade even though the source for it had changed (actually, I don't know if it is in make itself or a bug in my suite of make files - something to be investigated further - doing make realclean followed by make solved that).  My current problem is something else, probably due to how memory is handled either in GSL or in boost::shared_array, but maybe something more subtle.

 

I had left the comments in place to show alternatives I had tried.  It is no surprise that compiling would fail if they were uncommented.

 

I have now changed regressionPCA back to the way I had intended it to be in the first place.

 

Here is the header:

 

#ifndef REG_CLASS_H

#define REG_CLASS_H

 

#include <iosfwd>

#include <vector>

#include <gsl/gsl_linalg.h>

#include <boost/smart_ptr/shared_array.hpp>

 

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

 

And here is the implementation:

 

#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) {

  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;

  r = data.size();

  std::cout << "r: " << r << "\tc: " << c << "\tn: " << n << std::endl;

  c = 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;

  B.reset(new double[c]);

  Y.reset(new double[r]);

  U.reset(new double[n]);

  V.reset(new double[c*c]);

  S.reset(new double[c]);

  double* Btmp = B.get();

  double* Ytmp = Y.get();

  double* Utmp = U.get();

  double* Vtmp = V.get();

  double* Stmp = S.get();

  double *bptr = Utmp;

  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;

  }

  bptr = Ytmp;

  std::copy(Ydata.begin(),Ydata.end(),bptr);

  gsl_matrix_view Um = gsl_matrix_view_array(Utmp, getNrows(), getNcols());

  gsl_vector_view Ym = gsl_vector_view_array(Ytmp, getNrows());

  gsl_vector_view Bm = gsl_vector_view_array(Btmp, getNcols());

  gsl_vector_view Sm = gsl_vector_view_array(Stmp, getNcols());

  gsl_matrix_view Vm = gsl_matrix_view_array(Vtmp, getNcols(), getNcols());

  gsl_linalg_SV_decomp_jacobi(&Um.matrix,&Vm.matrix,&Sm.vector);

  gsl_linalg_SV_solve(&Um.matrix,&Vm.matrix,&Sm.vector,&Ym.vector,&Bm.vector);

  std::cout << std::endl << std::endl << "Sv = " << gsl_vector_get(&Sm.vector,0) << "\t" <<  gsl_vector_get(&Sm.vector,1) << std::endl;

  std::cout << std::endl << std::endl << "V = " << std::endl;

  std::cout << "\t" << gsl_matrix_get(&Vm.matrix,0,0) << "\t" << gsl_matrix_get(&Vm.matrix,0,1) << std::endl;

  std::cout << "\t" << gsl_matrix_get(&Vm.matrix,1,0) << "\t" << gsl_matrix_get(&Vm.matrix,1,1) << std::endl;

  std::cout << std::endl << std::endl << "Beta = " << gsl_vector_get(&Bm.vector,0) << "\t" <<  gsl_vector_get(&Bm.vector,1) << std::endl;

};

 

Obviously a work in progress as I need to create data members to store the result of the analysis, and member functions to pass the results back to the calling code.  I also need to get and analyse the residuals….

 

Here is function main:

 

#include <iostream>

#include <vector>

#include "data.generator.h"

#include "reg.class.h"

#include <boost/smart_ptr/shared_ptr.hpp>

 

int main(int argc, char* argv[]) {

  basicGenerator bg;

  std::cout << "Sample size: " << bg.get_sampleSize() << std::endl;

  bg.makeData();

  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 << "==================================================================" << std::endl;

  regressionPCA rpca(x,y);

  std::cout << "==================================================================" << std::endl;

  boost::shared_ptr<regressionPCA> prpca;

  for (unsigned int j = 0 ; j < 25 ; j++) {

    std::cout << std::endl << std::endl << "Run #: " << (j + 1) << std::endl;

    bg.makeData();

    bg.getDataForRegression(x,y);

    /*    for (unsigned int i = 0 ; i < imax ; i++) {

      std::cout << i << "\t" << y[i] << "\t" << x[i][0] << "\t" << x[i][1] << std::endl;

      }*/

    prpca.reset(new regressionPCA(x,y));

    std::cout << "==================================================================" << std::endl;

  }

  return 0;

}

 

The following output shows that cleaning up the declaration and definition has not changed anything:

 

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

==================================================================

r: 0    c: 0    n: 0

r: 100  c: 0    n: 0

r: 100  c: 2    n: 0

r: 100  c: 2    n: 200

 

Sv = 18.3225    4.69155

 

V =

        0.695362        -0.718659

        0.718659        0.695362

 

Beta = 0.693195 0.627794

==================================================================

Run #: 1

r: 0    c: 0    n: 0

r: 100  c: 0    n: 0

r: 100  c: 2    n: 0

r: 100  c: 2    n: 200

 

Sv = 14.1699    10.7091

 

V =

        0.49497 -0.86891

        0.86891 0.49497

 

Beta = 0.476181 0.391545

==================================================================

 

Run #: 2

r: 0    c: 0    n: 0

r: 100  c: 0    n: 0

r: 100  c: 2    n: 0

r: 100  c: 2    n: 200

 

Sv = 18.3225    4.69155

 

V =

        0.695362        -0.718659

        0.718659        0.695362

 

Beta = 0.693195 0.627794

Aborted (core dumped)

 

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

$

 

There are a couple things I see.

 

1) the last statement successfully executed was the last statement of my Reset function.  It then returns to function main, and it's nect statement produces the line of '+' symbols, and this is not successfully executed.  So, something bad is happening between the end of the function and the return to function main.

2) if I uncomment the inner loop, over I, and comment out ' prpca.reset(new regressionPCA(x,y));', that loop runs to completion.  Thus, we can exclude bg from having any role in this core dump.  It must be something in either gsl or boost::shared_array, or maybe boost::shared_ptr (but the core dump happens too soon to be in shared_ptr's reset, unless the output to std::cout is buffered and contents of the buffer are lost in the core dump before they can be printed, but this is unlikely).

3) It must be something subtle, as I get one execution of the analysis done in Reset, as shown in the first analysis reported by 'rpca', and then I get two complete passes through my loop over j before it crashes.  I am baffled as to even how I determine whether it is something I have done wrong with GSL or with boost::shared_array.

 

But here is an experiment I tried, let's see what you make of this.  I changed my main loop to use:

 

    rpca.Reset(x,y);

 

instead of

 

    prpca.reset(new regressionPCA(x,y));

 

That is, I used the object I'd made on the stack, instead of the one(s) made on the heap.  I then added " std::cout << "Flag 1" << std::endl << std::flush;" right after I invoked boost::shared_array's reset, and " std::cout << "Flag 2" << std::endl << std::flush;" right after I copy new data to these arrays.  Here is the result:

 

98      -0.508617       -0.370597       -0.703225

99      2.69086 2.75267 1.40633

==================================================================

r: 0    c: 0    n: 0

r: 100  c: 0    n: 0

r: 100  c: 2    n: 0

r: 100  c: 2    n: 200

Flag 1

Flag 2

 

 

Sv = 18.3225    4.69155

 

 

V =

        0.695362        -0.718659

        0.718659        0.695362

 

 

Beta = 0.693195 0.627794

==================================================================

 

 

Run #: 1

r: 0    c: 0    n: 0

r: 100  c: 0    n: 0

r: 100  c: 2    n: 0

r: 100  c: 2    n: 200

Flag 1

Flag 2

 

 

Sv = 14.1699    10.7091

 

 

V =

        0.49497 -0.86891

        0.86891 0.49497

 

 

Beta = 0.476181 0.391545

==================================================================

 

 

Run #: 2

r: 0    c: 0    n: 0

r: 100  c: 0    n: 0

r: 100  c: 2    n: 0

r: 100  c: 2    n: 200

      0 [main] test.pca.reg.gsl 7888 exception::handle: Error while dumping state (probably corrupted stack)

Segmentation fault (core dumped)

 

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

$

 

 

Do you notice the change?

 

Instead of finishing the second analysis in the loop, it crashes after I print the dimensions of the input matrix and vector, at some point in resetting one of the boost::shared_array objects.  How do I determine whether this is a bug in boost::shared_array or some bad interaction between boost::shared_array and GSL?

 

 

Cheers

 

Ted

 

From: boost-users-bounces@lists.boost.org [mailto:boost-users-bounces@lists.boost.org] On Behalf Of Sven Steckmann
Sent: March-25-12 4:13 AM
To: boost-users@lists.boost.org
Subject: Re: [Boost-users] Mystery exception.

 

Hi Ted,