diff --git a/examples/undocumented/libshogun/classifier_lda.cpp b/examples/undocumented/libshogun/classifier_lda.cpp index 7246f5ef8de..d4c4ba81981 100644 --- a/examples/undocumented/libshogun/classifier_lda.cpp +++ b/examples/undocumented/libshogun/classifier_lda.cpp @@ -23,7 +23,7 @@ using namespace shogun; -#define NUM 50 +#define NUM 50 #define DIMS 2 #define CLASSES 2 diff --git a/examples/undocumented/libshogun/classifier_qda.cpp b/examples/undocumented/libshogun/classifier_qda.cpp index 488a371877f..f81e2fdea88 100644 --- a/examples/undocumented/libshogun/classifier_qda.cpp +++ b/examples/undocumented/libshogun/classifier_qda.cpp @@ -18,12 +18,13 @@ using namespace shogun; -#define NUM 50 +#define NUM 50 #define DIMS 2 #define CLASSES 2 void test() { +#ifdef HAVE_EIGEN3 SGVector< float64_t > lab(CLASSES*NUM); SGMatrix< float64_t > feat(DIMS, CLASSES*NUM); @@ -51,6 +52,7 @@ void test() // Free memory SG_UNREF(output); SG_UNREF(qda); +#endif } int main(int argc, char ** argv) diff --git a/src/shogun/multiclass/QDA.cpp b/src/shogun/multiclass/QDA.cpp index 0f821f4f4c2..8de189ba287 100644 --- a/src/shogun/multiclass/QDA.cpp +++ b/src/shogun/multiclass/QDA.cpp @@ -10,7 +10,7 @@ #include -#ifdef HAVE_LAPACK +#ifdef HAVE_EIGEN3 #include #include @@ -18,7 +18,14 @@ #include #include #include -#include + +#include + +using namespace Eigen; + +typedef Matrix< float64_t, Dynamic, Dynamic, ColMajor > EMatrix; +typedef Matrix< float64_t, Dynamic, 1, ColMajor > EVector; +typedef Array< float64_t, Dynamic, 1 > EArray; using namespace shogun; @@ -83,44 +90,44 @@ CMulticlassLabels* CQDA::apply_multiclass(CFeatures* data) CDenseFeatures< float64_t >* rf = (CDenseFeatures< float64_t >*) m_features; - SGMatrix< float64_t > X(num_vecs, m_dim); - SGMatrix< float64_t > A(num_vecs, m_dim); - SGVector< float64_t > norm2(num_vecs*m_num_classes); - norm2.zero(); + EMatrix X(num_vecs, m_dim); + EMatrix A(num_vecs, m_dim); + EVector norm2(num_vecs*m_num_classes); + norm2.setZero(); - int i, j, k, vlen; + int32_t vlen; bool vfree; float64_t* vec; - for ( k = 0 ; k < m_num_classes ; ++k ) + for (int k = 0; k < m_num_classes; k++) { // X = features - means - for ( i = 0 ; i < num_vecs ; ++i ) + for (int i = 0; i < num_vecs; i++) { vec = rf->get_feature_vector(i, vlen, vfree); ASSERT(vec) - for ( j = 0 ; j < m_dim ; ++j ) - X[i + j*num_vecs] = vec[j] - m_means[k*m_dim + j]; + Eigen::Map< EVector > Evec(vec,vlen); + Eigen::Map< EVector > Em_means_col(m_means.get_column_vector(k), m_dim); - rf->free_feature_vector(vec, i, vfree); + X.row(i) = Evec - Em_means_col; + rf->free_feature_vector(vec, i, vfree); } - cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, num_vecs, m_dim, - m_dim, 1.0, X.matrix, num_vecs, m_M.get_matrix(k), m_dim, 0.0, - A.matrix, num_vecs); + Eigen::Map< EMatrix > Em_M(m_M.get_matrix(k), m_dim, m_dim); + A = X*Em_M; - for ( i = 0 ; i < num_vecs ; ++i ) - for ( j = 0 ; j < m_dim ; ++j ) - norm2[i + k*num_vecs] += CMath::sq(A[i + j*num_vecs]); + for (int i = 0; i < num_vecs; i++) + norm2(i + k*num_vecs) = A.row(i).array().square().sum(); #ifdef DEBUG_QDA - CMath::display_matrix(A.matrix, num_vecs, m_dim, "A"); + SG_PRINT("\n>>> Displaying A ...\n") + SGMatrix< float64_t >::display_matrix(A.data(), num_vecs, m_dim); #endif } - for ( i = 0 ; i < num_vecs ; ++i ) - for ( k = 0 ; k < m_num_classes ; ++k ) + for (int i = 0; i < num_vecs; i++) + for (int k = 0; k < m_num_classes; k++) { norm2[i + k*num_vecs] += m_slog[k]; norm2[i + k*num_vecs] *= -0.5; @@ -128,12 +135,14 @@ CMulticlassLabels* CQDA::apply_multiclass(CFeatures* data) CMulticlassLabels* out = new CMulticlassLabels(num_vecs); - for ( i = 0 ; i < num_vecs ; ++i ) - out->set_label(i, SGVector::arg_max(norm2.vector+i, num_vecs, m_num_classes)); + for (int i = 0 ; i < num_vecs; i++) + out->set_label(i, SGVector::arg_max(norm2.data()+i, num_vecs, m_num_classes)); #ifdef DEBUG_QDA - CMath::display_matrix(norm2.vector, num_vecs, m_num_classes, "norm2"); - CMath::display_vector(out->get_labels().vector, num_vecs, "Labels"); + SG_PRINT("\n>>> Displaying norm2 ...\n") + SGMatrix< float64_t >::display_matrix(norm2.data(), num_vecs, m_num_classes); + SG_PRINT("\n>>> Displaying out ...\n") + SGVector< float64_t >::display_vector(out->get_labels().vector, num_vecs); #endif return out; @@ -141,19 +150,23 @@ CMulticlassLabels* CQDA::apply_multiclass(CFeatures* data) bool CQDA::train_machine(CFeatures* data) { - if ( !m_labels ) + if (!m_labels) SG_ERROR("No labels allocated in QDA training\n") if ( data ) { - if ( !data->has_property(FP_DOT) ) + if (!data->has_property(FP_DOT)) SG_ERROR("Speficied features are not of type CDotFeatures\n") + set_features((CDotFeatures*) data); } - if ( !m_features ) + + if (!m_features) SG_ERROR("No features allocated in QDA training\n") + SGVector< int32_t > train_labels = ((CMulticlassLabels*) m_labels)->get_int_labels(); - if ( !train_labels.vector ) + + if (!train_labels.vector) SG_ERROR("No train_labels allocated in QDA training\n") cleanup(); @@ -161,20 +174,20 @@ bool CQDA::train_machine(CFeatures* data) m_num_classes = ((CMulticlassLabels*) m_labels)->get_num_classes(); m_dim = m_features->get_dim_feature_space(); int32_t num_vec = m_features->get_num_vectors(); - if ( num_vec != train_labels.vlen ) + + if (num_vec != train_labels.vlen) SG_ERROR("Dimension mismatch between features and labels in QDA training") - int32_t* class_idxs = SG_MALLOC(int32_t, num_vec*m_num_classes); - // number of examples of each class + int32_t* class_idxs = SG_MALLOC(int32_t, num_vec*m_num_classes); // number of examples of each class int32_t* class_nums = SG_MALLOC(int32_t, m_num_classes); memset(class_nums, 0, m_num_classes*sizeof(int32_t)); int32_t class_idx; - int32_t i, j, k; - for ( i = 0 ; i < train_labels.vlen ; ++i ) + + for (int i = 0; i < train_labels.vlen; i++) { class_idx = train_labels.vector[i]; - if ( class_idx < 0 || class_idx >= m_num_classes ) + if (class_idx < 0 || class_idx >= m_num_classes) { SG_ERROR("found label out of {0, 1, 2, ..., num_classes-1}...") return false; @@ -185,16 +198,16 @@ bool CQDA::train_machine(CFeatures* data) } } - for ( i = 0 ; i < m_num_classes ; ++i ) + for (int i = 0; i < m_num_classes; i++) { - if ( class_nums[i] <= 0 ) + if (class_nums[i] <= 0) { SG_ERROR("What? One class with no elements\n") return false; } } - if ( m_store_covs ) + if (m_store_covs) { // cov_dims will be free in m_covs.destroy_ndarray() index_t * cov_dims = SG_MALLOC(index_t, 3); @@ -221,63 +234,55 @@ bool CQDA::train_machine(CFeatures* data) int32_t vlen; bool vfree; float64_t* vec; - for ( k = 0 ; k < m_num_classes ; ++k ) + for (int k = 0; k < m_num_classes; k++) { - SGMatrix< float64_t > buffer(class_nums[k], m_dim); - for ( i = 0 ; i < class_nums[k] ; ++i ) + EMatrix buffer(class_nums[k], m_dim); + Eigen::Map< EVector > Em_means(m_means.get_column_vector(k), m_dim); + for (int i = 0; i < class_nums[k]; i++) { vec = rf->get_feature_vector(class_idxs[k*num_vec + i], vlen, vfree); ASSERT(vec) - for ( j = 0 ; j < vlen ; ++j ) - { - m_means[k*m_dim + j] += vec[j]; - buffer[i + j*class_nums[k]] = vec[j]; - } + Eigen::Map< EVector > Evec(vec, vlen); + Em_means += Evec; + buffer.row(i) = Evec; rf->free_feature_vector(vec, class_idxs[k*num_vec + i], vfree); } - for ( j = 0 ; j < m_dim ; ++j ) - m_means[k*m_dim + j] /= class_nums[k]; + Em_means /= class_nums[k]; - for ( i = 0 ; i < class_nums[k] ; ++i ) - for ( j = 0 ; j < m_dim ; ++j ) - buffer[i + j*class_nums[k]] -= m_means[k*m_dim + j]; + for (int i = 0; i < class_nums[k]; i++) + buffer.row(i) -= Em_means; - /* calling external lib, buffer = U * S * V^T, U is not interesting here */ - char jobu = 'N', jobvt = 'A'; - int m = class_nums[k], n = m_dim; - int lda = m, ldu = m, ldvt = n; - int info = -1; + // SVD float64_t * col = scalings.get_column_vector(k); float64_t * rot_mat = rotations.get_matrix(k); - wrap_dgesvd(jobu, jobvt, m, n, buffer.matrix, lda, col, NULL, ldu, - rot_mat, ldvt, &info); - ASSERT(info == 0) - buffer=SGMatrix(); + Eigen::JacobiSVD eSvd; + eSvd.compute(buffer,Eigen::ComputeFullV); + memcpy(col, eSvd.singularValues().data(), m_dim*sizeof(float64_t)); + memcpy(rot_mat, eSvd.matrixV().data(), m_dim*m_dim*sizeof(float64_t)); SGVector::vector_multiply(col, col, col, m_dim); - SGVector::scale_vector(1.0/(m-1), col, m_dim); + SGVector::scale_vector(1.0/(class_nums[k]-1), col, m_dim); rotations.transpose_matrix(k); - if ( m_store_covs ) + if (m_store_covs) { - SGMatrix< float64_t > M(n ,n); + Eigen::Map< EMatrix > EM(SGVector::clone_vector(rot_mat, m_dim*m_dim), m_dim, m_dim); + Eigen::Map< EArray > Escalings(scalings.get_column_vector(k), m_dim); + for (int i = 0; i < m_dim; i++) + EM.row(i) = ( (EM.row(i).array()) * Escalings ).matrix(); - M.matrix = SGVector::clone_vector(rot_mat, n*n); - for ( i = 0 ; i < m_dim ; ++i ) - for ( j = 0 ; j < m_dim ; ++j ) - M[i + j*m_dim] *= scalings[k*m_dim + j]; + Eigen::Map< EMatrix > Em_covs(m_covs.get_matrix(k), m_dim, m_dim); + Eigen::Map< EMatrix > Erot_mat(rot_mat, m_dim, m_dim); - cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans, n, n, n, 1.0, - M.matrix, n, rot_mat, n, 0.0, m_covs.get_matrix(k), n); + Em_covs = EM*Erot_mat; } } /* Computation of terms required for classification */ - SGVector< float32_t > sinvsqrt(m_dim); // M_dims will be freed in m_M.destroy_ndarray() @@ -291,16 +296,16 @@ bool CQDA::train_machine(CFeatures* data) m_slog.zero(); index_t idx = 0; - for ( k = 0 ; k < m_num_classes ; ++k ) + for (int k = 0; k < m_num_classes; k++) { - for ( j = 0 ; j < m_dim ; ++j ) + for (int j = 0; j < m_dim; j++) { sinvsqrt[j] = 1.0 / CMath::sqrt(scalings[k*m_dim + j]); m_slog[k] += CMath::log(scalings[k*m_dim + j]); } - for ( i = 0 ; i < m_dim ; ++i ) - for ( j = 0 ; j < m_dim ; ++j ) + for (int i = 0; i < m_dim; i++) + for (int j = 0; j < m_dim; j++) { idx = k*m_dim*m_dim + i + j*m_dim; m_M[idx] = rotations[idx] * sinvsqrt[j]; @@ -311,21 +316,21 @@ bool CQDA::train_machine(CFeatures* data) SG_PRINT(">>> QDA machine trained with %d classes\n", m_num_classes) SG_PRINT("\n>>> Displaying means ...\n") - CMath::display_matrix(m_means.matrix, m_dim, m_num_classes); + SGMatrix< float64_t >::display_matrix(m_means.matrix, m_dim, m_num_classes); SG_PRINT("\n>>> Displaying scalings ...\n") - CMath::display_matrix(scalings.matrix, m_dim, m_num_classes); + SGMatrix< float64_t >::display_matrix(scalings.matrix, m_dim, m_num_classes); SG_PRINT("\n>>> Displaying rotations ... \n") - for ( k = 0 ; k < m_num_classes ; ++k ) - CMath::display_matrix(rotations.get_matrix(k), m_dim, m_dim); + for (int k = 0; k < m_num_classes; k++) + SGMatrix< float64_t >::display_matrix(rotations.get_matrix(k), m_dim, m_dim); SG_PRINT("\n>>> Displaying sinvsqrt ... \n") sinvsqrt.display_vector(); SG_PRINT("\n>>> Diplaying m_M matrices ... \n") - for ( k = 0 ; k < m_num_classes ; ++k ) - CMath::display_matrix(m_M.get_matrix(k), m_dim, m_dim); + for (int k = 0; k < m_num_classes; k++) + SGMatrix< float64_t >::display_matrix(m_M.get_matrix(k), m_dim, m_dim); SG_PRINT("\n>>> Exit DEBUG_QDA\n") #endif @@ -335,4 +340,4 @@ bool CQDA::train_machine(CFeatures* data) return true; } -#endif /* HAVE_LAPACK */ +#endif /* HAVE_EIGEN3 */ diff --git a/src/shogun/multiclass/QDA.h b/src/shogun/multiclass/QDA.h index 446a83088a2..9667f0d217d 100644 --- a/src/shogun/multiclass/QDA.h +++ b/src/shogun/multiclass/QDA.h @@ -13,7 +13,7 @@ #include -#ifdef HAVE_LAPACK +#ifdef HAVE_EIGEN3 #include #include @@ -189,5 +189,5 @@ class CQDA : public CNativeMulticlassMachine }; /* class QDA */ } /* namespace shogun */ -#endif /* HAVE_LAPACK */ +#endif /* HAVE_EIGEN3 */ #endif /* _QDA_H__ */ diff --git a/tests/unit/multiclass/MCLDA_unittest.cc b/tests/unit/multiclass/MCLDA_unittest.cc index 8a09ec58301..d87f747b74a 100644 --- a/tests/unit/multiclass/MCLDA_unittest.cc +++ b/tests/unit/multiclass/MCLDA_unittest.cc @@ -5,7 +5,7 @@ #include #ifdef HAVE_EIGEN3 -#define NUM 50 +#define NUM 50 #define DIMS 2 #define CLASSES 2 @@ -30,7 +30,7 @@ TEST(MCLDA, train_and_apply) CMulticlassLabels* output=CLabelsFactory::to_multiclass(lda->apply()); SG_REF(output); - + // Test for ( index_t i = 0; i < CLASSES*NUM; ++i ) EXPECT_EQ(output->get_label(i), labels->get_label(i)); diff --git a/tests/unit/multiclass/QDA_unittest.cc b/tests/unit/multiclass/QDA_unittest.cc new file mode 100644 index 00000000000..3d28531ad38 --- /dev/null +++ b/tests/unit/multiclass/QDA_unittest.cc @@ -0,0 +1,41 @@ +#include +#include +#include +#include +#include +#ifdef HAVE_EIGEN3 + +#define NUM 50 +#define DIMS 2 +#define CLASSES 2 + +using namespace shogun; + +TEST(QDA, train_and_apply) +{ + SGVector< float64_t > lab(CLASSES*NUM); + SGMatrix< float64_t > feat(DIMS, CLASSES*NUM); + + feat = CDataGenerator::generate_gaussians(NUM,CLASSES,DIMS); + for( int i = 0 ; i < CLASSES ; ++i ) + for( int j = 0 ; j < NUM ; ++j ) + lab[i*NUM+j] = double(i); + + CMulticlassLabels* labels = new CMulticlassLabels(lab); + CDenseFeatures< float64_t >* features = new CDenseFeatures< float64_t >(feat); + + CQDA* qda = new CQDA(features, labels); + SG_REF(qda); + qda->train(); + + CMulticlassLabels* output=CLabelsFactory::to_multiclass(qda->apply()); + SG_REF(output); + + // Test + for ( index_t i = 0; i < CLASSES*NUM; ++i ) + EXPECT_EQ(output->get_label(i), labels->get_label(i)); + + SG_UNREF(output); + SG_UNREF(qda); +} +#endif //HAVE_EIGEN3