Skip to content

Commit

Permalink
Merge pull request #1299 from pickle27/develop
Browse files Browse the repository at this point in the history 
Changed input NDArray to pass by value for AJD class and fixed a memory leak in ICA unit tests and examples
  • Loading branch information
@iglesias
iglesias committed Jul 25, 2013
2 parents 6751620 + 27a4e9e commit da8857a
Show file tree
Hide file tree
Showing 12 changed files with 117 additions and 113 deletions.
2 changes: 2 additions & 0 deletions examples/undocumented/libshogun/converter_sobi_bss.cpp
View file
Expand Up @@ -91,7 +91,9 @@ void test()
std::cout << "Separation error: " << sep_error << std::endl;

SG_UNREF(sobi);
SG_UNREF(mixed_signals);
SG_UNREF(signals);

return;
}

Expand Down
4 changes: 2 additions & 2 deletions src/shogun/mathematics/ajd/ApproxJointDiagonalizer.h
View file
Expand Up @@ -54,8 +54,8 @@ class CApproxJointDiagonalizer : public CSGObject
* @param itermax maximum number of iterations
* @return V the matrix that best diagonalizes C
*/
virtual SGMatrix<float64_t> compute(SGNDArray<float64_t> &C,
SGMatrix<float64_t> V0= SGMatrix<float64_t>(NULL,0,0),
virtual SGMatrix<float64_t> compute(SGNDArray<float64_t> C,
SGMatrix<float64_t> V0 = SGMatrix<float64_t>(NULL,0,0),
double eps=CMath::MACHINE_EPSILON,
int itermax=200) = 0;

Expand Down
2 changes: 1 addition & 1 deletion src/shogun/mathematics/ajd/FFDiag.cpp
View file
Expand Up @@ -16,7 +16,7 @@ using namespace shogun;

void getW(double *C, int *ptN, int *ptK, double *W);

SGMatrix<float64_t> CFFDiag::diagonalize(SGNDArray<float64_t> &C0, SGMatrix<float64_t> V0,
SGMatrix<float64_t> CFFDiag::diagonalize(SGNDArray<float64_t> C0, SGMatrix<float64_t> V0,
double eps, int itermax)
{
int n = C0.dims[0];
Expand Down
8 changes: 4 additions & 4 deletions src/shogun/mathematics/ajd/FFDiag.h
View file
Expand Up @@ -53,8 +53,8 @@ class CFFDiag : public CApproxJointDiagonalizer
* @param itermax maximum number of iterations
* @return V the matrix that best diagonalizes C
*/
static SGMatrix<float64_t> diagonalize(SGNDArray<float64_t> &C,
SGMatrix<float64_t> V0= SGMatrix<float64_t>(NULL,0,0),
static SGMatrix<float64_t> diagonalize(SGNDArray<float64_t> C,
SGMatrix<float64_t> V0 = SGMatrix<float64_t>(NULL,0,0),
double eps=CMath::MACHINE_EPSILON,
int itermax=200);

Expand All @@ -65,8 +65,8 @@ class CFFDiag : public CApproxJointDiagonalizer
* @param itermax maximum number of iterations
* @return V the matrix that best diagonalizes C
*/
virtual SGMatrix<float64_t> compute(SGNDArray<float64_t> &C,
SGMatrix<float64_t> V0= SGMatrix<float64_t>(NULL,0,0),
virtual SGMatrix<float64_t> compute(SGNDArray<float64_t> C,
SGMatrix<float64_t> V0 = SGMatrix<float64_t>(NULL,0,0),
double eps=CMath::MACHINE_EPSILON,
int itermax=200)
{
Expand Down
2 changes: 1 addition & 1 deletion src/shogun/mathematics/ajd/JADiag.cpp
View file
Expand Up @@ -17,7 +17,7 @@ using namespace shogun;
void jadiagw(double c[], double w[], int *ptn, int *ptm, double a[],
double *logdet, double *decr, double *result);

SGMatrix<float64_t> CJADiag::diagonalize(SGNDArray<float64_t> &C, SGMatrix<float64_t> V0,
SGMatrix<float64_t> CJADiag::diagonalize(SGNDArray<float64_t> C, SGMatrix<float64_t> V0,
double eps, int itermax)
{
int d = C.dims[0];
Expand Down
8 changes: 4 additions & 4 deletions src/shogun/mathematics/ajd/JADiag.h
View file
Expand Up @@ -53,8 +53,8 @@ class CJADiag : public CApproxJointDiagonalizer
* @param itermax maximum number of iterations
* @return V the matrix that best diagonalizes C
*/
static SGMatrix<float64_t> diagonalize(SGNDArray<float64_t> &C,
SGMatrix<float64_t> V0= SGMatrix<float64_t>(NULL,0,0),
static SGMatrix<float64_t> diagonalize(SGNDArray<float64_t> C,
SGMatrix<float64_t> V0 = SGMatrix<float64_t>(NULL,0,0),
double eps=CMath::MACHINE_EPSILON,
int itermax=200);

Expand All @@ -65,8 +65,8 @@ class CJADiag : public CApproxJointDiagonalizer
* @param itermax maximum number of iterations
* @return V the matrix that best diagonalizes C
*/
virtual SGMatrix<float64_t> compute(SGNDArray<float64_t> &C,
SGMatrix<float64_t> V0= SGMatrix<float64_t>(NULL,0,0),
virtual SGMatrix<float64_t> compute(SGNDArray<float64_t> C,
SGMatrix<float64_t> V0 = SGMatrix<float64_t>(NULL,0,0),
double eps=CMath::MACHINE_EPSILON,
int itermax=200)
{
Expand Down
184 changes: 92 additions & 92 deletions src/shogun/mathematics/ajd/JADiagOrth.cpp
View file
Expand Up @@ -19,11 +19,11 @@ void left_rot_stack(double *A, int M, int N, int K, int p, int q, double c, doub
void right_rot_stack(double *A, int M, int N, int K, int p, int q, double c, double s);
void left_rot_simple(double *A, int m, int n, int p, int q, double c, double s);

SGMatrix<float64_t> CJADiagOrth::diagonalize(SGNDArray<float64_t> &C, SGMatrix<float64_t> V0,
SGMatrix<float64_t> CJADiagOrth::diagonalize(SGNDArray<float64_t> C, SGMatrix<float64_t> V0,
double eps, int itermax)
{
{
int m = C.dims[0];
int L = C.dims[2];
int L = C.dims[2];

SGMatrix<float64_t> V;
if (V0.num_rows != 0)
Expand All @@ -40,28 +40,28 @@ SGMatrix<float64_t> CJADiagOrth::diagonalize(SGNDArray<float64_t> &C, SGMatrix<f

while (more)
{
more = false;

for (int p = 0; p < m; p++)
{
for (int q = p+1; q < m; q++)
{
// computation of Givens angle
double theta = givens_stack(C.array, m, L, p, q);

// Givens update
if (fabs(theta) > eps)
{
double c = cos(theta);
double s = sin(theta);
left_rot_stack (C.array, m, m, L, p, q, c, s);
right_rot_stack(C.array, m, m, L, p, q, c, s);
left_rot_simple(V.matrix, m, m, p, q, c, s);
rots++;
more = true;
}
}
}
more = false;

for (int p = 0; p < m; p++)
{
for (int q = p+1; q < m; q++)
{
// computation of Givens angle
double theta = givens_stack(C.array, m, L, p, q);

// Givens update
if (fabs(theta) > eps)
{
double c = cos(theta);
double s = sin(theta);
left_rot_stack (C.array, m, m, L, p, q, c, s);
right_rot_stack(C.array, m, m, L, p, q, c, s);
left_rot_simple(V.matrix, m, m, p, q, c, s);
rots++;
more = true;
}
}
}
}

return V;
Expand All @@ -70,33 +70,33 @@ SGMatrix<float64_t> CJADiagOrth::diagonalize(SGNDArray<float64_t> &C, SGMatrix<f
/* Givens angle for the pair (p,q) of a stack of K M*M matrices */
double givens_stack(double *A, int M, int K, int p, int q)
{
int k;
double diff_on, sum_off, ton, toff;
double *cm; // A cumulant matrix
double G11 = 0.0;
double G12 = 0.0;
double G22 = 0.0;

int M2 = M*M;
int pp = p+p*M;
int pq = p+q*M;
int qp = q+p*M;
int qq = q+q*M;

for (k=0, cm=A; k<K; k++, cm+=M2)
{
diff_on = cm[pp] - cm[qq];
sum_off = cm[pq] + cm[qp];

G11 += diff_on * diff_on;
G22 += sum_off * sum_off;
G12 += diff_on * sum_off;
}
int k;
double diff_on, sum_off, ton, toff;
double *cm; // A cumulant matrix
double G11 = 0.0;
double G12 = 0.0;
double G22 = 0.0;

int M2 = M*M;
int pp = p+p*M;
int pq = p+q*M;
int qp = q+p*M;
int qq = q+q*M;

for (k=0, cm=A; k<K; k++, cm+=M2)
{
diff_on = cm[pp] - cm[qq];
sum_off = cm[pq] + cm[qp];

G11 += diff_on * diff_on;
G22 += sum_off * sum_off;
G12 += diff_on * sum_off;
}

ton = G11 - G22;
toff = 2.0 * G12;
ton = G11 - G22;
toff = 2.0 * G12;

return -0.5 * atan2 ( toff , ton+sqrt(ton*ton+toff*toff) );
return -0.5 * atan2 ( toff , ton+sqrt(ton*ton+toff*toff) );
}

/*
Expand All @@ -105,60 +105,60 @@ double givens_stack(double *A, int M, int K, int p, int q)
*/
void left_rot_stack(double *A, int M, int N, int K, int p, int q, double c, double s )
{
int k, ix, iy, cpt;
int MN = M*N;
int kMN;
double nx, ny;

for (k=0, kMN=0; k<K; k++, kMN+=MN)
{
for (cpt=0, ix=p+kMN, iy=q+kMN; cpt<N; cpt++, ix+=M, iy+=M)
{
nx = A[ix];
ny = A[iy];
A[ix] = c*nx - s*ny;
A[iy] = s*nx + c*ny;
}
}
int k, ix, iy, cpt;
int MN = M*N;
int kMN;
double nx, ny;

for (k=0, kMN=0; k<K; k++, kMN+=MN)
{
for (cpt=0, ix=p+kMN, iy=q+kMN; cpt<N; cpt++, ix+=M, iy+=M)
{
nx = A[ix];
ny = A[iy];
A[ix] = c*nx - s*ny;
A[iy] = s*nx + c*ny;
}
}
}

/* Ak(mxn) --> Ak(mxn) x R where R rotates the (p,q) columns R =[ c s ; -s c ]
and Ak is the k-th M*N matrix in the stack */
void right_rot_stack(double *A, int M, int N, int K, int p, int q, double c, double s )
{
int k, ix, iy, cpt, kMN;
int pM = p*M;
int qM = q*M;
double nx, ny;

for (k=0, kMN=0; k<K; k++, kMN+=M*N)
{
for (cpt=0, ix=pM+kMN, iy=qM+kMN; cpt<M; cpt++)
{
nx = A[ix];
ny = A[iy];
A[ix++] = c*nx - s*ny;
A[iy++] = s*nx + c*ny;
}
}
int k, ix, iy, cpt, kMN;
int pM = p*M;
int qM = q*M;
double nx, ny;

for (k=0, kMN=0; k<K; k++, kMN+=M*N)
{
for (cpt=0, ix=pM+kMN, iy=qM+kMN; cpt<M; cpt++)
{
nx = A[ix];
ny = A[iy];
A[ix++] = c*nx - s*ny;
A[iy++] = s*nx + c*ny;
}
}
}

/*
A(mxn) --> R * A(mxn) where R=[ c -s ; s c ] rotates the (p,q) rows of R
*/
void left_rot_simple(double *A, int m, int n, int p, int q, double c, double s)
{
int ix = p;
int iy = q;
double nx, ny;
int j;

for (j=0; j<n; j++, ix+=m, iy+=m)
{
nx = A[ix];
ny = A[iy];
A[ix] = c*nx - s*ny;
A[iy] = s*nx + c*ny;
}
int ix = p;
int iy = q;
double nx, ny;
int j;

for (j=0; j<n; j++, ix+=m, iy+=m)
{
nx = A[ix];
ny = A[iy];
A[ix] = c*nx - s*ny;
A[iy] = s*nx + c*ny;
}
}
#endif //HAVE_EIGEN3
8 changes: 4 additions & 4 deletions src/shogun/mathematics/ajd/JADiagOrth.h
View file
Expand Up @@ -53,8 +53,8 @@ class CJADiagOrth : public CApproxJointDiagonalizer
* @param itermax maximum number of iterations
* @return V the matrix that best diagonalizes C
*/
static SGMatrix<float64_t> diagonalize(SGNDArray<float64_t> &C,
SGMatrix<float64_t> V0= SGMatrix<float64_t>(NULL,0,0),
static SGMatrix<float64_t> diagonalize(SGNDArray<float64_t> C,
SGMatrix<float64_t> V0 = SGMatrix<float64_t>(NULL,0,0),
double eps=CMath::MACHINE_EPSILON,
int itermax=200);

Expand All @@ -65,8 +65,8 @@ class CJADiagOrth : public CApproxJointDiagonalizer
* @param itermax maximum number of iterations
* @return V the matrix that best diagonalizes C
*/
virtual SGMatrix<float64_t> compute(SGNDArray<float64_t> &C,
SGMatrix<float64_t> V0= SGMatrix<float64_t>(NULL,0,0),
virtual SGMatrix<float64_t> compute(SGNDArray<float64_t> C,
SGMatrix<float64_t> V0 = SGMatrix<float64_t>(NULL,0,0),
double eps=CMath::MACHINE_EPSILON,
int itermax=200)
{
Expand Down
2 changes: 1 addition & 1 deletion src/shogun/mathematics/ajd/UWedge.cpp
View file
Expand Up @@ -14,7 +14,7 @@ typedef Matrix< float64_t, Dynamic, Dynamic, ColMajor > EMatrix;

using namespace shogun;

SGMatrix<float64_t> CUWedge::diagonalize(SGNDArray<float64_t> &C, SGMatrix<float64_t> V0,
SGMatrix<float64_t> CUWedge::diagonalize(SGNDArray<float64_t> C, SGMatrix<float64_t> V0,
double eps, int itermax)
{
int d = C.dims[0];
Expand Down
8 changes: 4 additions & 4 deletions src/shogun/mathematics/ajd/UWedge.h
View file
Expand Up @@ -52,8 +52,8 @@ class CUWedge : public CApproxJointDiagonalizer
* @param itermax maximum number of iterations
* @return V the matrix that best diagonalizes C
*/
static SGMatrix<float64_t> diagonalize(SGNDArray<float64_t> &C,
SGMatrix<float64_t> V0= SGMatrix<float64_t>(NULL,0,0),
static SGMatrix<float64_t> diagonalize(SGNDArray<float64_t> C,
SGMatrix<float64_t> V0 = SGMatrix<float64_t>(NULL,0,0),
double eps=1e-12,
int itermax=200);

Expand All @@ -64,8 +64,8 @@ class CUWedge : public CApproxJointDiagonalizer
* @param itermax maximum number of iterations
* @return V the matrix that best diagonalizes C
*/
virtual SGMatrix<float64_t> compute(SGNDArray<float64_t> &C,
SGMatrix<float64_t> V0= SGMatrix<float64_t>(NULL,0,0),
virtual SGMatrix<float64_t> compute(SGNDArray<float64_t> C,
SGMatrix<float64_t> V0 = SGMatrix<float64_t>(NULL,0,0),
double eps=1e-12,
int itermax=200)
{
Expand Down
1 change: 1 addition & 0 deletions tests/unit/converter/ica/Jade_unittest.cc
View file
Expand Up @@ -61,6 +61,7 @@ TEST(CJade, blind_source_separation)
EXPECT_EQ(isperm,true);

SG_UNREF(jade);
SG_UNREF(mixed_signals);
SG_UNREF(signals);
}

Expand Down
1 change: 1 addition & 0 deletions tests/unit/converter/ica/SOBI_unittest.cc
View file
Expand Up @@ -56,6 +56,7 @@ TEST(CSOBI, blind_source_separation)
EXPECT_LE(sep_error, 1e-6);

SG_UNREF(sobi);
SG_UNREF(mixed_signals);
SG_UNREF(signals);
}

Expand Down

0 comments on commit da8857a

Please sign in to comment.