24 using namespace shogun;
25 using namespace Eigen;
29 m_store_covs(store_covs), m_num_classes(0), m_dim(0)
74 SG_ERROR(
"Specified features are not of type CDotFeatures\n")
88 MatrixXd X(num_vecs, m_dim);
89 MatrixXd A(num_vecs, m_dim);
90 VectorXd norm2(num_vecs*m_num_classes);
96 for (
int k = 0; k < m_num_classes; k++)
99 for (
int i = 0; i < num_vecs; i++)
101 vec = rf->get_feature_vector(i, vlen, vfree);
104 Map< VectorXd > Evec(vec,vlen);
107 X.row(i) = Evec - Em_means_col;
109 rf->free_feature_vector(vec, i, vfree);
112 Map< MatrixXd > Em_M(m_M.
get_matrix(k), m_dim, m_dim);
115 for (
int i = 0; i < num_vecs; i++)
116 norm2(i + k*num_vecs) = A.row(i).array().square().sum();
119 SG_PRINT(
"\n>>> Displaying A ...\n")
124 for (
int i = 0; i < num_vecs; i++)
125 for (
int k = 0; k < m_num_classes; k++)
127 norm2[i + k*num_vecs] += m_slog[k];
128 norm2[i + k*num_vecs] *= -0.5;
132 SG_PRINT(
"\n>>> Displaying norm2 ...\n")
138 for (
int i = 0 ; i < num_vecs; i++)
147 SG_ERROR(
"No labels allocated in QDA training\n")
152 SG_ERROR(
"Speficied features are not of type CDotFeatures\n")
158 SG_ERROR(
"No features allocated in QDA training\n")
163 SG_ERROR(
"No train_labels allocated in QDA training\n")
171 if (num_vec != train_labels.
vlen)
172 SG_ERROR(
"Dimension mismatch between features and labels in QDA training")
174 int32_t* class_idxs = SG_MALLOC(int32_t, num_vec*m_num_classes);
175 int32_t* class_nums = SG_MALLOC(int32_t, m_num_classes);
176 memset(class_nums, 0, m_num_classes*
sizeof(int32_t));
179 for (
int i = 0; i < train_labels.
vlen; i++)
181 class_idx = train_labels.
vector[i];
183 if (class_idx < 0 || class_idx >= m_num_classes)
185 SG_ERROR(
"found label out of {0, 1, 2, ..., num_classes-1}...")
190 class_idxs[ class_idx*num_vec + class_nums[class_idx]++ ] = i;
194 for (
int i = 0; i < m_num_classes; i++)
196 if (class_nums[i] <= 0)
198 SG_ERROR(
"What? One class with no elements\n")
209 cov_dims[2] = m_num_classes;
220 rot_dims[2] = m_num_classes;
230 for (
int k = 0; k < m_num_classes; k++)
232 MatrixXd buffer(class_nums[k], m_dim);
234 for (
int i = 0; i < class_nums[k]; i++)
236 vec = rf->get_feature_vector(class_idxs[k*num_vec + i], vlen, vfree);
239 Map< VectorXd > Evec(vec, vlen);
241 buffer.row(i) = Evec;
243 rf->free_feature_vector(vec, class_idxs[k*num_vec + i], vfree);
246 Em_means /= class_nums[k];
248 for (
int i = 0; i < class_nums[k]; i++)
249 buffer.row(i) -= Em_means;
255 Eigen::JacobiSVD<MatrixXd> eSvd;
256 eSvd.compute(buffer,Eigen::ComputeFullV);
257 memcpy(col, eSvd.singularValues().data(), m_dim*
sizeof(
float64_t));
258 memcpy(rot_mat, eSvd.matrixV().data(), m_dim*m_dim*
sizeof(
float64_t));
267 MatrixXd MEig = Map<MatrixXd>(rot_mat,m_dim,m_dim);
268 for (
int i = 0; i < m_dim; i++)
269 for (
int j = 0; j < m_dim; j++)
270 M(i,j)*=scalings[k*m_dim + j];
271 MatrixXd rotE = Map<MatrixXd>(rot_mat,m_dim,m_dim);
272 MatrixXd resE(m_dim,m_dim);
273 resE = MEig * rotE.transpose();
285 M_dims[2] = m_num_classes;
292 for (
int k = 0; k < m_num_classes; k++)
294 for (
int j = 0; j < m_dim; j++)
296 sinvsqrt[j] = 1.0 /
CMath::sqrt(scalings[k*m_dim + j]);
297 m_slog[k] +=
CMath::log(scalings[k*m_dim + j]);
300 for (
int i = 0; i < m_dim; i++)
301 for (
int j = 0; j < m_dim; j++)
303 idx = k*m_dim*m_dim + i + j*m_dim;
304 m_M[idx] = rotations[idx] * sinvsqrt[j];
309 SG_PRINT(
">>> QDA machine trained with %d classes\n", m_num_classes)
311 SG_PRINT(
"\n>>> Displaying means ...\n")
314 SG_PRINT(
"\n>>> Displaying scalings ...\n")
317 SG_PRINT(
"\n>>> Displaying rotations ... \n")
318 for (
int k = 0; k < m_num_classes; k++)
321 SG_PRINT(
"\n>>> Displaying sinvsqrt ... \n")
324 SG_PRINT(
"\n>>> Diplaying m_M matrices ... \n")
325 for (
int k = 0; k < m_num_classes; k++)