faster folders removed
This commit is contained in:
Parikshit Ram
@@ -1,26 +0,0 @@
|
||||
|
||||
librule(
|
||||
name = "mog_em", # this line can be safely omitted
|
||||
sources = ["mog.cc"], # files that must be compiled
|
||||
headers = ["mog.h","phi.h","math_functions.h"],# include files part of the 'lib'
|
||||
deplibs = ["fastlib:fastlib"], # depends on fastlib core
|
||||
#tests = ["mog_em_tests.cc"] # this file contains a main with test functions
|
||||
)
|
||||
|
||||
binrule(
|
||||
name = "mog_em_main", # the executable name
|
||||
sources = ["mog_em_main.cc"], # compile main.cc
|
||||
#headers = ["mog_em_main.h"], # no extra headers
|
||||
deplibs = [":mog_em", "fastlib:fastlib"] #
|
||||
)
|
||||
|
||||
# to build:
|
||||
# 1. make sure have environment variables set up:
|
||||
# $ source /full/path/to/fastlib/script/fl-env /full/path/to/fastlib
|
||||
# (you might want to put this in bashrc)
|
||||
# 2. fl-build main
|
||||
# - this automatically will assume --mode=check, the default
|
||||
# - type fl-build --help for help
|
||||
# 3. ./main
|
||||
# - to build same target again, type: make
|
||||
# - to force recompilation, type: make clean
|
||||
File diff suppressed because it is too large
Load Diff
@@ -1,138 +0,0 @@
|
||||
/**
|
||||
* @author Parikshit Ram (pram@cc.gatech.edu)
|
||||
* @file math_functions.h
|
||||
*
|
||||
* This file has certain functions that find the
|
||||
* highest or lowest element in an array or
|
||||
* in a row of a matrix and returns them
|
||||
*
|
||||
*/
|
||||
|
||||
#include "fastlib/fastlib.h"
|
||||
#include "fastlib/fastlib_int.h"
|
||||
|
||||
|
||||
/**
|
||||
* Finds the index of the minimum element
|
||||
* in each row of a matrix
|
||||
*
|
||||
* Example use:
|
||||
* @code
|
||||
* Matrix& mat;
|
||||
* index_t indices[mat.n_rows()];
|
||||
*
|
||||
* ...
|
||||
*
|
||||
* min_element(mat, indices);
|
||||
* @endcode
|
||||
*/
|
||||
|
||||
void min_element(Matrix& element, index_t *indices) {
|
||||
|
||||
index_t last = element.n_cols() - 1;
|
||||
index_t first, lowest;
|
||||
index_t i;
|
||||
|
||||
for (i = 0; i < element.n_rows(); i++) {
|
||||
|
||||
first = lowest = 0;
|
||||
if (first == last) {
|
||||
indices[i] = last;
|
||||
}
|
||||
while (++first <= last) {
|
||||
if (element.get(i, first) < element.get(i, lowest)) {
|
||||
lowest = first;
|
||||
}
|
||||
}
|
||||
indices[i] = lowest;
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the index of the maximum element
|
||||
* in a float array 'array' of length 'length'.
|
||||
*
|
||||
* Example use:
|
||||
* @code
|
||||
* index_t length, index;
|
||||
* float array[length];
|
||||
* ...
|
||||
* index = max_element_index(array, length);
|
||||
* @endcode
|
||||
*/
|
||||
int max_element_index(float *array, int length) {
|
||||
|
||||
int last = length - 1;
|
||||
int first = 0;
|
||||
int highest = 0;
|
||||
|
||||
if (first == last) {
|
||||
return last;
|
||||
}
|
||||
while (++first <= last) {
|
||||
if (array[first] > array[highest]) {
|
||||
highest = first;
|
||||
}
|
||||
}
|
||||
return highest;
|
||||
}
|
||||
|
||||
/**
|
||||
* Finds the index of the maximum element in an arraylist
|
||||
* of floats
|
||||
*
|
||||
* Example use:
|
||||
* @code
|
||||
* index_t index;
|
||||
* ArrayList<float> array;
|
||||
* ...
|
||||
* index = max_element_index(array);
|
||||
* @endcode
|
||||
*/
|
||||
int max_element_index(ArrayList<float>& array){
|
||||
|
||||
int last = array.size() - 1;
|
||||
int first = 0;
|
||||
int highest = 0;
|
||||
|
||||
if (first == last) {
|
||||
return last;
|
||||
}
|
||||
while (++first <= last) {
|
||||
if (array[first] > array[highest]) {
|
||||
highest = first;
|
||||
}
|
||||
}
|
||||
return highest;
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the index of the maximum element in
|
||||
* an arraylist of doubles
|
||||
*
|
||||
* Example use:
|
||||
* @code
|
||||
* index_t index;
|
||||
* ArrayList<double> array;
|
||||
* ...
|
||||
* index = max_element_index(array);
|
||||
* @endcode
|
||||
*/
|
||||
int max_element_index(ArrayList<double>& array) {
|
||||
|
||||
int last = array.size() - 1;
|
||||
int first = 0;
|
||||
int highest = 0;
|
||||
|
||||
if (first == last) {
|
||||
return last;
|
||||
}
|
||||
while (++first <= last) {
|
||||
if(array[first] > array[highest]) {
|
||||
highest = first;
|
||||
}
|
||||
}
|
||||
return highest;
|
||||
}
|
||||
|
||||
@@ -1,301 +0,0 @@
|
||||
/**
|
||||
* @author Parikshit Ram (pram@cc.gatech.edu)
|
||||
* @file mog.cc
|
||||
*
|
||||
* Implementation for the loglikelihood function, the EM algorithm
|
||||
* and also computes the K-means for getting an initial point
|
||||
*
|
||||
*/
|
||||
|
||||
#include "mog.h"
|
||||
#include "phi.h"
|
||||
#include "math_functions.h"
|
||||
|
||||
void MoGEM::ExpectationMaximization(Matrix& data_points) {
|
||||
|
||||
// Declaration of the variables */
|
||||
index_t num_points;
|
||||
index_t dim, num_gauss;
|
||||
double sum, tmp;
|
||||
ArrayList<Vector> mu_temp, mu;
|
||||
ArrayList<Matrix> sigma_temp, sigma;
|
||||
Vector omega_temp, omega, x;
|
||||
Matrix cond_prob;
|
||||
long double l, l_old, best_l, INFTY = 99999, TINY = 1.0e-10;
|
||||
|
||||
// Initializing values
|
||||
dim = dimension();
|
||||
num_gauss = number_of_gaussians();
|
||||
num_points = data_points.n_cols();
|
||||
|
||||
// Initializing the number of the vectors and matrices
|
||||
// according to the parameters input
|
||||
mu_temp.Init(num_gauss);
|
||||
mu.Init(num_gauss);
|
||||
sigma_temp.Init(num_gauss);
|
||||
sigma.Init(num_gauss);
|
||||
omega_temp.Init(num_gauss);
|
||||
omega.Init(num_gauss);
|
||||
|
||||
// Allocating size to the vectors and matrices
|
||||
// according to the dimensionality of the data
|
||||
for(index_t i = 0; i < num_gauss; i++) {
|
||||
mu_temp[i].Init(dim);
|
||||
mu[i].Init(dim);
|
||||
sigma_temp[i].Init(dim, dim);
|
||||
sigma[i].Init(dim, dim);
|
||||
}
|
||||
x.Init(dim);
|
||||
cond_prob.Init(num_gauss, num_points);
|
||||
|
||||
best_l = -INFTY;
|
||||
index_t restarts = 0;
|
||||
// performing 5 restarts and choosing the best from them
|
||||
while (restarts < 1) {
|
||||
|
||||
// assign initial values to 'mu', 'sig' and 'omega' using k-means
|
||||
KMeans(data_points, &mu_temp, &sigma_temp, &omega_temp, num_gauss);
|
||||
|
||||
l_old = -INFTY;
|
||||
|
||||
// calculates the loglikelihood value
|
||||
l = Loglikelihood(data_points, mu_temp, sigma_temp, omega_temp);
|
||||
|
||||
// added a check here to see if any
|
||||
// significant change is being made
|
||||
// at every iteration
|
||||
while (l - l_old > TINY) {
|
||||
// calculating the conditional probabilities
|
||||
// of choosing a particular gaussian given
|
||||
// the data and the present theta value
|
||||
/**********
|
||||
for (index_t j = 0; j < num_points; j++) {
|
||||
x.CopyValues(data_points.GetColumnPtr(j));
|
||||
sum = 0;
|
||||
for (index_t i = 0; i < num_gauss; i++) {
|
||||
tmp = phi(x, mu_temp[i], sigma_temp[i]) * omega_temp.get(i); // can be made faster by sending all the data at once instead of looping over
|
||||
cond_prob.set(i, j, tmp);
|
||||
sum += tmp;
|
||||
}
|
||||
for (index_t i = 0; i < num_gauss; i++) {
|
||||
tmp = cond_prob.get(i, j);
|
||||
cond_prob.set(i, j, tmp / sum);
|
||||
}
|
||||
}
|
||||
******/
|
||||
|
||||
/*******FIX THIS SOON********/
|
||||
Vector phi_val;
|
||||
phi_val.Init(num_points);
|
||||
|
||||
for(index_t i = 0; i < num_gauss; i++) {
|
||||
phi_val.CopyValues(phi(data_points, mu_temp[i], sigma_temp[i]));
|
||||
la::Scale(omega_temp.get(i), &phi_val);
|
||||
*/
|
||||
|
||||
// calculating the new value of the mu
|
||||
// using the updated conditional probabilities
|
||||
for (index_t i = 0; i < num_gauss; i++) {
|
||||
sum = 0;
|
||||
mu_temp[i].SetZero();
|
||||
for (index_t j = 0; j < num_points; j++) {
|
||||
x.CopyValues(data_points.GetColumnPtr(j));
|
||||
la::AddExpert(cond_prob.get(i, j), x, &mu_temp[i]);
|
||||
sum += cond_prob.get(i, j);
|
||||
}
|
||||
la::Scale((1.0 / sum), &mu_temp[i]);
|
||||
}
|
||||
|
||||
// calculating the new value of the sig
|
||||
// using the updated conditional probabilities
|
||||
// and the updated mu
|
||||
for (index_t i = 0; i < num_gauss; i++) {
|
||||
sum = 0;
|
||||
sigma_temp[i].SetZero();
|
||||
for (index_t j = 0; j < num_points; j++) {
|
||||
Matrix co, ro, c;
|
||||
c.Init(dim, dim);
|
||||
x.CopyValues(data_points.GetColumnPtr(j));
|
||||
la::SubFrom(mu_temp[i] , &x);
|
||||
co.AliasColVector(x);
|
||||
ro.AliasRowVector(x);
|
||||
la::MulOverwrite(co, ro, &c);
|
||||
la::AddExpert(cond_prob.get(i, j), c, &sigma_temp[i]);
|
||||
sum += cond_prob.get(i, j);
|
||||
}
|
||||
la::Scale((1.0 / sum), &sigma_temp[i]);
|
||||
}
|
||||
|
||||
// calculating the new values for omega
|
||||
// using the updated conditional probabilities
|
||||
Vector identity_vector;
|
||||
identity_vector.Init(num_points);
|
||||
identity_vector.SetAll(1.0 / num_points);
|
||||
la::MulOverwrite(cond_prob, identity_vector, &omega_temp);
|
||||
|
||||
l_old = l;
|
||||
l = Loglikelihood(data_points, mu_temp, sigma_temp, omega_temp);
|
||||
}
|
||||
|
||||
// putting a check to see if the best one is chosen
|
||||
if(l > best_l){
|
||||
best_l = l;
|
||||
for (index_t i = 0; i < num_gauss; i++) {
|
||||
mu[i].CopyValues(mu_temp[i]);
|
||||
sigma[i].CopyValues(sigma_temp[i]);
|
||||
}
|
||||
omega.CopyValues(omega_temp);
|
||||
}
|
||||
restarts++;
|
||||
}
|
||||
|
||||
for (index_t i = 0; i < num_gauss; i++) {
|
||||
set_mu(i, mu[i]);
|
||||
set_sigma(i, sigma[i]);
|
||||
}
|
||||
set_omega(omega);
|
||||
|
||||
NOTIFY("loglikelihood value of the estimated model: %Lf\n", best_l);
|
||||
return;
|
||||
}
|
||||
|
||||
long double MoGEM::Loglikelihood(Matrix& data_points, ArrayList<Vector>& means,
|
||||
ArrayList<Matrix>& covars, Vector& weights) {
|
||||
|
||||
index_t i, j;
|
||||
Vector x;
|
||||
long double likelihood, loglikelihood = 0;
|
||||
|
||||
x.Init(data_points.n_rows());
|
||||
|
||||
for (j = 0; j < data_points.n_cols(); j++) {
|
||||
x.CopyValues(data_points.GetColumnPtr(j));
|
||||
likelihood = 0;
|
||||
for(i = 0; i < number_of_gaussians() ; i++){
|
||||
likelihood += weights.get(i) * phi(x, means[i], covars[i]);
|
||||
}
|
||||
loglikelihood += log(likelihood);
|
||||
}
|
||||
return loglikelihood;
|
||||
}
|
||||
|
||||
void MoGEM::KMeans(Matrix& data, ArrayList<Vector> *means,
|
||||
ArrayList<Matrix> *covars, Vector *weights, index_t value_of_k){
|
||||
|
||||
|
||||
ArrayList<Vector> mu, mu_old;
|
||||
double* tmpssq;
|
||||
double* sig;
|
||||
double* sig_best;
|
||||
index_t *y;
|
||||
Vector x, diff;
|
||||
Matrix ssq;
|
||||
index_t i, j, k, n, t, dim;
|
||||
double score, score_old, sum;
|
||||
|
||||
n = data.n_cols();
|
||||
dim = data.n_rows();
|
||||
mu.Init(value_of_k);
|
||||
mu_old.Init(value_of_k);
|
||||
tmpssq = (double*)malloc(value_of_k * sizeof( double ));
|
||||
sig = (double*)malloc(value_of_k * sizeof( double ));
|
||||
sig_best = (double*)malloc(value_of_k * sizeof( double ));
|
||||
ssq.Init(n, value_of_k);
|
||||
|
||||
for( i = 0; i < value_of_k; i++){
|
||||
mu[i].Init(dim);
|
||||
mu_old[i].Init(dim);
|
||||
}
|
||||
x.Init(dim);
|
||||
y = (index_t*)malloc(n * sizeof(index_t));
|
||||
diff.Init(dim);
|
||||
|
||||
score_old = 999999;
|
||||
|
||||
// putting 5 random restarts to obtain the k-means
|
||||
for(i = 0; i < 5; i++){
|
||||
t = -1;
|
||||
for (k = 0; k < value_of_k; k++){
|
||||
t = (t + 1 + (rand()%((n - 1 - (value_of_k - k)) - (t + 1))));
|
||||
mu[k].CopyValues(data.GetColumnPtr(t));
|
||||
for(j = 0; j < n; j++){
|
||||
x.CopyValues( data.GetColumnPtr(j));
|
||||
la::SubOverwrite(mu[k], x, &diff);
|
||||
ssq.set( j, k, la::Dot(diff, diff));
|
||||
}
|
||||
}
|
||||
min_element(ssq, y);
|
||||
|
||||
do{
|
||||
for(k = 0; k < value_of_k; k++){
|
||||
mu_old[k].CopyValues(mu[k]);
|
||||
}
|
||||
|
||||
for(k = 0; k < value_of_k; k++){
|
||||
index_t p = 0;
|
||||
mu[k].SetZero();
|
||||
for(j = 0; j < n; j++){
|
||||
x.CopyValues(data.GetColumnPtr(j));
|
||||
if(y[j] == k){
|
||||
la::AddTo(x, &mu[k]);
|
||||
p++;
|
||||
}
|
||||
}
|
||||
|
||||
if(p == 0){
|
||||
}
|
||||
else{
|
||||
double sc = 1 ;
|
||||
sc = sc / p;
|
||||
la::Scale(sc , &mu[k]);
|
||||
}
|
||||
for(j = 0; j < n; j++){
|
||||
x.CopyValues(data.GetColumnPtr(j));
|
||||
la::SubOverwrite(mu[k], x, &diff);
|
||||
ssq.set(j, k, la::Dot(diff, diff));
|
||||
}
|
||||
}
|
||||
min_element(ssq, y);
|
||||
|
||||
sum = 0;
|
||||
for(k = 0; k < value_of_k; k++) {
|
||||
la::SubOverwrite(mu[k], mu_old[k], &diff);
|
||||
sum += la::Dot(diff, diff);
|
||||
}
|
||||
}while(sum != 0);
|
||||
|
||||
for(k = 0; k < value_of_k; k++){
|
||||
index_t p = 0;
|
||||
tmpssq[k] = 0;
|
||||
for(j = 0; j < n; j++){
|
||||
if(y[j] == k){
|
||||
tmpssq[k] += ssq.get(j, k);
|
||||
p++;
|
||||
}
|
||||
}
|
||||
sig[k] = sqrt(tmpssq[k] / p);
|
||||
}
|
||||
|
||||
score = 0;
|
||||
for(k = 0; k < value_of_k; k++){
|
||||
score += tmpssq[k];
|
||||
}
|
||||
score = score / n;
|
||||
|
||||
if (score < score_old) {
|
||||
score_old = score;
|
||||
for(k = 0; k < value_of_k; k++){
|
||||
(*means)[k].CopyValues(mu[k]);
|
||||
sig_best[k] = sig[k];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
for(k = 0; k < value_of_k; k++){
|
||||
x.SetAll(sig_best[k]);
|
||||
(*covars)[k].SetDiagonal(x);
|
||||
}
|
||||
double tmp = 1;
|
||||
(*weights).SetAll(tmp / value_of_k);
|
||||
return;
|
||||
}
|
||||
@@ -1,244 +0,0 @@
|
||||
/**
|
||||
* @author Parikshit Ram (pram@cc.gatech.edu)
|
||||
* @file mog.h
|
||||
*
|
||||
* Defines a Gaussian Mixture model and
|
||||
* estimates the parameters of the model
|
||||
*/
|
||||
|
||||
#ifndef MOGEM_H
|
||||
#define MOGEM_H
|
||||
|
||||
#include <fastlib/fastlib.h>
|
||||
|
||||
/**
|
||||
* A Gaussian mixture model class.
|
||||
*
|
||||
* This class uses maximum likelihood loss functions to
|
||||
* estimate the parameters of a gaussian mixture
|
||||
* model on a given data via the EM algorithm.
|
||||
*
|
||||
*
|
||||
* Example use:
|
||||
*
|
||||
* @code
|
||||
* MoGEM mog;
|
||||
* ArrayList<double> results;
|
||||
*
|
||||
* mog.Init(number_of_gaussians, dimension);
|
||||
* mog.ExpectationMaximization(data, &results, optim_flag);
|
||||
* @endcode
|
||||
*/
|
||||
class MoGEM {
|
||||
|
||||
private:
|
||||
|
||||
// The parameters of the mixture model
|
||||
ArrayList<Vector> mu_;
|
||||
ArrayList<Matrix> sigma_;
|
||||
Vector omega_;
|
||||
index_t number_of_gaussians_;
|
||||
index_t dimension_;
|
||||
|
||||
public:
|
||||
|
||||
MoGEM() {
|
||||
mu_.Init(0);
|
||||
sigma_.Init(0);
|
||||
}
|
||||
|
||||
~MoGEM() {
|
||||
}
|
||||
|
||||
void Init(index_t num_gauss, index_t dimension) {
|
||||
|
||||
// Initialize the private variables
|
||||
number_of_gaussians_ = num_gauss;
|
||||
dimension_ = dimension;
|
||||
|
||||
// Resize the ArrayList of Vectors and Matrices
|
||||
mu_.Resize(number_of_gaussians_);
|
||||
sigma_.Resize(number_of_gaussians_);
|
||||
}
|
||||
|
||||
void Init(datanode *mog_em_module) {
|
||||
|
||||
index_t num_gauss = fx_param_int_req(mog_em_module, "K");
|
||||
index_t dim = fx_param_int_req(mog_em_module, "D");
|
||||
Init(num_gauss, dim);
|
||||
}
|
||||
// The get functions
|
||||
|
||||
ArrayList<Vector>& mu() {
|
||||
return mu_;
|
||||
}
|
||||
|
||||
ArrayList<Matrix>& sigma() {
|
||||
return sigma_;
|
||||
}
|
||||
|
||||
Vector& omega() {
|
||||
return omega_;
|
||||
}
|
||||
|
||||
index_t number_of_gaussians() {
|
||||
return number_of_gaussians_;
|
||||
}
|
||||
|
||||
index_t dimension() {
|
||||
return dimension_;
|
||||
}
|
||||
|
||||
Vector& mu(index_t i) {
|
||||
return mu_[i] ;
|
||||
}
|
||||
|
||||
Matrix& sigma(index_t i) {
|
||||
return sigma_[i];
|
||||
}
|
||||
|
||||
double omega(index_t i) {
|
||||
return omega_.get(i);
|
||||
}
|
||||
|
||||
// The set functions
|
||||
|
||||
void set_mu(index_t i, Vector& mu) {
|
||||
DEBUG_ASSERT(i < number_of_gaussians());
|
||||
DEBUG_ASSERT(mu.length() == dimension());
|
||||
mu_[i].Copy(mu);
|
||||
return;
|
||||
}
|
||||
|
||||
void set_mu(index_t i, index_t length, const double *mu) {
|
||||
DEBUG_ASSERT(i < number_of_gaussians());
|
||||
DEBUG_ASSERT(length == dimension());
|
||||
mu_[i].Copy(mu, length);
|
||||
return;
|
||||
}
|
||||
|
||||
void set_sigma(index_t i, Matrix& sigma) {
|
||||
DEBUG_ASSERT(i < number_of_gaussians());
|
||||
DEBUG_ASSERT(sigma.n_rows() == dimension());
|
||||
DEBUG_ASSERT(sigma.n_cols() == dimension());
|
||||
sigma_[i].Copy(sigma);
|
||||
return;
|
||||
}
|
||||
|
||||
void set_omega(Vector& omega) {
|
||||
DEBUG_ASSERT(omega.length() == number_of_gaussians());
|
||||
omega_.Copy(omega);
|
||||
return;
|
||||
}
|
||||
|
||||
void set_omega(index_t length, const double *omega) {
|
||||
DEBUG_ASSERT(length == number_of_gaussians());
|
||||
omega_.Copy(omega, length);
|
||||
return;
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* This function outputs the parameters of the model
|
||||
* to an arraylist of doubles
|
||||
*
|
||||
* @code
|
||||
* ArrayList<double> results;
|
||||
* mog.OutputResults(&results);
|
||||
* @endcode
|
||||
*/
|
||||
void OutputResults(ArrayList<double> *results) {
|
||||
|
||||
// Initialize the size of the output array
|
||||
(*results).Init(number_of_gaussians_ * (1 + dimension_*(1 + dimension_)));
|
||||
|
||||
// Copy values to the array from the private variables of the class
|
||||
for (index_t i = 0; i < number_of_gaussians_; i++) {
|
||||
(*results)[i] = omega(i);
|
||||
for (index_t j = 0; j < dimension_; j++) {
|
||||
(*results)[number_of_gaussians_ + i*dimension_ + j] = mu(i).get(j);
|
||||
for (index_t k = 0; k < dimension_; k++) {
|
||||
(*results)[number_of_gaussians_*(1 + dimension_)
|
||||
+ i*dimension_*dimension_ + j*dimension_
|
||||
+ k] = sigma(i).get(j, k);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* This function prints the parameters of the model
|
||||
*
|
||||
* @code
|
||||
* mog.Display();
|
||||
* @endcode
|
||||
*/
|
||||
void Display(){
|
||||
|
||||
// Output the model parameters as the omega, mu and sigma
|
||||
printf(" Omega : [ ");
|
||||
for (index_t i = 0; i < number_of_gaussians_; i++) {
|
||||
printf("%lf ", omega(i));
|
||||
}
|
||||
printf("]\n");
|
||||
printf(" Mu : \n[");
|
||||
for (index_t i = 0; i < number_of_gaussians_; i++) {
|
||||
for (index_t j = 0; j < dimension_ ; j++) {
|
||||
printf("%lf ", mu(i).get(j));
|
||||
}
|
||||
printf(";");
|
||||
if (i == (number_of_gaussians_ - 1)) {
|
||||
printf("\b]\n");
|
||||
}
|
||||
}
|
||||
printf("Sigma : ");
|
||||
for (index_t i = 0; i < number_of_gaussians_; i++) {
|
||||
printf("\n[");
|
||||
for (index_t j = 0; j < dimension_ ; j++) {
|
||||
for(index_t k = 0; k < dimension_ ; k++) {
|
||||
printf("%lf ",sigma(i).get(j, k));
|
||||
}
|
||||
printf(";");
|
||||
}
|
||||
printf("\b]");
|
||||
}
|
||||
printf("\n");
|
||||
}
|
||||
|
||||
/**
|
||||
* This function calculates the parameters of the model
|
||||
* using the Maximum Likelihood function via the
|
||||
* Expectation Maximization (EM) Algorithm.
|
||||
*
|
||||
* @code
|
||||
* MoG mog;
|
||||
* Matrix data = "the data on which you want to fit the model";
|
||||
* ArrayList<double> results;
|
||||
* mog.ExpectationMaximization(data, &results);
|
||||
* @endcode
|
||||
*/
|
||||
void ExpectationMaximization(Matrix& data_points);
|
||||
|
||||
/**
|
||||
* This function computes the loglikelihood of model.
|
||||
* This function is used by the 'ExpectationMaximization'
|
||||
* function.
|
||||
*
|
||||
*/
|
||||
long double Loglikelihood(Matrix& data_points, ArrayList<Vector>& means,
|
||||
ArrayList<Matrix>& covars, Vector& weights);
|
||||
|
||||
/**
|
||||
* This function computes the k-means of the data and stores
|
||||
* the calculated means and covariances in the ArrayList
|
||||
* of Vectors and Matrices passed to it. It sets the weights
|
||||
* uniformly.
|
||||
*
|
||||
* This function is used to obtain a starting point for
|
||||
* the optimization
|
||||
*/
|
||||
void KMeans(Matrix& data, ArrayList<Vector> *means,
|
||||
ArrayList<Matrix> *covars, Vector *weights, index_t value_of_k);
|
||||
};
|
||||
|
||||
#endif
|
||||
@@ -1,72 +0,0 @@
|
||||
/**
|
||||
* @author Parikshit Ram (pram@cc.gatech.edu)
|
||||
* @file mog_l2e_main.cc
|
||||
*
|
||||
* This program test drives the L2 estimation
|
||||
* of a Gaussian Mixture model.
|
||||
*
|
||||
* PARAMETERS TO BE INPUT:
|
||||
*
|
||||
* --data
|
||||
* This is the file that contains the data on which
|
||||
* the model is to be fit
|
||||
*
|
||||
* --mog_em/K
|
||||
* This is the number of gaussians we want to fit
|
||||
* on the data, defaults to '1'
|
||||
*
|
||||
* --output
|
||||
* This file will contain the parameters estimated,
|
||||
* defaults to 'ouotput.csv'
|
||||
*
|
||||
*/
|
||||
|
||||
#include "mog.h"
|
||||
|
||||
|
||||
int main(int argc, char* argv[]) {
|
||||
|
||||
fx_init(argc, argv);
|
||||
|
||||
////// READING PARAMETERS AND LOADING DATA //////
|
||||
|
||||
const char *data_filename = fx_param_str_req(NULL, "data");
|
||||
|
||||
Matrix data_points;
|
||||
data::Load(data_filename, &data_points);
|
||||
|
||||
////// MIXTURE OF GAUSSIANS USING EM //////
|
||||
|
||||
MoGEM mog;
|
||||
|
||||
struct datanode* mog_em_module = fx_submodule(NULL, "mog_em", "mog_em");
|
||||
fx_param_int(mog_em_module, "K", 1);
|
||||
fx_format_param(mog_em_module, "D", "%d", data_points.n_rows());
|
||||
|
||||
////// Timing the initialization of the mixture model //////
|
||||
fx_timer_start(mog_em_module, "model_init");
|
||||
mog.Init(mog_em_module);
|
||||
fx_timer_stop(mog_em_module, "model_init");
|
||||
|
||||
////// Computing the parameters of the model using the EM algorithm //////
|
||||
ArrayList<double> results;
|
||||
|
||||
fx_timer_start(mog_em_module, "EM");
|
||||
mog.ExpectationMaximization(data_points);
|
||||
fx_timer_stop(mog_em_module, "EM");
|
||||
|
||||
mog.Display();
|
||||
mog.OutputResults(&results);
|
||||
|
||||
////// OUTPUT RESULTS //////
|
||||
|
||||
const char *output_filename = fx_param_str(NULL, "output", "output.csv");
|
||||
|
||||
FILE *output_file = fopen(output_filename, "w");
|
||||
|
||||
ot::Print(results, output_file);
|
||||
fclose(output_file);
|
||||
fx_done();
|
||||
|
||||
return 1;
|
||||
}
|
||||
@@ -1,142 +0,0 @@
|
||||
/**
|
||||
* @author Parikshit Ram (pram@cc.gatech.edu)
|
||||
* @file phi.h
|
||||
*
|
||||
* This file computes the Gaussian probability
|
||||
* density function
|
||||
*/
|
||||
#include "fastlib/fastlib.h"
|
||||
#include "fastlib/fastlib_int.h"
|
||||
#include <cmath>
|
||||
|
||||
/**
|
||||
* Calculates the multivariate Gaussian probability density function
|
||||
*
|
||||
* Example use:
|
||||
* @code
|
||||
* Vector x, mean;
|
||||
* Matrix cov;
|
||||
* ....
|
||||
* long double f = phi(x, mean, cov);
|
||||
* @endcode
|
||||
*/
|
||||
|
||||
long double phi(Vector& x , Vector& mean , Matrix& cov) {
|
||||
|
||||
long double det, f;
|
||||
double exponent;
|
||||
index_t dim;
|
||||
Matrix inv;
|
||||
Vector diff, tmp;
|
||||
|
||||
dim = x.length();
|
||||
la::InverseInit(cov, &inv);
|
||||
det = la::Determinant(cov);
|
||||
|
||||
if( det < 0){
|
||||
det = -det;
|
||||
}
|
||||
la::SubInit(mean,x,&diff);
|
||||
la::MulInit(inv, diff, &tmp);
|
||||
exponent = la::Dot(diff, tmp);
|
||||
long double tmp1, tmp2, tmp3;
|
||||
tmp1 = 1;
|
||||
tmp2 = dim;
|
||||
tmp2 = tmp2/2;
|
||||
tmp2 = pow((2*(math::PI)),tmp2);
|
||||
tmp1 = tmp1/tmp2;
|
||||
tmp3 = 1;
|
||||
tmp2 = sqrt(det);
|
||||
tmp3 = tmp3/tmp2;
|
||||
tmp2 = -exponent;
|
||||
tmp2 = tmp2 / 2;
|
||||
f = (tmp1*tmp3*exp(tmp2));
|
||||
|
||||
return f;
|
||||
}
|
||||
|
||||
/**
|
||||
* Calculates the univariate Gaussian probability density function
|
||||
*
|
||||
* Example use:
|
||||
* @code
|
||||
* double x, mean, var;
|
||||
* ....
|
||||
* long double f = phi(x, mean, var);
|
||||
* @endcode
|
||||
*/
|
||||
|
||||
long double phi(double x, double mean, double var) {
|
||||
|
||||
long double f;
|
||||
|
||||
f = exp(-1.0*((x-mean)*(x-mean)/(2*var)))/sqrt(2*math::PI*var);
|
||||
return f;
|
||||
}
|
||||
|
||||
/**
|
||||
* Calculates the multivariate Gaussian probability density function
|
||||
* and also the gradients with respect to the mean and the variance
|
||||
*
|
||||
* Example use:
|
||||
* @code
|
||||
* Vector x, mean, g_mean, g_cov;
|
||||
* ArrayList<Matrix> d_cov; // the dSigma
|
||||
* ....
|
||||
* long double f = phi(x, mean, cov, d_cov, &g_mean, &g_cov);
|
||||
* @endcode
|
||||
*/
|
||||
|
||||
long double phi(Vector& x, Vector& mean, Matrix& cov, ArrayList<Matrix>& d_cov, Vector *g_mean, Vector *g_cov){
|
||||
|
||||
long double det, f;
|
||||
double exponent;
|
||||
index_t dim;
|
||||
Matrix inv;
|
||||
Vector diff, tmp;
|
||||
|
||||
dim = x.length();
|
||||
la::InverseInit(cov, &inv);
|
||||
det = la::Determinant(cov);
|
||||
|
||||
if( det < 0){
|
||||
det = -det;
|
||||
}
|
||||
la::SubInit(mean,x,&diff);
|
||||
la::MulInit(inv, diff, &tmp);
|
||||
exponent = la::Dot(diff, tmp);
|
||||
long double tmp1, tmp2, tmp3;
|
||||
tmp1 = 1;
|
||||
tmp2 = dim;
|
||||
tmp2 = tmp2/2;
|
||||
tmp2 = pow((2*(math::PI)),tmp2);
|
||||
tmp1 = tmp1/tmp2;
|
||||
tmp3 = 1;
|
||||
tmp2 = sqrt(det);
|
||||
tmp3 = tmp3/tmp2;
|
||||
tmp2 = -exponent;
|
||||
tmp2 = tmp2 / 2;
|
||||
f = (tmp1*tmp3*exp(tmp2));
|
||||
|
||||
// Calculating the g_mean values which would be a (1 X dim) vector
|
||||
la::ScaleInit(f,tmp,g_mean);
|
||||
|
||||
// Calculating the g_cov values which would be a (1 X (dim*(dim+1)/2)) vector
|
||||
double *g_cov_tmp;
|
||||
g_cov_tmp = (double*)malloc(d_cov.size()*sizeof(double));
|
||||
for(index_t i = 0; i < d_cov.size(); i++){
|
||||
Vector tmp_d;
|
||||
Matrix inv_d;
|
||||
long double tmp_d_cov_d_r;
|
||||
|
||||
la::MulInit(d_cov[i],tmp,&tmp_d);
|
||||
tmp_d_cov_d_r = la::Dot(tmp_d,tmp);
|
||||
la::MulInit(inv,d_cov[i],&inv_d);
|
||||
for(index_t j = 0; j < dim; j++)
|
||||
tmp_d_cov_d_r += inv_d.get(j,j);
|
||||
g_cov_tmp[i] = f*tmp_d_cov_d_r/2;
|
||||
}
|
||||
g_cov->Copy(g_cov_tmp,d_cov.size());
|
||||
|
||||
return f;
|
||||
}
|
||||
@@ -1,26 +0,0 @@
|
||||
|
||||
librule(
|
||||
name = "mog_l2e", # this line can be safely omitted
|
||||
sources = ["mog.cc"], # files that must be compiled
|
||||
headers = ["mog.h","phi.h"],# include files part of the 'lib'
|
||||
deplibs = ["fastlib:fastlib"] # depends on fastlib core
|
||||
#tests = ["mog_l2e_tests.cc"]
|
||||
)
|
||||
|
||||
binrule(
|
||||
name = "mog_l2e_main", # the executable name
|
||||
sources = ["mog_l2e_main.cc"], # compile main.cc
|
||||
headers = ["../opt/optimizers.h"], # no extra headers
|
||||
deplibs = [":mog_l2e","fastlib:fastlib"] #
|
||||
)
|
||||
|
||||
# to build:
|
||||
# 1. make sure have environment variables set up:
|
||||
# $ source /full/path/to/fastlib/script/fl-env /full/path/to/fastlib
|
||||
# (you might want to put this in bashrc)
|
||||
# 2. fl-build main
|
||||
# - this automatically will assume --mode=check, the default
|
||||
# - type fl-build --help for help
|
||||
# 3. ./main
|
||||
# - to build same target again, type: make
|
||||
# - to force recompilation, type: make clean
|
||||
@@ -1,502 +0,0 @@
|
||||
/**
|
||||
* @author Parikshit Ram (pram@cc.gatech.edu)
|
||||
* @file mog.cc
|
||||
*
|
||||
* Implementation for L2 loss function, and
|
||||
* also some initial points generator
|
||||
*
|
||||
*/
|
||||
|
||||
#include "mog.h"
|
||||
#include "phi.h"
|
||||
|
||||
long double MoGL2E::L2Error(const Matrix& data, Vector *gradients) {
|
||||
|
||||
long double reg, fit, l2e;
|
||||
index_t number_of_points = data.n_cols();
|
||||
|
||||
|
||||
if (gradients != NULL) {
|
||||
Vector g_reg, g_fit;
|
||||
reg = RegularizationTerm_(&g_reg);
|
||||
fit = GoodnessOfFitTerm_(data, &g_fit);
|
||||
DEBUG_ASSERT(gradients->length() == (number_of_gaussians()*
|
||||
(dimension()+1)*(dimension()+2)/2 - 1)
|
||||
|
||||
);
|
||||
gradients->SetAll(0.0);
|
||||
la::AddTo(g_reg, gradients);
|
||||
la::AddExpert((-2.0/number_of_points), g_fit, gradients);
|
||||
}
|
||||
else {
|
||||
reg = RegularizationTerm_();
|
||||
fit = GoodnessOfFitTerm_(data);
|
||||
}
|
||||
l2e = reg - 2*fit / number_of_points;
|
||||
return l2e;
|
||||
}
|
||||
|
||||
long double MoGL2E::RegularizationTerm_(Vector *g_reg){
|
||||
Matrix phi_mu, sum_covar;
|
||||
Vector x, y;
|
||||
long double reg, tmpVal;
|
||||
index_t num_gauss, dim;
|
||||
|
||||
Vector df_dw, g_omega;
|
||||
ArrayList<Vector> g_mu, g_sigma;
|
||||
ArrayList<ArrayList<Vector> > dp_d_mu, dp_d_sigma;
|
||||
|
||||
|
||||
num_gauss = number_of_gaussians();
|
||||
dim = dimension();
|
||||
|
||||
phi_mu.Init(num_gauss, num_gauss);
|
||||
sum_covar.Init(dim, dim);
|
||||
x.Copy(omega());
|
||||
|
||||
if (g_reg != NULL) {
|
||||
g_mu.Init(num_gauss);
|
||||
g_sigma.Init(num_gauss);
|
||||
dp_d_mu.Init(num_gauss);
|
||||
dp_d_sigma.Init(num_gauss);
|
||||
for(index_t k = 0; k < num_gauss; k++){
|
||||
dp_d_mu[k].Init(num_gauss);
|
||||
dp_d_sigma[k].Init(num_gauss);
|
||||
}
|
||||
}
|
||||
else {
|
||||
g_mu.Init(0);
|
||||
g_sigma.Init(0);
|
||||
dp_d_mu.Init(0);
|
||||
dp_d_sigma.Init(0);
|
||||
df_dw.Init(0);
|
||||
g_omega.Init(0);
|
||||
}
|
||||
|
||||
for(index_t k = 1; k < num_gauss; k++) {
|
||||
for(index_t j = 0; j < k; j++) {
|
||||
la::AddOverwrite(sigma(k), sigma(j), &sum_covar);
|
||||
|
||||
if (g_reg != NULL) {
|
||||
ArrayList<Matrix> tmp_d_cov;
|
||||
Vector tmp_dp_d_sigma;
|
||||
|
||||
tmp_d_cov.Init(dim*(dim+1));
|
||||
for(index_t i = 0; i < (dim*(dim + 1) / 2); i++){
|
||||
tmp_d_cov[i].Copy(d_sigma(k)[i]);
|
||||
tmp_d_cov[(dim*(dim+1)/2)+i].Copy(d_sigma(j)[i]);
|
||||
}
|
||||
|
||||
//tmpVal = phi(mu(k),mu(j),sum_covar,
|
||||
// tmp_d_cov,&dp_d_mu[j][k],&tmp_dp_d_sigma);
|
||||
tmpVal = phi(mu(k),mu(j),sum_covar,tmp_d_cov, &dp_d_mu[k][j],
|
||||
&tmp_dp_d_sigma);
|
||||
|
||||
phi_mu.set(j, k, tmpVal);
|
||||
phi_mu.set(k, j, tmpVal);
|
||||
|
||||
// la::ScaleInit(-1.0, dp_d_mu[j][k], &dp_d_mu[k][j]);
|
||||
la::ScaleInit(-1.0, dp_d_mu[k][j], &dp_d_mu[j][k]);
|
||||
|
||||
double *tmp_dp, *tmp_dp_1, *tmp_dp_2;
|
||||
tmp_dp = tmp_dp_d_sigma.ptr();
|
||||
tmp_dp_1 = (double*)malloc((tmp_dp_d_sigma.length()/2) * sizeof(double));
|
||||
tmp_dp_2 = (double*)malloc((tmp_dp_d_sigma.length()/2) * sizeof(double));
|
||||
for(index_t i = 0; i < (tmp_dp_d_sigma.length()/2); i++){
|
||||
tmp_dp_1[i] = tmp_dp[i];
|
||||
tmp_dp_2[i] = tmp_dp[(dim*(dim + 1) / 2) + i];
|
||||
}
|
||||
dp_d_sigma[j][k].Copy(tmp_dp_1, (dim*(dim + 1) / 2));
|
||||
dp_d_sigma[k][j].Copy(tmp_dp_2, (dim*(dim + 1) / 2));
|
||||
}
|
||||
else {
|
||||
tmpVal = phi(mu(k), mu(j), sum_covar);
|
||||
phi_mu.set(j, k, tmpVal);
|
||||
phi_mu.set(k, j, tmpVal);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
for(index_t k = 0; k < num_gauss; k++) {
|
||||
la::ScaleOverwrite(2, sigma(k), &sum_covar);
|
||||
|
||||
if (g_reg != NULL) {
|
||||
Vector junk;
|
||||
tmpVal = phi(mu(k), mu(k), sum_covar,
|
||||
d_sigma(k), &junk, &dp_d_sigma[k][k]);
|
||||
phi_mu.set(k, k, tmpVal);
|
||||
dp_d_mu[k][k].Init(dim);
|
||||
dp_d_mu[k][k].SetZero();
|
||||
}
|
||||
else {
|
||||
tmpVal = phi(mu(k), mu(k), sum_covar);
|
||||
phi_mu.set(k, k, tmpVal);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
// Calculating the reg value
|
||||
la::MulInit( x, phi_mu, &y );
|
||||
reg = la::Dot( x, y );
|
||||
|
||||
if (g_reg != NULL) {
|
||||
// Calculating the g_omega values - a vector of size K-1
|
||||
la::ScaleInit(2.0,y,&df_dw);
|
||||
la::MulInit(d_omega(),df_dw,&g_omega);
|
||||
|
||||
// Calculating the g_mu values - K vectors of size D
|
||||
for(index_t k = 0; k < num_gauss; k++){
|
||||
g_mu[k].Init(dim);
|
||||
g_mu[k].SetZero();
|
||||
for(index_t j = 0; j < num_gauss; j++) {
|
||||
la::AddExpert(x.get(j), dp_d_mu[j][k], &g_mu[k]);
|
||||
}
|
||||
la::Scale((2.0 * x.get(k)), &g_mu[k]);
|
||||
}
|
||||
|
||||
// Calculating the g_sigma values - K vectors of size D(D+1)/2
|
||||
for(index_t k = 0; k < num_gauss; k++){
|
||||
g_sigma[k].Init((dim*(dim + 1)) / 2);
|
||||
g_sigma[k].SetZero();
|
||||
for(index_t j = 0; j < num_gauss; j++) {
|
||||
la::AddExpert(x.get(j), dp_d_sigma[j][k], &g_sigma[k]);
|
||||
}
|
||||
la::Scale((2.0 * x.get(k)), &g_sigma[k]);
|
||||
}
|
||||
|
||||
// Making the single gradient vector of size K*(D+1)*(D+2)/2 - 1
|
||||
double *tmp_g_reg;
|
||||
tmp_g_reg = (double*)malloc(((num_gauss*(dim + 1)*(dim + 2) / 2) - 1)
|
||||
*sizeof(double));
|
||||
index_t j = 0;
|
||||
for(index_t k = 0; k < g_omega.length(); k++) {
|
||||
tmp_g_reg[k] = g_omega.get(k);
|
||||
}
|
||||
j = g_omega.length();
|
||||
for(index_t k = 0; k < num_gauss; k++){
|
||||
for(index_t i = 0; i < dim; i++){
|
||||
tmp_g_reg[j + k*(dim) + i] = g_mu[k].get(i);
|
||||
}
|
||||
for(index_t i = 0; i < (dim*(dim+1)/2); i++){
|
||||
tmp_g_reg[j + num_gauss*dim
|
||||
+ k*(dim*(dim+1) / 2)
|
||||
+ i] = g_sigma[k].get(i);
|
||||
}
|
||||
}
|
||||
g_reg->Copy(tmp_g_reg, ((num_gauss*(dim+1)*(dim+2) / 2) - 1));
|
||||
}
|
||||
|
||||
return reg;
|
||||
}
|
||||
|
||||
long double MoGL2E::GoodnessOfFitTerm_(const Matrix& data, Vector *g_fit) {
|
||||
long double fit;
|
||||
Matrix phi_x;
|
||||
Vector weights, x, y, identity_vector;
|
||||
index_t num_gauss, num_points, dim;
|
||||
long double tmpVal;
|
||||
Vector g_omega,tmp_g_omega;
|
||||
ArrayList<Vector> g_mu, g_sigma;
|
||||
|
||||
num_gauss = number_of_gaussians();
|
||||
num_points = data.n_cols();
|
||||
dim = data.n_rows();
|
||||
phi_x.Init(num_gauss, num_points);
|
||||
weights.Copy(omega());
|
||||
x.Init(data.n_rows());
|
||||
identity_vector.Init(num_points);
|
||||
identity_vector.SetAll(1);
|
||||
|
||||
if(g_fit != NULL) {
|
||||
g_mu.Init(num_gauss);
|
||||
g_sigma.Init(num_gauss);
|
||||
}
|
||||
else {
|
||||
g_mu.Init(0);
|
||||
g_sigma.Init(0);
|
||||
g_omega.Init(0);
|
||||
tmp_g_omega.Init(0);
|
||||
}
|
||||
|
||||
for(index_t k = 0; k < num_gauss; k++) {
|
||||
if (g_fit != NULL) {
|
||||
g_mu[k].Init(dim);
|
||||
g_mu[k].SetZero();
|
||||
g_sigma[k].Init((dim * (dim+1) / 2));
|
||||
g_sigma[k].SetZero();
|
||||
}
|
||||
for(index_t i = 0; i < num_points; i++) {
|
||||
if (g_fit != NULL) {
|
||||
Vector tmp_g_mu, tmp_g_sigma;
|
||||
x.CopyValues(data.GetColumnPtr(i));
|
||||
tmpVal = phi(x, mu(k), sigma(k),
|
||||
d_sigma(k), &tmp_g_mu, &tmp_g_sigma);
|
||||
phi_x.set(k, i, tmpVal);
|
||||
la::AddTo(tmp_g_mu, &g_mu[k]);
|
||||
la::AddTo(tmp_g_sigma, &g_sigma[k]);
|
||||
}
|
||||
else {
|
||||
x.CopyValues(data.GetColumnPtr(i));
|
||||
phi_x.set(k, i, phi(x, mu(k), sigma(k)));
|
||||
}
|
||||
}
|
||||
if (g_fit != NULL) {
|
||||
la::Scale(weights.get(k), &g_mu[k]);
|
||||
la::Scale(weights.get(k), &g_sigma[k]);
|
||||
}
|
||||
}
|
||||
|
||||
la::MulInit(weights, phi_x, &y);
|
||||
fit = la::Dot(y, identity_vector);
|
||||
|
||||
if (g_fit != NULL) {
|
||||
// Calculating the g_omega
|
||||
la::MulInit(phi_x, identity_vector, &tmp_g_omega);
|
||||
la::MulInit(d_omega(), tmp_g_omega, &g_omega);
|
||||
|
||||
// Making the single gradient vector of size K*(D+1)*(D+2)/2
|
||||
double *tmp_g_fit;
|
||||
tmp_g_fit = (double*)malloc(((num_gauss * (dim+1)*(dim+2) / 2) - 1)
|
||||
*sizeof(double));
|
||||
index_t j = 0;
|
||||
for(index_t k = 0; k < g_omega.length(); k++)
|
||||
tmp_g_fit[k] = g_omega.get(k);
|
||||
j = g_omega.length();
|
||||
for(index_t k = 0; k < num_gauss; k++){
|
||||
for(index_t i = 0; i < dim; i++){
|
||||
tmp_g_fit[j + k*dim + i] = g_mu[k].get(i);
|
||||
}
|
||||
for(index_t i = 0; i < (dim * (dim+1) / 2); i++){
|
||||
tmp_g_fit[j + num_gauss*dim
|
||||
+ k*(dim * (dim+1) / 2)
|
||||
+ i] = g_sigma[k].get(i);
|
||||
}
|
||||
}
|
||||
g_fit->Copy(tmp_g_fit, ((num_gauss*(dim+1)*(dim+2) / 2) - 1));
|
||||
}
|
||||
return fit;
|
||||
}
|
||||
|
||||
void MoGL2E::MultiplePointsGenerator(double **points,
|
||||
index_t number_of_points,
|
||||
const Matrix& d,
|
||||
index_t number_of_components) {
|
||||
|
||||
index_t dim, n, i, j, x;
|
||||
|
||||
dim = d.n_rows();
|
||||
n = d.n_cols();
|
||||
|
||||
for( i = 0; i < number_of_points; i++) {
|
||||
for(j = 0; j < number_of_components - 1; j++) {
|
||||
points[i][j] = (rand() % 20001)/1000 - 10;
|
||||
}
|
||||
}
|
||||
|
||||
for(i = 0; i < number_of_points; i++){
|
||||
for(j = 0; j < number_of_components; j++){
|
||||
Vector tmp_mu;
|
||||
tmp_mu.Init(dim);
|
||||
tmp_mu.CopyValues(d.GetColumnPtr((rand() % n)));
|
||||
for(x = 0; x < dim; x++)
|
||||
points[i][number_of_components - 1 + j * dim + x] = tmp_mu.get(x);
|
||||
}
|
||||
}
|
||||
|
||||
for(i = 0; i < number_of_points; i++)
|
||||
for(j = 0; j < number_of_components; j++)
|
||||
for(x = 0 ; x < (dim * (dim + 1) / 2); x++)
|
||||
points[i][(number_of_components * (dim + 1) - 1)
|
||||
+ (j * (dim * (dim + 1) / 2)) + x] = (rand() % 501)/100;
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
void MoGL2E::InitialPointGenerator(double *theta, const Matrix& data,
|
||||
index_t k_comp) {
|
||||
|
||||
ArrayList<Vector> means;
|
||||
ArrayList<Matrix> covars;
|
||||
Vector weights;
|
||||
double temp, noise;
|
||||
index_t dim;
|
||||
|
||||
weights.Init(k_comp);
|
||||
means.Init(k_comp);
|
||||
covars.Init(k_comp);
|
||||
dim = data.n_rows();
|
||||
|
||||
for (index_t i = 0; i < k_comp; i++) {
|
||||
means[i].Init(dim);
|
||||
covars[i].Init(dim, dim);
|
||||
}
|
||||
|
||||
KMeans_(data, &means, &covars, &weights, k_comp);
|
||||
|
||||
for(index_t k = 0; k < k_comp - 1; k++){
|
||||
temp = weights[k] / weights[k_comp - 1];
|
||||
noise = (double)(rand() % 10000) / (double)1000;
|
||||
theta[k] = noise - 5;
|
||||
}
|
||||
for(index_t k = 0; k < k_comp; k++){
|
||||
for(index_t j = 0; j < dim; j++)
|
||||
theta[k_comp - 1 + k * dim + j] = means[k].get(j);
|
||||
|
||||
Matrix U, U_tran;
|
||||
la::CholeskyInit(covars[k], &U);
|
||||
la::TransposeInit(U, &U_tran);
|
||||
for(index_t j = 0; j < dim; j++) {
|
||||
for(index_t i = 0; i < j + 1; i++) {
|
||||
noise = (rand() % 501) / 100;
|
||||
theta[k_comp - 1 + k_comp * dim
|
||||
+ k * dim * (dim + 1) / 2
|
||||
+ j * (j + 1) / 2 + i] = U_tran.get(j, i) + noise;
|
||||
}
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
void MoGL2E::KMeans_(const Matrix& data, ArrayList<Vector> *means,
|
||||
ArrayList<Matrix> *covars, Vector *weights,
|
||||
index_t value_of_k){
|
||||
|
||||
ArrayList<Vector> mu, mu_old;
|
||||
double* tmpssq;
|
||||
double* sig;
|
||||
double* sig_best;
|
||||
index_t *y;
|
||||
Vector x, diff;
|
||||
Matrix ssq;
|
||||
index_t i, j, k, n, t, dim;
|
||||
double score, score_old, sum;
|
||||
|
||||
n = data.n_cols();
|
||||
dim = data.n_rows();
|
||||
mu.Init(value_of_k);
|
||||
mu_old.Init(value_of_k);
|
||||
tmpssq = (double*)malloc(value_of_k * sizeof( double ));
|
||||
sig = (double*)malloc(value_of_k * sizeof( double ));
|
||||
sig_best = (double*)malloc(value_of_k * sizeof( double ));
|
||||
ssq.Init(n, value_of_k);
|
||||
|
||||
for( i = 0; i < value_of_k; i++){
|
||||
mu[i].Init(dim);
|
||||
mu_old[i].Init(dim);
|
||||
}
|
||||
x.Init(dim);
|
||||
y = (index_t*)malloc(n * sizeof(index_t));
|
||||
diff.Init(dim);
|
||||
|
||||
score_old = 999999;
|
||||
|
||||
// putting 5 random restarts to obtain the k-means
|
||||
for(i = 0; i < 5; i++){
|
||||
t = -1;
|
||||
for (k = 0; k < value_of_k; k++){
|
||||
t = (t + 1 + (rand()%((n - 1 - (value_of_k - k)) - (t + 1))));
|
||||
mu[k].CopyValues(data.GetColumnPtr(t));
|
||||
for(j = 0; j < n; j++){
|
||||
x.CopyValues( data.GetColumnPtr(j));
|
||||
la::SubOverwrite(mu[k], x, &diff);
|
||||
ssq.set( j, k, la::Dot(diff, diff));
|
||||
}
|
||||
}
|
||||
min_element(ssq, y);
|
||||
|
||||
do{
|
||||
for(k = 0; k < value_of_k; k++){
|
||||
mu_old[k].CopyValues(mu[k]);
|
||||
}
|
||||
|
||||
for(k = 0; k < value_of_k; k++){
|
||||
index_t p = 0;
|
||||
mu[k].SetZero();
|
||||
for(j = 0; j < n; j++){
|
||||
x.CopyValues(data.GetColumnPtr(j));
|
||||
if(y[j] == k){
|
||||
la::AddTo(x, &mu[k]);
|
||||
p++;
|
||||
}
|
||||
}
|
||||
|
||||
if(p == 0){
|
||||
}
|
||||
else{
|
||||
double sc = 1 ;
|
||||
sc = sc / p;
|
||||
la::Scale(sc , &mu[k]);
|
||||
}
|
||||
for(j = 0; j < n; j++){
|
||||
x.CopyValues(data.GetColumnPtr(j));
|
||||
la::SubOverwrite(mu[k], x, &diff);
|
||||
ssq.set(j, k, la::Dot(diff, diff));
|
||||
}
|
||||
}
|
||||
min_element(ssq, y);
|
||||
|
||||
sum = 0;
|
||||
for(k = 0; k < value_of_k; k++) {
|
||||
la::SubOverwrite(mu[k], mu_old[k], &diff);
|
||||
sum += la::Dot(diff, diff);
|
||||
}
|
||||
}while(sum != 0);
|
||||
|
||||
for(k = 0; k < value_of_k; k++){
|
||||
index_t p = 0;
|
||||
tmpssq[k] = 0;
|
||||
for(j = 0; j < n; j++){
|
||||
if(y[j] == k){
|
||||
tmpssq[k] += ssq.get(j, k);
|
||||
p++;
|
||||
}
|
||||
}
|
||||
sig[k] = sqrt(tmpssq[k] / p);
|
||||
}
|
||||
|
||||
score = 0;
|
||||
for(k = 0; k < value_of_k; k++){
|
||||
score += tmpssq[k];
|
||||
}
|
||||
score = score / n;
|
||||
|
||||
if (score < score_old) {
|
||||
score_old = score;
|
||||
for(k = 0; k < value_of_k; k++){
|
||||
(*means)[k].CopyValues(mu[k]);
|
||||
sig_best[k] = sig[k];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
for(k = 0; k < value_of_k; k++){
|
||||
x.SetAll(sig_best[k]);
|
||||
(*covars)[k].SetDiagonal(x);
|
||||
}
|
||||
double tmp = 1;
|
||||
weights->SetAll(tmp / value_of_k);
|
||||
return;
|
||||
}
|
||||
|
||||
void MoGL2E::min_element( Matrix& element, index_t *indices ){
|
||||
|
||||
index_t last = element.n_cols() - 1;
|
||||
index_t first, lowest;
|
||||
index_t i;
|
||||
|
||||
for( i = 0; i < element.n_rows(); i++ ){
|
||||
|
||||
first = lowest = 0;
|
||||
if(first == last){
|
||||
indices[ i ] = last;
|
||||
}
|
||||
while(++first <= last){
|
||||
if( element.get( i , first ) < element.get( i , lowest ) ){
|
||||
lowest = first;
|
||||
}
|
||||
}
|
||||
indices[ i ] = lowest;
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
@@ -1,595 +0,0 @@
|
||||
/**
|
||||
* @author Parikshit Ram (pram@cc.gatech.edu)
|
||||
* @file mog.h
|
||||
*
|
||||
* Defines a Gaussian Mixture model and
|
||||
* estimates the parameters of the model
|
||||
*
|
||||
*/
|
||||
|
||||
#ifndef MOGL2E_H
|
||||
#define MOGL2E_H
|
||||
|
||||
#include <fastlib/fastlib.h>
|
||||
|
||||
|
||||
/**
|
||||
* A Gaussian mixture model class.
|
||||
*
|
||||
* This class uses L2 loss function to
|
||||
* estimate the parameters of a gaussian mixture
|
||||
* model on a given data.
|
||||
*
|
||||
* The parameters are converted for optimization
|
||||
* to maintain the following facts:
|
||||
* - the weights sum to one
|
||||
* - for this, the weights were parameterized using
|
||||
* the logistic function
|
||||
* - the covariance matrix is always positive definite
|
||||
* - for this, the Cholesky decomposition is used
|
||||
*
|
||||
* Example use:
|
||||
*
|
||||
* @code
|
||||
* MoGL2E mog;
|
||||
* ArrayList<double> results;
|
||||
* double *params;
|
||||
*
|
||||
* mog.MakeModel(number_of_gaussians, dimension, params);
|
||||
* mog.L2Error(data);
|
||||
* mog.OutputResults(&results);
|
||||
* @endcode
|
||||
*/
|
||||
class MoGL2E {
|
||||
|
||||
private:
|
||||
|
||||
// The parameters of the Mixture model
|
||||
ArrayList<Vector> mu_;
|
||||
ArrayList<Matrix> sigma_;
|
||||
Vector omega_;
|
||||
index_t number_of_gaussians_;
|
||||
index_t dimension_;
|
||||
|
||||
// The differential for the paramterization
|
||||
// for optimization
|
||||
Matrix d_omega_;
|
||||
ArrayList<ArrayList<Matrix> > d_sigma_;
|
||||
|
||||
public:
|
||||
|
||||
MoGL2E() {
|
||||
mu_.Init(0);
|
||||
sigma_.Init(0);
|
||||
d_sigma_.Init(0);
|
||||
d_omega_.Init(0, 0);
|
||||
}
|
||||
|
||||
~MoGL2E() {
|
||||
}
|
||||
|
||||
void Init(index_t num_gauss, index_t dimension) {
|
||||
|
||||
// Destruct everything to initialize afresh
|
||||
mu_.Clear();
|
||||
sigma_.Clear();
|
||||
d_sigma_.Clear();
|
||||
// Initialize the private variables
|
||||
number_of_gaussians_ = num_gauss;
|
||||
dimension_ = dimension;
|
||||
|
||||
// Resize the ArrayList of Vectors and Matrices
|
||||
mu_.Resize(number_of_gaussians_);
|
||||
sigma_.Resize(number_of_gaussians_);
|
||||
}
|
||||
|
||||
void Resize_d_sigma_() {
|
||||
d_sigma_.Resize(number_of_gaussians());
|
||||
for(index_t i =0; i < number_of_gaussians(); i++) {
|
||||
d_sigma_[i].Init(dimension()*(dimension()+1)/2);
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
*
|
||||
* This function uses the parameters used for optimization
|
||||
* and converts it into athe parameters of a Gaussian
|
||||
* mixture model. This is to be used when you do not want
|
||||
* the gradient values.
|
||||
*
|
||||
* Example use:
|
||||
*
|
||||
* @code
|
||||
* MoGL2E mog;
|
||||
* mog.MakeModel(number_of_gaussians, dimension,
|
||||
* parameters_for_optimization);
|
||||
* @endcode
|
||||
*/
|
||||
|
||||
void MakeModel(index_t num_mods, index_t dimension, double* theta) {
|
||||
double *temp_mu;
|
||||
Matrix lower_triangle_matrix, upper_triangle_matrix;
|
||||
double sum, s_min = 0.01;
|
||||
|
||||
Init(num_mods, dimension);
|
||||
temp_mu = (double*) malloc (dimension * sizeof(double)) ;
|
||||
lower_triangle_matrix.Init(dimension, dimension);
|
||||
upper_triangle_matrix.Init(dimension, dimension);
|
||||
|
||||
// calculating the omega values
|
||||
sum = 0;
|
||||
double *temp_array;
|
||||
temp_array = (double*) malloc (num_mods * sizeof(double));
|
||||
for(index_t i = 0; i < num_mods - 1; i++) {
|
||||
temp_array[i] = exp(theta[i]) ;
|
||||
sum += temp_array[i] ;
|
||||
}
|
||||
temp_array[num_mods - 1] = 1 ;
|
||||
++sum ;
|
||||
la::Scale(num_mods, (1.0 / sum), temp_array);
|
||||
set_omega(dimension, temp_array);
|
||||
|
||||
// calculating the mu values
|
||||
for(index_t k = 0; k < num_mods; k++) {
|
||||
for(index_t j = 0; j < dimension; j++) {
|
||||
temp_mu[j] = theta[num_mods + k*dimension + j - 1];
|
||||
}
|
||||
set_mu(k, dimension, temp_mu);
|
||||
}
|
||||
|
||||
// calculating the sigma values
|
||||
// using a lower triangular matrix and its transpose
|
||||
// to obtain a positive definite symmetric matrix
|
||||
Matrix sigma_temp;
|
||||
sigma_temp.Init(dimension, dimension);
|
||||
for(index_t k = 0; k < num_mods; k++) {
|
||||
lower_triangle_matrix.SetAll(0.0);
|
||||
for(index_t j = 0; j < dimension; j++) {
|
||||
for(index_t i = 0; i < j; i++) {
|
||||
lower_triangle_matrix.set(j, i,
|
||||
theta[(num_mods - 1)
|
||||
+ num_mods*dimension
|
||||
+ k*(dimension*(dimension + 1) / 2)
|
||||
+ (j*(j + 1) / 2) + i]);
|
||||
}
|
||||
lower_triangle_matrix.set(j, j,
|
||||
theta[(num_mods - 1)
|
||||
+ num_mods*dimension
|
||||
+ k*(dimension*(dimension + 1) / 2)
|
||||
+ (j*(j + 1) / 2) + j] + s_min);
|
||||
|
||||
}
|
||||
la::TransposeOverwrite(lower_triangle_matrix, &upper_triangle_matrix);
|
||||
la::MulOverwrite(lower_triangle_matrix, upper_triangle_matrix, &sigma_temp);
|
||||
set_sigma(k, sigma_temp);
|
||||
}
|
||||
}
|
||||
|
||||
void MakeModel(datanode *mog_l2e_module, double* theta) {
|
||||
index_t num_gauss = fx_param_int_req(mog_l2e_module, "K");
|
||||
index_t dimension = fx_param_int_req(mog_l2e_module, "D");
|
||||
MakeModel(num_gauss, dimension, theta);
|
||||
}
|
||||
|
||||
/**
|
||||
*
|
||||
* This function uses the parameters used for optimization
|
||||
* and converts it into athe parameters of a Gaussian
|
||||
* mixture model. This is to be used when you want
|
||||
* the gradient values.
|
||||
*
|
||||
* Example use:
|
||||
*
|
||||
* @code
|
||||
* MoGL2E mog;
|
||||
* mog.MakeModelWithGradients(number_of_gaussians, dimension,
|
||||
* parameters_for_optimization);
|
||||
* @endcode
|
||||
*/
|
||||
|
||||
void MakeModelWithGradients(index_t num_mods, index_t dimension, double* theta) {
|
||||
double *temp_mu;
|
||||
Matrix lower_triangle_matrix, upper_triangle_matrix;
|
||||
double sum, s_min = 0.01;
|
||||
|
||||
Init(num_mods, dimension);
|
||||
temp_mu = (double*) malloc (dimension * sizeof(double));
|
||||
lower_triangle_matrix.Init(dimension, dimension);
|
||||
upper_triangle_matrix.Init(dimension, dimension);
|
||||
|
||||
// calculating the omega values
|
||||
sum = 0;
|
||||
double *temp_array;
|
||||
temp_array = (double*) malloc (num_mods * sizeof(double));
|
||||
for(index_t i = 0; i < num_mods - 1; i++) {
|
||||
temp_array[i] = exp(theta[i]) ;
|
||||
sum += temp_array[i] ;
|
||||
}
|
||||
temp_array[num_mods - 1] = 1 ;
|
||||
++sum ;
|
||||
la::Scale(num_mods, (1.0 / sum), temp_array);
|
||||
set_omega(num_mods, temp_array);
|
||||
|
||||
// calculating the d_omega values
|
||||
Matrix d_omega_temp;
|
||||
d_omega_temp.Init(num_mods - 1, num_mods);
|
||||
d_omega_temp.SetAll(0.0);
|
||||
for(index_t i = 0; i < num_mods - 1; i++) {
|
||||
for(index_t j = 0; j < i; j++) {
|
||||
d_omega_temp.set(i,j,-(omega(i)*omega(j)));
|
||||
d_omega_temp.set(j,i,-(omega(i)*omega(j)));
|
||||
}
|
||||
d_omega_temp.set(i,i,omega(i)*(1-omega(i)));
|
||||
}
|
||||
for(index_t i = 0; i < num_mods - 1; i++) {
|
||||
d_omega_temp.set(i, num_mods - 1, -(omega(i)*omega(num_mods - 1)));
|
||||
}
|
||||
set_d_omega(d_omega_temp);
|
||||
|
||||
// calculating the mu values
|
||||
for(index_t k = 0; k < num_mods; k++) {
|
||||
for(index_t j = 0; j < dimension; j++) {
|
||||
temp_mu[j] = theta[num_mods + k*dimension + j - 1];
|
||||
}
|
||||
set_mu(k, dimension, temp_mu);
|
||||
}
|
||||
// d_mu is not computed because it is implicitly known
|
||||
// since no parameterization is applied on them
|
||||
|
||||
// using a lower triangular matrix and its transpose
|
||||
// to obtain a positive definite symmetric matrix
|
||||
|
||||
// initializing the d_sigma values
|
||||
|
||||
Matrix d_sigma_temp;
|
||||
d_sigma_temp.Init(dimension, dimension);
|
||||
Resize_d_sigma_();
|
||||
|
||||
// calculating the sigma values
|
||||
Matrix sigma_temp;
|
||||
sigma_temp.Init(dimension, dimension);
|
||||
for(index_t k = 0; k < num_mods; k++) {
|
||||
lower_triangle_matrix.SetAll(0.0);
|
||||
for(index_t j = 0; j < dimension; j++) {
|
||||
for(index_t i = 0; i < j; i++) {
|
||||
lower_triangle_matrix.set( j, i,
|
||||
theta[(num_mods - 1)
|
||||
+ num_mods*dimension
|
||||
+ k*(dimension*(dimension + 1) / 2)
|
||||
+ (j*(j + 1) / 2) + i]) ;
|
||||
}
|
||||
lower_triangle_matrix.set(j, j,
|
||||
theta[(num_mods - 1)
|
||||
+ num_mods*dimension
|
||||
+ k*(dimension*(dimension + 1) / 2)
|
||||
+ (j*(j + 1) / 2) + j] + s_min);
|
||||
}
|
||||
la::TransposeOverwrite(lower_triangle_matrix, &upper_triangle_matrix);
|
||||
la::MulOverwrite(lower_triangle_matrix, upper_triangle_matrix, &sigma_temp);
|
||||
set_sigma(k, sigma_temp);
|
||||
|
||||
// calculating the d_sigma values
|
||||
for(index_t i = 0; i < dimension; i++){
|
||||
for(index_t in = 0; in < i+1; in++){
|
||||
Matrix d_sigma_d_r,d_sigma_d_r_t,temp_matrix_1,temp_matrix_2;
|
||||
d_sigma_d_r.Init(dimension, dimension);
|
||||
d_sigma_d_r_t.Init(dimension, dimension);
|
||||
d_sigma_d_r.SetAll(0.0);
|
||||
d_sigma_d_r_t.SetAll(0.0);
|
||||
d_sigma_d_r.set(i,in,1.0);
|
||||
d_sigma_d_r_t.set(in,i,1.0);
|
||||
|
||||
la::MulInit(d_sigma_d_r,upper_triangle_matrix,&temp_matrix_1);
|
||||
la::MulInit(lower_triangle_matrix,d_sigma_d_r_t,&temp_matrix_2);
|
||||
la::AddOverwrite(temp_matrix_1,temp_matrix_2,&d_sigma_temp);
|
||||
set_d_sigma(k, (i*(i+1)/2)+in, d_sigma_temp);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
////// THE GET FUNCTIONS //////
|
||||
ArrayList<Vector>& mu() {
|
||||
return mu_;
|
||||
}
|
||||
|
||||
ArrayList<Matrix>& sigma() {
|
||||
return sigma_;
|
||||
}
|
||||
|
||||
Vector& omega() {
|
||||
return omega_;
|
||||
}
|
||||
|
||||
index_t number_of_gaussians() {
|
||||
return number_of_gaussians_;
|
||||
}
|
||||
|
||||
index_t dimension() {
|
||||
return dimension_;
|
||||
}
|
||||
|
||||
Vector& mu(index_t i) {
|
||||
return mu_[i] ;
|
||||
}
|
||||
|
||||
Matrix& sigma(index_t i) {
|
||||
return sigma_[i];
|
||||
}
|
||||
|
||||
double omega(index_t i) {
|
||||
return omega_.get(i);
|
||||
}
|
||||
|
||||
Matrix& d_omega(){
|
||||
return d_omega_;
|
||||
}
|
||||
|
||||
ArrayList<ArrayList<Matrix> >& d_sigma(){
|
||||
return d_sigma_;
|
||||
}
|
||||
|
||||
ArrayList<Matrix>& d_sigma(index_t i){
|
||||
return d_sigma_[i];
|
||||
}
|
||||
|
||||
////// THE SET FUNCTIONS //////
|
||||
|
||||
void set_mu(index_t i, Vector& mu) {
|
||||
DEBUG_ASSERT(i < number_of_gaussians());
|
||||
DEBUG_ASSERT(mu.length() == dimension());
|
||||
mu_[i].Copy(mu);
|
||||
return;
|
||||
}
|
||||
|
||||
void set_mu(index_t i, index_t length, const double *mu) {
|
||||
DEBUG_ASSERT(i < number_of_gaussians());
|
||||
DEBUG_ASSERT(length == dimension());
|
||||
mu_[i].Copy(mu, length);
|
||||
return;
|
||||
}
|
||||
|
||||
void set_sigma(index_t i, Matrix& sigma) {
|
||||
DEBUG_ASSERT(i < number_of_gaussians());
|
||||
DEBUG_ASSERT(sigma.n_rows() == dimension());
|
||||
DEBUG_ASSERT(sigma.n_cols() == dimension());
|
||||
sigma_[i].Copy(sigma);
|
||||
return;
|
||||
}
|
||||
|
||||
void set_omega(Vector& omega) {
|
||||
DEBUG_ASSERT(omega.length() == number_of_gaussians());
|
||||
omega_.Copy(omega);
|
||||
return;
|
||||
}
|
||||
|
||||
void set_omega(index_t length, const double *omega) {
|
||||
DEBUG_ASSERT(length == number_of_gaussians());
|
||||
omega_.Copy(omega, length);
|
||||
return;
|
||||
}
|
||||
|
||||
void set_d_omega(Matrix& d_omega) {
|
||||
d_omega_.Destruct();
|
||||
d_omega_.Copy(d_omega);
|
||||
return;
|
||||
}
|
||||
|
||||
void set_d_sigma(index_t i, index_t j, Matrix& d_sigma_i_j) {
|
||||
d_sigma_[i][j].Copy(d_sigma_i_j);
|
||||
return;
|
||||
}
|
||||
|
||||
/**
|
||||
* This function outputs the parameters of the model
|
||||
* to an arraylist of doubles
|
||||
*
|
||||
* @code
|
||||
* ArrayList<double> results;
|
||||
* mog.OutputResults(&results);
|
||||
* @endcode
|
||||
*/
|
||||
void OutputResults(ArrayList<double> *results) {
|
||||
|
||||
// Initialize the size of the output array
|
||||
(*results).Init(number_of_gaussians_ * (1 + dimension_*(1 + dimension_)));
|
||||
|
||||
// Copy values to the array from the private variables of the class
|
||||
for (index_t i = 0; i < number_of_gaussians_; i++) {
|
||||
(*results)[i] = omega(i);
|
||||
for (index_t j = 0; j < dimension_; j++) {
|
||||
(*results)[number_of_gaussians_ + i*dimension_ + j] = mu(i).get(j);
|
||||
for (index_t k = 0; k < dimension_; k++) {
|
||||
(*results)[number_of_gaussians_*(1 + dimension_)
|
||||
+ i*dimension_*dimension_ + j*dimension_
|
||||
+ k] = sigma(i).get(j, k);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* This function prints the parameters of the model
|
||||
*
|
||||
* @code
|
||||
* mog.Display();
|
||||
* @endcode
|
||||
*/
|
||||
void Display(){
|
||||
|
||||
// Output the model parameters as the omega, mu and sigma
|
||||
printf(" Omega : [ ");
|
||||
for (index_t i = 0; i < number_of_gaussians_; i++) {
|
||||
printf("%lf ", omega(i));
|
||||
}
|
||||
printf("]\n");
|
||||
printf(" Mu : \n[");
|
||||
for (index_t i = 0; i < number_of_gaussians_; i++) {
|
||||
for (index_t j = 0; j < dimension_ ; j++) {
|
||||
printf("%lf ", mu(i).get(j));
|
||||
}
|
||||
printf(";");
|
||||
if (i == (number_of_gaussians_ - 1)) {
|
||||
printf("\b]\n");
|
||||
}
|
||||
}
|
||||
printf("Sigma : ");
|
||||
for (index_t i = 0; i < number_of_gaussians_; i++) {
|
||||
printf("\n[");
|
||||
for (index_t j = 0; j < dimension_ ; j++) {
|
||||
for(index_t k = 0; k < dimension_ ; k++) {
|
||||
printf("%lf ",sigma(i).get(j, k));
|
||||
}
|
||||
printf(";");
|
||||
}
|
||||
printf("\b]");
|
||||
}
|
||||
printf("\n");
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* This function calculates the L2 error and
|
||||
* the gradient of the error with respect to the
|
||||
* parameters given the data and the parameterized
|
||||
* mixture
|
||||
*
|
||||
* Example use:
|
||||
*
|
||||
* @code
|
||||
* const Matrix data;
|
||||
* MoGL2E mog;
|
||||
* index_t num_gauss, dimension;
|
||||
* double *params; // get the parameters
|
||||
*
|
||||
* mog.MakeModel(num_gauss, dimension, params);
|
||||
* mog.L2Error(data);
|
||||
* @endcode
|
||||
*/
|
||||
long double L2Error(const Matrix&, Vector* = NULL);
|
||||
|
||||
/**
|
||||
* Calculates the regularization value for a
|
||||
* Gaussian mixture and its gradient with
|
||||
* respect to the parameters
|
||||
*
|
||||
* Used by the 'L2Error' function to calculate
|
||||
* the regularization part of the error
|
||||
*/
|
||||
long double RegularizationTerm_(Vector* = NULL);
|
||||
|
||||
/**
|
||||
* Calculates the goodness-of-fit value for a
|
||||
* Gaussian mixture and its gradient with
|
||||
* respect to the parameters
|
||||
*
|
||||
* Used by the 'L2Error' function to calculate
|
||||
* the goodness-of-fit part of the error
|
||||
*/
|
||||
long double GoodnessOfFitTerm_(const Matrix&, Vector* = NULL);
|
||||
|
||||
/**
|
||||
* This function computes multiple number of starting points
|
||||
* required for the Nelder Mead method
|
||||
*
|
||||
* Example use:
|
||||
* @code
|
||||
* double **p;
|
||||
* index_t n, num_gauss;
|
||||
* const Matrix data;
|
||||
*
|
||||
* MoGL2E::MultiplePointsGeneratot(p, n, data, num_gauss);
|
||||
* @endcode
|
||||
*/
|
||||
static void MultiplePointsGenerator(double**, index_t,
|
||||
const Matrix&, index_t);
|
||||
|
||||
/**
|
||||
* This function parameterizes the starting point obtained
|
||||
* from the 'k_means" for optimization purposes using the
|
||||
* Quasi Newton method
|
||||
*
|
||||
* Example use:
|
||||
* @code
|
||||
* double *p;
|
||||
* index_t num_gauss;
|
||||
* const Matrix data;
|
||||
*
|
||||
* MoGL2E::InitialPointGeneratot(p, data, num_gauss);
|
||||
* @endcode
|
||||
*/
|
||||
static void InitialPointGenerator(double*, const Matrix&, index_t);
|
||||
|
||||
/**
|
||||
* This function computes the k-means of the data and stores
|
||||
* the calculated means and covariances in the ArrayList
|
||||
* of Vectors and Matrices passed to it. It sets the weights
|
||||
* uniformly.
|
||||
*
|
||||
* This function is used to obtain a starting point for
|
||||
* the optimization
|
||||
*
|
||||
* Example use:
|
||||
*
|
||||
* @code
|
||||
* const Matrix data;
|
||||
* ArrayList<Vector> *means;
|
||||
* ArrayList<Matrix> *covars;
|
||||
* Vector *weights;
|
||||
* index_t num_gauss;
|
||||
*
|
||||
* ...
|
||||
* MoGL2E::KMeans(data, means, covars, weights, num_gauss);
|
||||
*@endcode
|
||||
*/
|
||||
static void KMeans_(const Matrix&, ArrayList<Vector>*,
|
||||
ArrayList<Matrix>*, Vector*, index_t);
|
||||
|
||||
/**
|
||||
* This function returns the indices of the minimum
|
||||
* element in each row of a matrix
|
||||
*/
|
||||
static void min_element(Matrix&, index_t*);
|
||||
|
||||
/**
|
||||
* This is the function which would be used for
|
||||
* optimization. It creates its own object of
|
||||
* class MoGL2E and returns the L2 error
|
||||
* and the gradient which are computed by
|
||||
* the functions of the class
|
||||
*
|
||||
*/
|
||||
static long double L2ErrorForOpt(Vector& params,
|
||||
const Matrix& data,
|
||||
Vector *gradient) {
|
||||
|
||||
MoGL2E model;
|
||||
index_t dimension = data.n_rows();
|
||||
index_t num_gauss;
|
||||
|
||||
num_gauss = (params.length() + 1)*2 / ((dimension+1)*(dimension+2));
|
||||
// This check added here to see if
|
||||
// the gradient is actually demanded here
|
||||
if (gradient != NULL) {
|
||||
model.MakeModelWithGradients(num_gauss, dimension, params.ptr());
|
||||
return model.L2Error(data, gradient);
|
||||
}
|
||||
else {
|
||||
model.MakeModel(num_gauss, dimension, params.ptr());
|
||||
return model.L2Error(data);
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* This is the function which should be used for
|
||||
* optimization when there is no need to compute
|
||||
* any gradients
|
||||
*
|
||||
*/
|
||||
static long double L2ErrorForOpt(Vector& params, const Matrix& data) {
|
||||
return L2ErrorForOpt(params, data, NULL);
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
#endif
|
||||
File diff suppressed because it is too large
Load Diff
@@ -1,145 +0,0 @@
|
||||
/**
|
||||
* @author Parikshit Ram (pram@cc.gatech.edu)
|
||||
* @file mog_l2e_main.cc
|
||||
*
|
||||
* This program test drives the L2 estimation
|
||||
* of a Gaussian Mixture model.
|
||||
*
|
||||
* PARAMETERS TO BE INPUT:
|
||||
*
|
||||
* --data
|
||||
* This is the file that contains the data on which
|
||||
* the model is to be fit
|
||||
*
|
||||
* --mog_l2e/K
|
||||
* This is the number of gaussians we want to fit
|
||||
* on the data, defaults to '1'
|
||||
*
|
||||
* --output
|
||||
* This file will contain the parameters estimated,
|
||||
* defaults to 'output.csv'
|
||||
*
|
||||
*/
|
||||
|
||||
#include "mog.h"
|
||||
#include "../opt/optimizers.h"
|
||||
//#include <time.h>
|
||||
|
||||
|
||||
int main(int argc, char* argv[]) {
|
||||
|
||||
fx_init(argc, argv);
|
||||
//srand(time(NULL));
|
||||
////// READING PARAMETERS AND LOADING DATA //////
|
||||
|
||||
const char *data_filename = fx_param_str_req(NULL, "data");
|
||||
|
||||
Matrix data_points;
|
||||
data::Load(data_filename, &data_points);
|
||||
|
||||
////// MIXTURE OF GAUSSIANS USING L2 ESTIMATION //////
|
||||
|
||||
datanode *mog_l2e_module = fx_submodule(NULL, "mog_l2e", "mog_l2e");
|
||||
index_t number_of_gaussians = fx_param_int(mog_l2e_module, "K", 1);
|
||||
fx_format_param(mog_l2e_module, "D", "%d", data_points.n_rows());
|
||||
index_t dimension = fx_param_int_req(mog_l2e_module, "D");;
|
||||
|
||||
////// RUNNING AN OPTIMIZER TO MINIMIZE THE L2 ERROR //////
|
||||
|
||||
datanode *opt_module = fx_submodule(NULL, "opt", "opt");
|
||||
const char *opt_method = fx_param_str(opt_module, "method", "QuasiNewton");
|
||||
index_t param_dim = (number_of_gaussians*(dimension+1)*(dimension+2)/2 - 1);
|
||||
fx_param_int(opt_module, "param_space_dim", param_dim);
|
||||
|
||||
index_t optim_flag = (strcmp(opt_method, "NelderMead") == 0 ? 1 : 0);
|
||||
MoGL2E mog;
|
||||
|
||||
|
||||
if (optim_flag == 1) {
|
||||
|
||||
////// OPTIMIZER USING NELDER MEAD METHOD //////
|
||||
|
||||
NelderMead opt;
|
||||
|
||||
////// Initializing the optimizer //////
|
||||
fx_timer_start(opt_module, "init_opt");
|
||||
opt.Init(MoGL2E::L2ErrorForOpt, data_points, opt_module);
|
||||
fx_timer_stop(opt_module, "init_opt");
|
||||
|
||||
////// Getting starting points for the optimization //////
|
||||
double **pts;
|
||||
pts = (double**)malloc((param_dim+1)*sizeof(double*));
|
||||
for(index_t i = 0; i < param_dim+1; i++) {
|
||||
pts[i] = (double*)malloc(param_dim*sizeof(double));
|
||||
}
|
||||
|
||||
fx_timer_start(opt_module, "get_init_pts");
|
||||
MoGL2E::MultiplePointsGenerator(pts, param_dim+1,
|
||||
data_points, number_of_gaussians);
|
||||
fx_timer_stop(opt_module, "get_init_pts");
|
||||
|
||||
////// The optimization //////
|
||||
|
||||
fx_timer_start(opt_module, "optimizing");
|
||||
opt.Eval(pts);
|
||||
fx_timer_stop(opt_module, "optimizing");
|
||||
|
||||
////// Making model with the optimal parameters //////
|
||||
mog.MakeModel(mog_l2e_module, pts[0]);
|
||||
|
||||
}
|
||||
else {
|
||||
|
||||
////// OPTIMIZER USING QUASI NEWTON METHOD //////
|
||||
|
||||
QuasiNewton opt;
|
||||
|
||||
////// Initializing the optimizer //////
|
||||
fx_timer_start(opt_module, "init_opt");
|
||||
opt.Init(MoGL2E::L2ErrorForOpt, data_points, opt_module);
|
||||
fx_timer_stop(opt_module, "init_opt");
|
||||
|
||||
////// Getting starting point for the optimization //////
|
||||
double *pt;
|
||||
pt = (double*)malloc(param_dim*sizeof(double));
|
||||
|
||||
//index_t rs = 0;
|
||||
//while ( rs < 5) {
|
||||
|
||||
fx_timer_start(opt_module, "get_init_pt");
|
||||
MoGL2E::InitialPointGenerator(pt, data_points, number_of_gaussians);
|
||||
fx_timer_stop(opt_module, "get_init_pt");
|
||||
|
||||
////// The optimization //////
|
||||
|
||||
fx_timer_start(opt_module, "optimizing");
|
||||
opt.Eval(pt);
|
||||
fx_timer_stop(opt_module, "optimizing");
|
||||
|
||||
////// Making model with optimal parameters //////
|
||||
mog.MakeModel(mog_l2e_module, pt);
|
||||
//printf("minimum achieved: %Lf\n", mog.L2Error(data_points));
|
||||
//}
|
||||
|
||||
}
|
||||
|
||||
long double error = mog.L2Error(data_points);
|
||||
NOTIFY("Minimum L2 error achieved: %Lf", error);
|
||||
mog.Display();
|
||||
|
||||
ArrayList<double> results;
|
||||
mog.OutputResults(&results);
|
||||
|
||||
|
||||
////// OUTPUT RESULTS //////
|
||||
|
||||
const char *output_filename = fx_param_str(NULL, "output", "output.csv");
|
||||
|
||||
FILE *output_file = fopen(output_filename, "w");
|
||||
|
||||
ot::Print(results, output_file);
|
||||
fclose(output_file);
|
||||
fx_done();
|
||||
|
||||
return 1;
|
||||
}
|
||||
@@ -1,156 +0,0 @@
|
||||
/**
|
||||
* @author Parikshit Ram (pram@cc.gatech.edu)
|
||||
* @file phi.h
|
||||
*
|
||||
* This file computes the Gaussian probability
|
||||
* density function
|
||||
*/
|
||||
#include "fastlib/fastlib.h"
|
||||
#include "fastlib/fastlib_int.h"
|
||||
#include <cmath>
|
||||
|
||||
/**
|
||||
* Calculates the multivariate Gaussian probability density function
|
||||
*
|
||||
* Example use:
|
||||
* @code
|
||||
* Vector x, mean;
|
||||
* Matrix cov;
|
||||
* ....
|
||||
* long double f = phi(x, mean, cov);
|
||||
* @endcode
|
||||
*/
|
||||
|
||||
long double phi(Vector& x , Vector& mean , Matrix& cov) {
|
||||
|
||||
long double det, f;
|
||||
double exponent;
|
||||
index_t dim;
|
||||
Matrix inv;
|
||||
Vector diff, tmp;
|
||||
|
||||
dim = x.length();
|
||||
la::InverseInit(cov, &inv);
|
||||
det = la::Determinant(cov);
|
||||
|
||||
if( det < 0){
|
||||
det = -det;
|
||||
}
|
||||
la::SubInit(mean,x,&diff);
|
||||
la::MulInit(inv, diff, &tmp);
|
||||
exponent = la::Dot(diff, tmp);
|
||||
long double tmp1, tmp2, tmp3;
|
||||
tmp1 = 1;
|
||||
tmp2 = dim;
|
||||
tmp2 = tmp2/2;
|
||||
tmp2 = pow((2*(math::PI)),tmp2);
|
||||
tmp1 = tmp1/tmp2;
|
||||
tmp3 = 1;
|
||||
tmp2 = sqrt(det);
|
||||
tmp3 = tmp3/tmp2;
|
||||
tmp2 = -exponent;
|
||||
tmp2 = tmp2 / 2;
|
||||
f = (tmp1*tmp3*exp(tmp2));
|
||||
|
||||
return f;
|
||||
}
|
||||
|
||||
/**
|
||||
* Calculates the univariate Gaussian probability density function
|
||||
*
|
||||
* Example use:
|
||||
* @code
|
||||
* double x, mean, var;
|
||||
* ....
|
||||
* long double f = phi(x, mean, var);
|
||||
* @endcode
|
||||
*/
|
||||
|
||||
long double phi(double x, double mean, double var) {
|
||||
|
||||
long double f;
|
||||
|
||||
f = exp(-1.0*((x-mean)*(x-mean)/(2*var)))/sqrt(2*math::PI*var);
|
||||
return f;
|
||||
}
|
||||
|
||||
/**
|
||||
* Calculates the multivariate Gaussian probability density function
|
||||
* and also the gradients with respect to the mean and the variance
|
||||
*
|
||||
* Example use:
|
||||
* @code
|
||||
* Vector x, mean, g_mean, g_cov;
|
||||
* ArrayList<Matrix> d_cov; // the dSigma
|
||||
* ....
|
||||
* long double f = phi(x, mean, cov, d_cov, &g_mean, &g_cov);
|
||||
* @endcode
|
||||
*/
|
||||
|
||||
long double phi(Vector& x, Vector& mean, Matrix& cov, ArrayList<Matrix>& d_cov, Vector *g_mean, Vector *g_cov){
|
||||
|
||||
long double det, f;
|
||||
double exponent;
|
||||
index_t dim;
|
||||
Matrix inv;
|
||||
Vector diff, tmp;
|
||||
|
||||
dim = x.length();
|
||||
la::InverseInit(cov, &inv);
|
||||
det = la::Determinant(cov);
|
||||
|
||||
if( det < 0){
|
||||
det = -det;
|
||||
}
|
||||
la::SubInit(mean,x,&diff);
|
||||
la::MulInit(inv, diff, &tmp);
|
||||
exponent = la::Dot(diff, tmp);
|
||||
long double tmp1, tmp2, tmp3;
|
||||
tmp1 = 1;
|
||||
tmp2 = dim;
|
||||
tmp2 = tmp2/2;
|
||||
tmp2 = pow((2*(math::PI)),tmp2);
|
||||
tmp1 = tmp1/tmp2;
|
||||
tmp3 = 1;
|
||||
tmp2 = sqrt(det);
|
||||
tmp3 = tmp3/tmp2;
|
||||
tmp2 = -exponent;
|
||||
tmp2 = tmp2 / 2;
|
||||
f = (tmp1*tmp3*exp(tmp2));
|
||||
|
||||
// Calculating the g_mean values which would be a (1 X dim) vector
|
||||
la::ScaleInit(f,tmp,g_mean);
|
||||
|
||||
// Calculating the g_cov values which would be a (1 X (dim*(dim+1)/2)) vector
|
||||
double *g_cov_tmp;
|
||||
g_cov_tmp = (double*)malloc(d_cov.size()*sizeof(double));
|
||||
for(index_t i = 0; i < d_cov.size(); i++){
|
||||
Vector tmp_d;
|
||||
// Matrix inv_d;
|
||||
Matrix inv_d, tmp_mat_1, tmp_mat_2;
|
||||
double tmp_d_cov_d_r;
|
||||
|
||||
la::MulInit(d_cov[i],inv,&tmp_mat_1);
|
||||
la::MulInit(inv, tmp_mat_1, &tmp_mat_2);
|
||||
la::MulInit(tmp_mat_2, diff, &tmp_d);
|
||||
tmp_d_cov_d_r = la::Dot(diff, tmp_d);
|
||||
|
||||
// la::MulInit(d_cov[i], tmp, &tmp_d);
|
||||
// tmp_d_cov_d_r = la::Dot(tmp_d,tmp);
|
||||
|
||||
la::MulInit(inv,d_cov[i],&inv_d);
|
||||
|
||||
double trace = 0;
|
||||
for(index_t j = 0; j < dim; j++) {
|
||||
trace += inv_d.get(j,j);
|
||||
}
|
||||
|
||||
tmp_d_cov_d_r -= trace;
|
||||
//printf("trace = %lf\n", trace);
|
||||
|
||||
g_cov_tmp[i] = f*tmp_d_cov_d_r/2;
|
||||
}
|
||||
g_cov->Copy(g_cov_tmp,d_cov.size());
|
||||
|
||||
return f;
|
||||
}
|
||||
Reference in New Issue
Block a user