mixtureDST.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 FASTLIB_MIXTURE_GAUSSIAN_H
#define FASTLIB_MIXTURE_GAUSSIAN_H
#include "fastlib/fastlib.h"
#include "support.h"
class MixtureGauss {
 private:
  ArrayList<Vector> means;

  ArrayList<Matrix> covs;

  Vector prior;

  ArrayList<Matrix> inv_covs;

  Vector det_covs;

  ArrayList<Vector> ACC_means;

  ArrayList<Matrix> ACC_covs;

  Vector ACC_prior;

  double total;

  OT_DEF(MixtureGauss) {
    OT_MY_OBJECT(means);
    OT_MY_OBJECT(covs);
    OT_MY_OBJECT(prior);
    OT_MY_OBJECT(inv_covs);
    OT_MY_OBJECT(det_covs);
    OT_MY_OBJECT(ACC_means);
    OT_MY_OBJECT(ACC_covs);
    OT_MY_OBJECT(ACC_prior);
    OT_MY_OBJECT(total);
  }
 public:
  void InitFromFile(const char* mean_fn, const char* covs_fn = NULL, const char* prior_fn = NULL);

  void InitFromProfile(const ArrayList<Matrix>& matlst, int start, int N);

  void Init(int K, int N);

  void Init(int K, const Matrix& data, const ArrayList<int>& labels);

  void print_mixture(const char* s) const;

  void generate(Vector* v) const;

  double getPDF(const Vector& v) const;

  double getPDF(int cluster, const Vector& v) const;

  const Vector& get_prior() const { return prior; }

  const Vector& get_mean(int k) const { return means[k]; }

  const Matrix& get_cov(int k) const { return covs[k]; }

  int n_clusters() const { return means.size(); }
  
  int v_length() const { return means[0].length(); }

  
  void start_accumulate() {
    total = 0;
    for (int i = 0; i < means.size(); i++) {
      ACC_means[i].SetZero();
      ACC_covs[i].SetZero();
      ACC_prior.SetZero();
    }
  }

  void accumulate(const Vector& v) {
    double s = getPDF(v);
    for (int i = 0; i < means.size(); i++) {
      double p = getPDF(i, v) / s;
      ACC_prior[i] += p;
      la::AddExpert(p, v, &ACC_means[i]);
      Vector d;
      la::SubInit(v, means[i], &d);
      Matrix D;
      D.AliasColVector(d);
      la::MulExpert(p, false, D, true, D, 1.0, &ACC_covs[i]);
    }
    total ++;
  }

  void accumulate_cluster(int i, const Vector& v) {
    la::AddTo(v, &ACC_means[i]);
    Matrix V;
    V.AliasColVector(v);
    la::MulExpert(1.0, false, V, true, V, 1.0, &ACC_covs[i]);
    ACC_prior[i]++;
    total++;
  }

  void accumulate(double p, int i, const Vector& v) {
    la::AddExpert(p, v, &ACC_means[i]);
    Matrix V;
    V.AliasColVector(v);
    la::MulExpert(p, false, V, true, V, 1.0, &ACC_covs[i]);
    ACC_prior[i] += p;
    total += p;
  }

  void end_accumulate_cluster() {
    for (int i = 0; i < means.size(); i++) 
      if (ACC_prior[i] != 0) {
        la::ScaleOverwrite(1.0/ACC_prior[i], ACC_means[i], &means[i]);
        Matrix M;
        M.AliasColVector(means[i]);
        la::MulExpert(-1.0, false, M, true, M, 1.0/ACC_prior[i], &ACC_covs[i]);
        covs[i].CopyValues(ACC_covs[i]);
        prior[i] = ACC_prior[i]/total;

        double det = la::Determinant(covs[i]);
        la::InverseOverwrite(covs[i], &inv_covs[i]);
        det_covs[i] = pow(2.0*math::PI, -means[i].length()/2.0) * pow(det, -0.5);
      }
  }

  void end_accumulate() {
    for (int i = 0; i < means.size(); i++) {
      if (ACC_prior[i] != 0) {
        la::ScaleOverwrite(1.0/ACC_prior[i], ACC_means[i], &means[i]);
        Matrix M;
        M.AliasColVector(means[i]);
        la::MulExpert(-1.0, false, M, true, M, 1.0/ACC_prior[i], &ACC_covs[i]);
        covs[i].CopyValues(ACC_covs[i]);
        prior[i] = ACC_prior[i]/total;

        double det = la::Determinant(covs[i]);
        la::InverseOverwrite(covs[i], &inv_covs[i]);
        det_covs[i] = pow(2.0*math::PI, -means[i].length()/2.0) * pow(det, -0.5);
      }
    }
  }
};
#endif