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,