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");
int size() const
Returns the number of data elements in the array object.
void set_size(int n, bool copy=false)
Resizing an Array<T>.
support class for MOG_diag_ML() and MOG_diag_MAP()
int max_iter
Maximum number of iterations.
double var_floor
ADD DOCUMENTATION HERE.
void sanitise_params()
ADD DOCUMENTATION HERE.
double weight_floor
ADD DOCUMENTATION HERE.
int N
number of training vectors
void update_internals()
ADD DOCUMENTATION HERE.
bool verbose
Whether we print the progress.
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.
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 ml_update_params()
ADD DOCUMENTATION HERE.
void ml_iterate()
ADD DOCUMENTATION HERE.
double ** c_X
'C' pointers to training vectors
Diagonal Mixture of Gaussians (MOG) class.
double ** c_diag_covs_inv_etc
pointers to the inverted covariance vectors
double log_lhood_single_gaus_internal(const double *c_x_in, const int k) const
ADD DOCUMENTATION HERE.
double ** c_means
pointers to the mean vectors
double ** disable_c_access(double **A_in)
Disable C style access to an Array of vectors (vec)
double ** enable_c_access(Array< vec > &A_in)
Enable C style access to an Array of vectors (vec)
void cleanup()
Release memory used by the model. The model will be empty.
double * c_log_det_etc
pointer to the log_det_etc vector
double * c_log_weights
pointer to the log version of the weight vector
double * c_weights
pointer to the weight vector
double ** c_diag_covs
pointers to the covariance vectors
bool check_array_uniformity(const Array< vec > &A) const
Check if all vectors in Array X_in have the same dimensionality.
void init()
Initialise the model to be empty.
vec get_weights() const
Obtain a copy of the weight vector.
Array< vec > diag_covs
diagonal covariance matrices, stored as vectors
Array< vec > get_diag_covs() const
Obtain a copy of the array of diagonal covariance vectors.
bool is_valid() const
Returns true if the model's parameters are valid.
Array< vec > get_means() const
Obtain a copy of the array of mean vectors.
double log_max_K
Pre-calcualted std::log(std::numeric_limits<double>::max() / K), where K is the number of Gaussians.
bool paranoid
indicates whether we are paranoid about numerical stability
double toc(void)
Returns the elapsed time since last tic()
void tic(void)
Resets the timer and starts it.
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)
#define it_error(s)
Abort unconditionally.
#define it_warning(s)
Display a warning message.
#define it_assert(t, s)
Abort if t is not true.
it_file & flush(it_file &f)
Flush operator.
double log_add(double log_a, double log_b)
Safe substitute for log(exp(log_a) + exp(log_b))
vec log(const vec &x)
The natural logarithm of the elements.
double trunc_exp(double x)
Truncated exponential function.
vec exp(const vec &x)
Exp of the elements of a vector x.
T sum(const Vec< T > &v)
Sum of all elements in the vector.
T min(const Vec< T > &in)
Minimum value of vector.
T max(const Vec< T > &v)
Maximum value of vector.
Logarithmic and exponenential functions - header file.
Expectation Maximisation (EM) based optimisers for MOG - header file.
const double m_2pi
Constant 2*Pi.
void MOG_diag_MAP(MOG_diag &, MOG_diag &, Array< vec > &, int, double, double, double, bool)
Definitions of Timing classes.