mog_diag_kmeans.cpp
Go to the documentation of this file.
#include <itpp/stat/mog_diag_kmeans.h>
#include <iostream>
namespace itpp
{
inline double MOG_diag_kmeans_sup::dist(const double * x, const double * y) const
{
double acc = 0.0;
for (int d = 0;d < D;d++) { double tmp = x[d] - y[d]; acc += tmp * tmp; }
return(acc);
}
void MOG_diag_kmeans_sup::assign_to_means()
{
for (int k = 0;k < K;k++) c_count[k] = 0;
for (int n = 0;n < N;n++) {
int k = 0;
double min_dist = dist(c_means[k], c_X[n]);
int k_winner = k;
for (int k = 1;k < K;k++) {
double tmp_dist = dist(c_means[k], c_X[n]);
if (tmp_dist < min_dist) { min_dist = tmp_dist; k_winner = k; }
}
c_partitions[ k_winner ][ count[k_winner] ] = n;
c_count[k_winner]++;
}
}
void MOG_diag_kmeans_sup::recalculate_means()
{
for (int k = 0;k < K;k++) {
for (int d = 0;d < D;d++) c_tmpvec[d] = 0.0;
int Nk = c_count[k];
for (int n = 0;n < Nk;n++) {
double * x = c_X[ c_partitions[k][n] ];
for (int d = 0;d < D;d++) c_tmpvec[d] += x[d];
}
if (Nk > 0) {
double * c_mean = c_means[k];
for (int d = 0;d < D;d++) c_mean[d] = c_tmpvec[d] / Nk;
}
}
}
bool MOG_diag_kmeans_sup::dezombify_means()
{
static int counter = 0;
bool zombie_mean = false;
int k = 0;
int max_count = count[k];
int k_hog = k;
for (int k = 1;k < K;k++) if (c_count[k] > max_count) { max_count = c_count[k]; k_hog = k; }
for (int k = 0;k < K;k++) {
if (c_count[k] == 0) {
zombie_mean = true;
if (verbose) {
it_warning("MOG_diag_kmeans_sup::dezombify_means(): detected zombie mean");
}
if (k_hog == k) {
it_warning("MOG_diag_kmeans_sup::dezombify_means(): weirdness: k_hog == k");
return(false);
}
if (counter >= c_count[k_hog]) counter = 0;
double * c_mean = c_means[k];
double * c_x = c_X[ c_partitions[k_hog][counter] ];
for (int d = 0;d < D;d++) c_mean[d] = 0.5 * (c_means[k_hog][d] + c_x[d]);
counter++;
}
}
if (zombie_mean) assign_to_means();
return(true);
}
double MOG_diag_kmeans_sup::measure_change() const
{
double tmp_dist = 0.0;
for (int k = 0;k < K;k++) tmp_dist += dist(c_means[k], c_means_old[k]);
return(tmp_dist);
}
void MOG_diag_kmeans_sup::initial_means()
{
for (int d = 0;d < D;d++) c_tmpvec[d] = 0.0;
for (int n = 0;n < N;n++) {
double * c_x = c_X[n];
for (int d = 0;d < D;d++) c_tmpvec[d] += c_x[d];
}
for (int d = 0;d < D;d++) c_tmpvec[d] /= N;
int step = int(floor(double(N) / double(K)));
for (int k = 0;k < K;k++) {
double * c_mean = c_means[k];
double * c_x = c_X[k*step];
for (int d = 0;d < D;d++) c_mean[d] = 0.5 * (c_tmpvec[d] + c_x[d]);
}
}
void MOG_diag_kmeans_sup::iterate()
{
for (int k = 0;k < K;k++) for (int d = 0;d < D;d++) c_means_old[k][d] = c_means[k][d];
for (int i = 0;i < max_iter;i++) {
assign_to_means();
if (!dezombify_means()) return;
recalculate_means();
double change = measure_change();
if (verbose) std::cout << "MOG_diag_kmeans_sup::iterate(): iteration = " << i << " change = " << change << std::endl;
if (change == 0) break;
for (int k = 0;k < K;k++) for (int d = 0;d < D;d++) c_means_old[k][d] = c_means[k][d];
}
}
void MOG_diag_kmeans_sup::calc_means()
{
initial_means();
iterate();
}
void MOG_diag_kmeans_sup::calc_covs()
{
for (int k = 0;k < K;k++) {
int Nk = c_count[k];
if (Nk >= 2) {
double * c_mean = c_means[k];
for (int d = 0;d < D;d++) c_tmpvec[d] = 0.0;
for (int n = 0;n < Nk;n++) {
double * c_x = c_X[ c_partitions[k][n] ];
for (int d = 0;d < D;d++) { double tmp = c_x[d] - c_mean[d]; c_tmpvec[d] += tmp * tmp; }
}
for (int d = 0;d < D;d++) c_diag_covs[k][d] = trust * (c_tmpvec[d] / (Nk - 1.0)) + (1.0 - trust) * (1.0);
}
else {
for (int d = 0;d < D;d++) c_diag_covs[k][d] = 1.0;
}
}
}
void MOG_diag_kmeans_sup::calc_weights()
{
for (int k = 0;k < K;k++) c_weights[k] = trust * (c_count[k] / double(N)) + (1.0 - trust) * (1.0 / K);
}
void MOG_diag_kmeans_sup::normalise_vectors()
{
for (int d = 0;d < D;d++) {
double acc = 0.0;
for (int n = 0;n < N;n++) acc += c_X[n][d];
c_norm_mu[d] = acc / N;
}
for (int d = 0;d < D;d++) {
double acc = 0.0;
for (int n = 0;n < N;n++) { double tmp = c_X[n][d] - c_norm_mu[d]; acc += tmp * tmp; }
c_norm_sd[d] = std::sqrt(acc / (N - 1));
}
for (int n = 0;n < N;n++) for (int d = 0;d < D;d++) {
c_X[n][d] -= c_norm_mu[d];
if (c_norm_sd[d] > 0.0) c_X[n][d] /= c_norm_sd[d];
}
}
void MOG_diag_kmeans_sup::unnormalise_vectors()
{
for (int n = 0;n < N;n++) for (int d = 0;d < D;d++) {
if (c_norm_sd[d] > 0.0) c_X[n][d] *= c_norm_sd[d];
c_X[n][d] += c_norm_mu[d];
}
}
void MOG_diag_kmeans_sup::unnormalise_means()
{
for (int k = 0;k < K;k++) for (int d = 0;d < D;d++) {
if (norm_sd[d] > 0.0) c_means[k][d] *= c_norm_sd[d];
c_means[k][d] += norm_mu[d];
}
}
void MOG_diag_kmeans_sup::run(MOG_diag &model_in, Array<vec> &X_in, int max_iter_in, double trust_in, bool normalise_in, bool verbose_in)
{
it_assert(model_in.is_valid(), "MOG_diag_kmeans_sup::run(): given model is not valid");
it_assert((max_iter_in > 0), "MOG_diag_kmeans_sup::run(): 'max_iter' needs to be greater than zero");
it_assert(((trust_in >= 0.0) && (trust_in <= 1.0)), "MOG_diag_kmeans_sup::run(): 'trust' must be between 0 and 1 (inclusive)");
verbose = verbose_in;
Array<vec> means_in = model_in.get_means();
Array<vec> diag_covs_in = model_in.get_diag_covs();
vec weights_in = model_in.get_weights();
init(means_in, diag_covs_in, weights_in);
means_in.set_size(0);
diag_covs_in.set_size(0);
weights_in.set_size(0);
it_assert(check_size(X_in), "MOG_diag_kmeans_sup::run(): 'X' is empty or contains vectors of wrong dimensionality");
N = X_in.size();
if (K > N) {
it_warning("MOG_diag_kmeans_sup::run(): K > N");
}
else {
if (K > N / 10) {
it_warning("MOG_diag_kmeans_sup::run(): K > N/10");
}
}
max_iter = max_iter_in;
trust = trust_in;
means_old.set_size(K);
for (int k = 0;k < K;k++) means_old(k).set_size(D);
partitions.set_size(K);
for (int k = 0;k < K;k++) partitions(k).set_size(N);
count.set_size(K);
tmpvec.set_size(D);
norm_mu.set_size(D);
norm_sd.set_size(D);
c_X = enable_c_access(X_in);
c_means_old = enable_c_access(means_old);
c_partitions = enable_c_access(partitions);
c_count = enable_c_access(count);
c_tmpvec = enable_c_access(tmpvec);
c_norm_mu = enable_c_access(norm_mu);
c_norm_sd = enable_c_access(norm_sd);
if (normalise_in) normalise_vectors();
calc_means();
if (normalise_in) { unnormalise_vectors(); unnormalise_means(); }
calc_covs();
calc_weights();
model_in.init(means, diag_covs, weights);
disable_c_access(c_X);
disable_c_access(c_means_old);
disable_c_access(c_partitions);
disable_c_access(c_count);
disable_c_access(c_tmpvec);
disable_c_access(c_norm_mu);
disable_c_access(c_norm_sd);
means_old.set_size(0);
partitions.set_size(0);
count.set_size(0);
tmpvec.set_size(0);
norm_mu.set_size(0);
norm_sd.set_size(0);
cleanup();
}
//
// convenience functions
void MOG_diag_kmeans(MOG_diag &model_in, Array<vec> &X_in, int max_iter_in, double trust_in, bool normalise_in, bool verbose_in)
{
MOG_diag_kmeans_sup km;
km.run(model_in, X_in, max_iter_in, trust_in, normalise_in, verbose_in);
}
}
itpp::Array
General array class.
Definition array.h:105
itpp::Array::size
int size() const
Returns the number of data elements in the array object.
Definition array.h:155
itpp::Array::set_size
void set_size(int n, bool copy=false)
Resizing an Array<T>.
Definition array.h:257
itpp::MOG_diag_kmeans_sup
support class for MOG_diag_kmeans()
Definition mog_diag_kmeans.h:45
itpp::MOG_diag_kmeans_sup::measure_change
double measure_change() const
ADD DOCUMENTATION HERE.
Definition mog_diag_kmeans.cpp:130
itpp::MOG_diag_kmeans_sup::run
void run(MOG_diag &model_in, Array< vec > &X_in, int max_iter_in=10, double trust_in=0.5, bool normalise_in=true, bool verbose_in=false)
ADD DOCUMENTATION HERE.
Definition mog_diag_kmeans.cpp:264
itpp::MOG_diag_kmeans_sup::unnormalise_vectors
void unnormalise_vectors()
ADD DOCUMENTATION HERE.
Definition mog_diag_kmeans.cpp:244
itpp::MOG_diag_kmeans_sup::calc_means
void calc_means()
ADD DOCUMENTATION HERE.
Definition mog_diag_kmeans.cpp:183
itpp::MOG_diag_kmeans_sup::dezombify_means
bool dezombify_means()
ADD DOCUMENTATION HERE.
Definition mog_diag_kmeans.cpp:88
itpp::MOG_diag_kmeans_sup::c_count
int * c_count
'C' pointer to the count vector
Definition mog_diag_kmeans.h:119
itpp::MOG_diag_kmeans_sup::N
int N
number of training vectors
Definition mog_diag_kmeans.h:98
itpp::MOG_diag_kmeans_sup::c_X
double ** c_X
'C' pointers to training vectors
Definition mog_diag_kmeans.h:101
itpp::MOG_diag_kmeans_sup::initial_means
void initial_means()
ADD DOCUMENTATION HERE.
Definition mog_diag_kmeans.cpp:139
itpp::MOG_diag_kmeans_sup::normalise_vectors
void normalise_vectors()
ADD DOCUMENTATION HERE.
Definition mog_diag_kmeans.cpp:222
itpp::MOG_diag_kmeans_sup::means_old
Array< vec > means_old
means from the previous iteration, used to measure progress
Definition mog_diag_kmeans.h:104
itpp::MOG_diag_kmeans_sup::calc_covs
void calc_covs()
ADD DOCUMENTATION HERE.
Definition mog_diag_kmeans.cpp:190
itpp::MOG_diag_kmeans_sup::c_partitions
int ** c_partitions
'C' pointers to partition vectors
Definition mog_diag_kmeans.h:113
itpp::MOG_diag_kmeans_sup::assign_to_means
void assign_to_means()
ADD DOCUMENTATION HERE.
Definition mog_diag_kmeans.cpp:43
itpp::MOG_diag_kmeans_sup::unnormalise_means
void unnormalise_means()
ADD DOCUMENTATION HERE.
Definition mog_diag_kmeans.cpp:254
itpp::MOG_diag_kmeans_sup::c_means_old
double ** c_means_old
'C' pointers to old means
Definition mog_diag_kmeans.h:107
itpp::MOG_diag_kmeans_sup::dist
double dist(const double *x, const double *y) const
ADD DOCUMENTATION HERE.
Definition mog_diag_kmeans.cpp:36
itpp::MOG_diag_kmeans_sup::iterate
void iterate()
ADD DOCUMENTATION HERE.
Definition mog_diag_kmeans.cpp:161
itpp::MOG_diag_kmeans_sup::max_iter
int max_iter
Maximum number of iterations.
Definition mog_diag_kmeans.h:87
itpp::MOG_diag_kmeans_sup::trust
double trust
trust factor, where 0 <= trust <= 1.
Definition mog_diag_kmeans.h:92
itpp::MOG_diag_kmeans_sup::count
ivec count
keeps a count of the number of vectors assigned to each mean
Definition mog_diag_kmeans.h:116
itpp::MOG_diag_kmeans_sup::partitions
Array< ivec > partitions
contains indices of vectors assigned to each mean
Definition mog_diag_kmeans.h:110
itpp::MOG_diag_kmeans_sup::recalculate_means
void recalculate_means()
ADD DOCUMENTATION HERE.
Definition mog_diag_kmeans.cpp:65
itpp::MOG_diag_kmeans_sup::verbose
bool verbose
Whether we print the progress.
Definition mog_diag_kmeans.h:95
itpp::MOG_diag_kmeans_sup::calc_weights
void calc_weights()
ADD DOCUMENTATION HERE.
Definition mog_diag_kmeans.cpp:216
itpp::MOG_diag
Diagonal Mixture of Gaussians (MOG) class.
Definition mog_diag.h:56
itpp::MOG_diag::c_means
double ** c_means
pointers to the mean vectors
Definition mog_diag.h:197
itpp::MOG_diag::disable_c_access
double ** disable_c_access(double **A_in)
Disable C style access to an Array of vectors (vec)
Definition mog_diag.cpp:300
itpp::MOG_diag::enable_c_access
double ** enable_c_access(Array< vec > &A_in)
Enable C style access to an Array of vectors (vec)
Definition mog_diag.cpp:284
itpp::MOG_diag::cleanup
void cleanup()
Release memory used by the model. The model will be empty.
Definition mog_diag.h:111
itpp::MOG_diag::c_weights
double * c_weights
pointer to the weight vector
Definition mog_diag.h:206
itpp::MOG_diag::c_diag_covs
double ** c_diag_covs
pointers to the covariance vectors
Definition mog_diag.h:200
itpp::MOG_generic::check_size
bool check_size(const vec &x_in) const
Check if vector x_in has the same dimensionality as the model.
Definition mog_generic.cpp:130
itpp::MOG_generic::init
void init()
Initialise the model to be empty.
Definition mog_generic.cpp:43
itpp::MOG_generic::K
int K
number of gaussians
Definition mog_generic.h:292
itpp::MOG_generic::weights
vec weights
weights
Definition mog_generic.h:307
itpp::MOG_generic::diag_covs
Array< vec > diag_covs
diagonal covariance matrices, stored as vectors
Definition mog_generic.h:301
itpp::MOG_generic::D
int D
dimensionality
Definition mog_generic.h:295
itpp::MOG_generic::means
Array< vec > means
means
Definition mog_generic.h:298
itpp::MOG_diag_kmeans
void MOG_diag_kmeans(MOG_diag &model_in, Array< vec > &X_in, int max_iter_in, double trust_in, bool normalise_in, bool verbose_in)
Definition mog_diag_kmeans.cpp:347
it_warning
#define it_warning(s)
Display a warning message.
Definition itassert.h:173
it_assert
#define it_assert(t, s)
Abort if t is not true.
Definition itassert.h:94
mog_diag_kmeans.h
K-means based optimiser for Mixture of Gaussians - header file.
itpp
itpp namespace
Definition itmex.h:37
itpp::floor
vec floor(const vec &x)
Round to nearest lower integer.
Definition converters.h:346

Generated on Mon Apr 22 2024 09:45:27 for IT++ by Doxygen 1.9.8