48 for (
int k = 0;k <
K;k++) {
53 for (
int d = 0;d <
D;d++) {
54 double tmp = c_diag_cov[d];
55 c_diag_cov_inv_etc[d] = 1.0 / (2.0 * tmp);
70 for (
int k = 0;k <
K;k++) {
77 for (
int k = 0;k <
K;k++)
78 for (
int d = 0;d <
D;d++)
87 double acc_loglhood = 0.0;
89 for (
int k = 0;k <
K;k++) {
90 c_acc_loglhood_K[k] = 0.0;
92 double * c_acc_mean = c_acc_means[k];
93 double * c_acc_cov = c_acc_covs[k];
95 for (
int d = 0;d <
D;d++) { c_acc_mean[d] = 0.0; c_acc_cov[d] = 0.0; }
98 for (
int n = 0;n <
N;n++) {
99 double * c_x =
c_X[n];
102 for (
int k = 0;k <
K;k++) {
110 double log_sum = c_tmpvecK[0];
111 for (
int k = 1;k <
K;k++) log_sum =
log_add(log_sum, c_tmpvecK[k]);
112 acc_loglhood += log_sum;
114 for (
int k = 0;k <
K;k++) {
116 double * c_acc_mean = c_acc_means[k];
117 double * c_acc_cov = c_acc_covs[k];
119 double tmp_k =
trunc_exp(c_tmpvecK[k] - log_sum);
120 acc_loglhood_K[k] += tmp_k;
122 for (
int d = 0;d <
D;d++) {
123 double tmp_x = c_x[d];
124 c_acc_mean[d] += tmp_k * tmp_x;
125 c_acc_cov[d] += tmp_k * tmp_x * tmp_x;
132 for (
int k = 0;k <
K;k++) {
double tmp =
std::exp(c_tmpvecK[k]); c_tmpvecK[k] = tmp; sum += tmp; }
135 for (
int k = 0;k <
K;k++) {
137 double * c_acc_mean = c_acc_means[k];
138 double * c_acc_cov = c_acc_covs[k];
140 double tmp_k = c_tmpvecK[k] /
sum;
141 c_acc_loglhood_K[k] += tmp_k;
143 for (
int d = 0;d <
D;d++) {
144 double tmp_x = c_x[d];
145 c_acc_mean[d] += tmp_k * tmp_x;
146 c_acc_cov[d] += tmp_k * tmp_x * tmp_x;
152 for (
int k = 0;k <
K;k++) {
157 double * c_acc_mean = c_acc_means[k];
158 double * c_acc_cov = c_acc_covs[k];
160 double tmp_k = c_acc_loglhood_K[k];
164 for (
int d = 0;d <
D;d++) {
165 double tmp_mean = c_acc_mean[d] / tmp_k;
166 c_mean[d] = tmp_mean;
167 c_diag_cov[d] = c_acc_cov[d] / tmp_k - tmp_mean * tmp_mean;
171 return(acc_loglhood / N);
182 using std::noshowpos;
183 using std::scientific;
186 using std::setprecision;
193 cout <<
"MOG_diag_EM_sup::ml_iterate()" << endl;
194 cout << setw(14) <<
"iteration";
195 cout << setw(14) <<
"avg_loglhood";
196 cout << setw(14) <<
"delta";
197 cout << setw(10) <<
"toc";
201 for (
int i = 0; i <
max_iter; i++) {
209 double delta = avg_log_lhood_new - avg_log_lhood_old;
211 cout << noshowpos << fixed;
212 cout << setw(14) << i;
213 cout << showpos << scientific << setprecision(3);
214 cout << setw(14) << avg_log_lhood_new;
215 cout << setw(14) << delta;
216 cout << noshowpos << fixed;
217 cout << setw(10) << tt.
toc();
218 cout << endl <<
flush;
221 if (avg_log_lhood_new <= avg_log_lhood_old)
break;
223 avg_log_lhood_old = avg_log_lhood_new;
231 it_assert(model_in.
is_valid(),
"MOG_diag_EM_sup::ml(): initial model not valid");
233 it_assert((max_iter_in > 0),
"MOG_diag_EM_sup::ml(): 'max_iter' needs to be greater than zero");
243 init(means_in, diag_covs_in, weights_in);
247 weights_in.set_size(0);
250 it_warning(
"MOG_diag_EM_sup::ml(): WARNING: K > N");
254 it_warning(
"MOG_diag_EM_sup::ml(): WARNING: K > N/10");
270 acc_loglhood_K.set_size(
K);
273 for (
int k = 0;k <
K;k++) acc_means(k).
set_size(
D);
275 for (
int k = 0;k <
K;k++) acc_covs(k).
set_size(
D);
298 acc_loglhood_K.set_size(0);
307 double,
double,
bool)
309 it_error(
"MOG_diag_EM_sup::map(): not implemented yet");
319 EM.
ml(model_in, X_in, max_iter_in, var_floor_in, weight_floor_in, verbose_in);
325 it_error(
"MOG_diag_MAP(): not implemented yet");
const double m_2pi
Constant 2*Pi.
double log_add(double log_a, double log_b)
Safe substitute for log(exp(log_a) + exp(log_b))
double ml_update_params()
ADD DOCUMENTATION HERE.
int size() const
Returns the number of data elements in the array object.
Array< vec > get_means() const
Obtain a copy of the array of mean vectors.
Array< vec > get_diag_covs() const
Obtain a copy of the array of diagonal covariance vectors.
Expectation Maximisation (EM) based optimisers for MOG - header file.
double * c_log_weights
pointer to the log version of the weight vector
double * c_log_det_etc
pointer to the log_det_etc vector
double toc(void)
Returns the elapsed time since last tic()
double ** c_diag_covs
pointers to the covariance vectors
void MOG_diag_MAP(MOG_diag &, MOG_diag &, Array< vec > &, int, double, double, double, bool)
T sum(const Vec< T > &v)
Sum of all elements in the vector.
bool is_valid() const
Returns true if the model's parameters are valid.
#define it_assert(t, s)
Abort if t is not true.
Logarithmic and exponenential functions - header file.
void set_size(int n, bool copy=false)
Resizing an Array<T>.
double ** c_means
pointers to the mean vectors
vec log(const vec &x)
The natural logarithm of the elements.
T min(const Vec< T > &in)
Minimum value of vector.
Diagonal Mixture of Gaussians (MOG) class.
vec exp(const vec &x)
Exp of the elements of a vector x.
T max(const Vec< T > &v)
Maximum value of vector.
double ** disable_c_access(double **A_in)
Disable C style access to an Array of vectors (vec)
double log_lhood_single_gaus_internal(const double *c_x_in, const int k) const
ADD DOCUMENTATION HERE.
void MOG_diag_ML(MOG_diag &model_in, Array< vec > &X_in, int max_iter_in, double var_floor_in, double weight_floor_in, bool verbose_in)
double var_floor
ADD DOCUMENTATION HERE.
it_file & flush(it_file &f)
Flush operatorFlushes the data. Usage:
void cleanup()
Release memory used by the model. The model will be empty.
void map(MOG_diag &model_in, MOG_diag &prior_model, Array< vec > &X_in, int max_iter_in=10, double alpha_in=0.5, double var_floor_in=0.0, double weight_floor_in=0.0, bool verbose_in=false)
ADD DOCUMENTATION HERE.
A real time timer classMeasures real time.
double ** c_diag_covs_inv_etc
pointers to the inverted covariance vectors
double ** enable_c_access(Array< vec > &A_in)
Enable C style access to an Array of vectors (vec)
void sanitise_params()
ADD DOCUMENTATION HERE.
bool verbose
Whether we print the progress.
void update_internals()
ADD DOCUMENTATION HERE.
void init()
Initialise the model to be empty.
double weight_floor
ADD DOCUMENTATION HERE.
#define it_warning(s)
Display a warning message.
int max_iter
Maximum number of iterations.
double trunc_exp(double x)
Truncated exponential function.
void ml_iterate()
ADD DOCUMENTATION HERE.
bool paranoid
indicates whether we are paranoid about numerical stability
vec get_weights() const
Obtain a copy of the weight vector.
#define it_error(s)
Abort unconditionally.
support class for MOG_diag_ML() and MOG_diag_MAP()
int N
number of training vectors
Array< vec > diag_covs
diagonal covariance matrices, stored as vectors
double log_max_K
Pre-calcualted std::log(std::numeric_limits<double>::max() / K), where K is the number of Gaussians...
void ml(MOG_diag &model_in, Array< vec > &X_in, int max_iter_in=10, double var_floor_in=0.0, double weight_floor_in=0.0, bool verbose_in=false)
ADD DOCUMENTATION HERE.
double * c_weights
pointer to the weight vector
bool check_array_uniformity(const Array< vec > &A) const
Check if all vectors in Array X_in have the same dimensionality.
double ** c_X
'C' pointers to training vectors
void tic(void)
Resets the timer and starts it.
Definitions of Timing classes.