mog_l2e.h

Go to the documentation of this file.

/* MLPACK 0.2
 *
 * Copyright (c) 2008, 2009 Alexander Gray,
 *                          Garry Boyer,
 *                          Ryan Riegel,
 *                          Nikolaos Vasiloglou,
 *                          Dongryeol Lee,
 *                          Chip Mappus, 
 *                          Nishant Mehta,
 *                          Hua Ouyang,
 *                          Parikshit Ram,
 *                          Long Tran,
 *                          Wee Chin Wong
 *
 * Copyright (c) 2008, 2009 Georgia Institute of Technology
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License as
 * published by the Free Software Foundation; either version 2 of the
 * License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
 * 02110-1301, USA.
 */
#ifndef MOG_L2E_H
#define MOG_L2E_H

#include <fastlib/fastlib.h>

const fx_entry_doc mog_l2e_entries[] = {
  {"K", FX_PARAM, FX_INT, NULL,
   " The number of Gaussians in the mixture model."
   " (defaults to 1)\n"},
  {"D", FX_RESERVED, FX_INT, NULL,
   " The number of dimensions of the data on which the"
   " the mixture model is to be fit.\n"},
  FX_ENTRY_DOC_DONE
};

const fx_module_doc mog_l2e_doc = {
  mog_l2e_entries, NULL,
  " This program defines a Gaussian mixture model"
  " and calculates the L2 error for the present"
  " parameter setting to be used by an optimizer.\n"
};

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);
    }
  }

  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);
  }

  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);
        }
      }
    }
  }                     
  
  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];
  }


  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;
  }

  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);
        }
      }
    }
  }

  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");
  }

  
  long double L2Error(const Matrix&, Vector* = NULL);
  
  long double RegularizationTerm_(Vector* = NULL);
   
  long double GoodnessOfFitTerm_(const Matrix&, Vector* = NULL);

  static void MultiplePointsGenerator(double**, index_t, 
                                      const Matrix&, index_t);

  static void InitialPointGenerator(double*, const Matrix&, index_t);

  static void KMeans_(const Matrix&, ArrayList<Vector>*,
                      ArrayList<Matrix>*, Vector*, index_t);

  static void min_element(Matrix&, index_t*);

  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);
    }
  }

  static long double L2ErrorForOpt(Vector& params, const Matrix& data) {
    return L2ErrorForOpt(params, data, NULL);
  }

};

#endif