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?