LDA.cpp

/*
This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
  Multimediale Systeme der Informationsverarbeitung
  Institut fuer Informatik
  Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS 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 Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include <MathAlgo/LDA.hh>
#include <MathAlgo/PCA.hh>
#include <MathAlgo/SVD.hh>
#include <Base/Math/Random.hh>

#include <Base/Common/BIASpragma.hh>
#include <Base/Common/FileHandling.hh>

using namespace BIAS;
using namespace std;


#define REDUCTION_ZERO_THRESH 1E-8

template<typename T>
void outveclda(BIAS::Vector<T> vec){
  for (int i=0;i<vec.size();i++){
    cout << vec[i] << " ";
  }
  cout << endl;
}

//only for debugging:
void LDA::AnalyzeData(const std::vector<std::vector<BIAS::Vector<LDAType> > > &vec){

  BIASASSERT(!vec.empty());
  BIASASSERT(!vec[0].empty());
  BIAS::Vector<LDAType> classMean;

  int dimensions=vec[0][0].size();

  Matrix<double> classCov(dimensions, dimensions,0.0),
                 cov(dimensions, dimensions, 0.0);

  classMean.newsize(dimensions);

  double error=0,error2=0, error3=0, maxerror=0, maxerror2=0;

  int counter=0, counter2=0;

  ComputeMeans(vec);

  for (int classNumber=0;classNumber<(int)vec.size();classNumber++){
    cov.SetZero();
    // save offset classNumber
    error=0;
    const std::vector<BIAS::Vector<LDAType> > & class_i = vec[classNumber];
    for (int classVecNumber=0; classVecNumber<(int)class_i.size();
         classVecNumber++){
      // save offset l
      const BIAS::Vector<LDAType> & classVector = class_i[classVecNumber];
      //outveclda(classVector);
      for (int dim=0; dim<dimensions; dim++){
        // compute kli once for all j
        const LDAType classVectordim =  classVector[dim];
        const LDAType classMean = classMeans_[classNumber][dim];
        for (int dim2=0;dim2<=dim;dim2++){
          cov[dim][dim2] += (classVectordim-classMean)*(classVector[dim2]-classMeans_[classNumber][dim2]);
        }
        error3=fabs(classVectordim-classMean);
        if (error3>maxerror) maxerror=error3;
        error+=error3;

      }
      counter2++;
    }
    cov/=double(class_i.size());
    for (int i=0; i<dimensions; i++){
      for (int j=0; j<=i; j++){
        classCov[i][j] += cov[i][j];
      }
    }
    error=error/(double)(128*vec[classNumber].size());
    if (error>maxerror2) maxerror2=error;
    counter++;
    error2+=error;

  }

  for (int i=0; i<dimensions; i++){
    for (int j=0; j<=i; j++){
      classCov[j][i] = classCov[i][j];
    }
  }
  classCov/=double(vec.size());


  SVD svd(false); // no absolute zero threshold
  Matrix<double> U;
  Vector<double> S;

  //get sure that matrix as full rank
  svd.Compute(classCov, REDUCTION_ZERO_THRESH);
  S = svd.GetS();

  error2/=(double)counter;

  cout << "Trailset has Trails = " << vec.size() << endl;
  cout << "Mean-Error = " << error2 << endl;
  cout << "Max-Single-Error = " << maxerror << endl;
  cout << "Max-Trail-Error = " << maxerror2 << endl;

  cout << "Eigenvalues from S_w = ";
  outveclda(S);

}

void LDA::GenerateRandomTestData(int vectorSize, int numberOfClasses,
  int classSizeMin, int classSizeMax, int vectorEntryMin, int vectorEntryMax,
  double stdDeviation, double mainDirectionStdDeviation,
  std::vector<std::vector<BIAS::Vector<LDAType> > > &vec,
  std::vector<BIAS::Vector<LDAType> > &means, bool Normalize){

  int classSize;
  double l;
  vec.resize(numberOfClasses);
  means.resize(numberOfClasses);
  SVD svd(false);// no absolute zero thresh
  Matrix<double> U(vectorSize,vectorSize,0.0),S(vectorSize,vectorSize,0.0),
                 UMain(vectorSize,vectorSize,0.0), SMain(vectorSize,vectorSize,0.0);
  Vector<double> T(vectorSize,0.0), SVec(vectorSize,0.0);
  Random random;
#ifndef BIAS_DEBUG
  random.Reset();
#endif

  //generate Eigenvalues (in matrix SMain) of covariance-matrix for class-means
  SMain.SetZero();
  for (int k=0;k<vectorSize;k++){
    T[k]=random.GetUniformDistributed(-1,1);
    SMain[k][k]=fabs(random.GetNormalDistributed(0,mainDirectionStdDeviation*mainDirectionStdDeviation));
  }
  //cout << "Haupteigenwerte=" << SMain << endl;

  //generate eigenvectors (matrix UMain) of covariance-matrix for class-means
  //we spend power to compute more svds so there is more random-factor
  UMain.SetZero();
  for (int k=0;k<vectorSize;k++){
    T/=T.NormL2();
    UMain.SetRow(k, T);
    if (k==vectorSize-1) break;
    svd.Compute(UMain);
    for (int j=k+1;j<vectorSize;j++){
      T+=svd.GetVT().GetRow(j)*random.GetUniformDistributed(-1,1);
    }
  }
/*
  UMain.SetZero();
  T/=T.NormL2();
  UMain.SetRow(0, T);
  svd.Compute(UMain);
  */
  //compute covariance-matrix (now UMain) for class-means
  UMain=svd.GetU()*SMain*svd.GetVT();
  //cout << "HauptKovarianzmatrix=" << UMain << endl;

  //for each class do:
  for (int i=0;i<numberOfClasses;i++){
    #ifdef BIAS_DEBUG
    if ((GetDebugLevel()&1)!=0)
      cout << i << " " << flush;
    #endif
    means[i].newsize(vectorSize);
    S.SetZero();
    //class means are generated, and eigenvalues (S) for class-covariance-matrices
    for (int k=0;k<vectorSize;k++){
      //Mean:
      means[i][k] = LDAType( random.GetUniformDistributedInt(vectorEntryMin,vectorEntryMax) );
      //Eigenvalues:
      S[k][k]=fabs(random.GetNormalDistributed(0,stdDeviation*stdDeviation));
      //init first eigenvector
      T[k]=random.GetUniformDistributed(-1,1);
    }


    //Eigenvectors:
    /*
    U.SetZero();
    for (int k=0;k<vectorSize;k++){
      T/=T.NormL2();
      U.SetRow(k, T);
      if (k==vectorSize-1) break;
      svd.Compute(U);
      VT=svd.GetVT();
      for (int j=k+1;j<vectorSize;j++){
        T+=VT.GetRow(j);
      }
    }
    */

    //cout << "Klassen-Eigenwerte=" << S << endl;
    //generate eigenvectors for each class
    U.SetZero();
    T/=T.NormL2();
    U.SetRow(0, T);
    svd.Compute(U);
    U=svd.GetU()*S*svd.GetVT();
    //cout << "Klassen-Kovarianzmatrix=" << U << endl;
    U+=UMain;
    //cout << "Addierte Kovarianzmatrix=" << U << endl;
    U/=2.0;
    svd.Compute(U);
    SVec=svd.GetS();
    //cout << "Neue Eigenwerte=" << endl;outveclda(SVec);
    U=svd.GetVT();



    //generate class-size
    classSize=random.GetUniformDistributedInt(classSizeMin,classSizeMax);
    vec[i].resize(classSize);
    //cout << "S" << S << "class " << i << " mean " << means[i] << endl ;
    for (int j=0;j<classSize;j++){

      vec[i][j]=means[i];
      Vector<LDAType> *vec_ij=&vec[i][j];
      //for each dimension get vectors position in class
      for (int k=0;k<vectorSize;k++){
        //get distance to class mean:
        l=random.GetNormalDistributed(0, sqrt(SVec[k]));

        //cout << "class " << i << " vec " << j << " add l=" << l << " * Evec(" << k << ")" << endl;
        //multiply distance with eigenvector:
        for (int n=0; n<int(vectorSize); n++){
          (*vec_ij)[n] += LDAType(l*U[k][n]);
        }
      }
      for (int n=0; n<int(vectorSize); n++){
        (*vec_ij)[n]=min(LDAType(vectorEntryMax),max(LDAType(vectorEntryMin),(*vec_ij)[n]));
      }
      if (Normalize) {
        (*vec_ij) /= (*vec_ij).NormL2()*255.0;
      }
      //cout << "class " << i << " vec " << j << " " <<  vec[i][j] << endl;
    }
  }
  #ifdef BIAS_DEBUG
  if ((GetDebugLevel()&1)!=0)
    cout << endl;
  #endif

}

void LDA::
ComputeReductionMatrix(const vector<vector<BIAS::Vector<LDAType> > > &vec,
                       BIAS::Matrix<LDAType> &matrix, int reductionSize,
                       std::vector<BIAS::Matrix<LDAType> > *covs){

  AddDebugLevel(1);
  //init
  if (reductionSize==-1)
    reductionSize = reductionSize_;
  BIASASSERT(reductionSize>0);

  BIASASSERT(!vec.empty());
  BIASASSERT(!vec[0].empty());
  const int dimensions = vec[0][0].size();
  if (covs!=NULL) {
    covs->resize(vec.size(), Matrix<float>(dimensions, dimensions,MatrixZero));
  }

  // compute means

  BIASDOUT(1, "compute means "<<FileHandling::CurrentDateTimeString());

  ComputeMeans(vec);

  //covariance-matrix
  Matrix<LDAType> withinclassCov(dimensions, dimensions,0.0),
    interclassCov(dimensions, dimensions, 0.0);

  vector<double> classWeights;
  classWeights.reserve(vec.size());

  // for each class compute its covariance matrix and sum them up weighted in classCov
  BIASDOUT(1, "compute within covariances"<<FileHandling::CurrentDateTimeString() );
  for (int classNumber=0;classNumber<(int)vec.size();classNumber++){
    Matrix<LDAType> &cov = (covs==NULL)? interclassCov : (*covs)[classNumber];
    cov.SetZero();
    const std::vector<BIAS::Vector<LDAType> > & class_i = vec[classNumber];
    if (equalClassWeights_) classWeights.push_back(1.0);
    else classWeights.push_back(class_i.size());
    // cov approximation only possible when more than one measurement available
    if (class_i.size()>1) {
      // save offset classNumber
      for (int classVecNumber=0; classVecNumber<(int)class_i.size(); classVecNumber++){
        // save offset l
        const BIAS::Vector<LDAType> & classVector = class_i[classVecNumber];
        for (int dim=0; dim<dimensions; dim++){
          // compute kli once for all j
          const LDAType classVectordim =  classVector[dim];
          const LDAType classMean = classMeans_[classNumber][dim];
          const LDAType diff  = (classVectordim-classMean);
          for (int dim2=0;dim2<=dim;dim2++){
            cov[dim][dim2] += diff *
              (classVector[dim2]-classMeans_[classNumber][dim2]);
          }
        }
      }
      // all classes with equal probability ?
      if (equalClassWeights_)
        cov *= LDAType(1.0/LDAType(class_i.size()-1));

      // just one copy of other triangular matrix at the end
      for (int i=0; i<dimensions; i++){
        for (int j=0; j<=i; j++){
          withinclassCov[i][j] += cov[i][j];
        }
      }
    }
  }
  // copy symmetic part
  for (int i=0; i<dimensions; i++){
    for (int j=0; j<i; j++){
      withinclassCov[j][i] = withinclassCov[i][j];
    }
  }

  BIASDOUT(2, "within class cov = "<<withinclassCov);

  // compute second central moment of class means
  BIASDOUT(1, "compute inter-class covariances "<<FileHandling::CurrentDateTimeString());
  interclassCov.SetZero();
  for (int classNumber=0;classNumber<(int)classMeans_.size();classNumber++){
    Vector<LDAType> & classMean = classMeans_[classNumber];
    const  LDAType w = classWeights[classNumber];
    for (int i=0;i<classMean.size();i++){
      const LDAType weight = w * (classMean[i] - mean_[i]);
      for (int j=0; j<=i; j++){
        // weight classes ?
        interclassCov[i][j] += weight * (classMean[j]-mean_[j]);
      }
    }
  }

  for (int i=0; i<dimensions; i++){
    for (int j=0;j<i;j++) {
      interclassCov[j][i] = interclassCov[i][j];
    }
  }


  // within and inter class cov are now scaled by number of classes



  BIASDOUT(2, "inter class cov = "<<interclassCov);
  ComputeReductionMatrix(withinclassCov, interclassCov, matrix, reductionSize);

}




void LDA::ComputeReductionMatrix(Matrix<LDAType>& withinclassCov,
                                 Matrix<LDAType>& interclassCov,
                                 BIAS::Matrix<LDAType> &matrix,
                                 int reductionSize) {

  const int dimensions = interclassCov.num_rows();

  // make symmetric double matrix, enforce positive diagonal and symmetry
  Matrix<double> interclassCovd(dimensions, dimensions);
  for (int i=0; i<dimensions; i++){
    interclassCovd[i][i] = 2.0 * fabs(interclassCov[i][i]);
    for (int j=0;j<i;j++) {
      interclassCovd[j][i] = interclassCovd[i][j] = double(interclassCov[i][j])
        + double(interclassCov[j][i]) ;
    }
  }
  BIASWARNONCE("Copying from float to double ...");


  // make SVD and decide if to go on with PCA
  SVD svd(false);// no absolute zero threshold !
  Matrix<double> U, withinclassCovd(dimensions, dimensions);
  for (int i=0; i<dimensions; i++){
    withinclassCovd[i][i] = 2.0*fabs(withinclassCov[i][i]);
    for (int j=0;j<i;j++) {
      withinclassCovd[i][j] = withinclassCovd[j][i] =
        double(withinclassCov[i][j])+double(withinclassCov[j][i]);
    }
  }

  BIASDOUT(1, "classcov svd."<<FileHandling::CurrentDateTimeString());
  int rank = 0;
  // make sure that matrix has full rank
  if (svd.Compute(withinclassCovd, REDUCTION_ZERO_THRESH)==0) {
    rank=svd.Rank();
  } else {
    BIASERR("SVD failed");
  }
  BIASDOUT(1, "rank=" << svd.Rank() << ", dimensions=" << dimensions);

  if (rank==0){
    // make PCA on means if rank=0
    BIASDOUT(1, " rank=0, make PCA");
    svd.Compute(interclassCovd);
    U = svd.GetVT();
  } else {
    // try regularization if rank is not full:
    if (rank<dimensions){
      double ev = svd.GetS()[rank-1];
      BIASDOUT(1, "Singularity, s="<< svd.GetS()<<" rank="<< svd.Rank()
               <<". Regularizing and trying again."<<FileHandling::CurrentDateTimeString());
      for (int i=0;i<dimensions;i++){
        withinclassCovd[i][i] += ev;
      }
      svd.Compute(withinclassCovd, REDUCTION_ZERO_THRESH);

      BIASDOUT(1, "SECOND TRY yielded rank="
               << svd.Rank() << ", added " << ev);

      if ((int)svd.Rank()<dimensions) {
        BIASERR("REGULARIZATION FAILED - ABORTING");
        BIASABORT;
      }
    }
    //cout << "InnerClassCov=" << withinclassCov << "BetweenClassCov=" << interclassCov ;
    //cout << "U=" << svd.GetU() << "V=" << svd.GetV() << "S=" <<  svd.GetS();

    // last steps for LDA-matrix, see paper
    // from Max Welling: "Fisher Linear Discriminant Analysis"
    BIASDOUT(1, "Inverting within class cov ... "<<FileHandling::CurrentDateTimeString());
    withinclassCovd = svd.Invert();

    //cout << "INV(InnerClassCov)=" << withinclassCov ;
    BIASDOUT(1, "Factorizing inter class cov  ... "<<FileHandling::CurrentDateTimeString());
    svd.Compute(interclassCovd);
    rank = svd.Rank();
    if (rank<dimensions){
      BIASDOUT(1, "inter class cov is singular, regularizing and trying "
               "again ... "<<FileHandling::CurrentDateTimeString());
      for (int i=0;i<dimensions;i++){
        interclassCovd[i][i] += svd.GetS()[rank-1];
      }
      svd.Compute(interclassCovd);
    }
    BIASDOUT(1, "Computing inter class root "<<FileHandling::CurrentDateTimeString());
    //cout << "U=" << svd.GetU() << "V=" << svd.GetV() << "S=" <<  svd.GetS();
    interclassCovd=svd.SqrtT();

    //cout <<  "SqrtT(BetweenClassCov)=" << interclassCovd ;
    BIASDOUT(1, "composing and solving transformation matrix "<<FileHandling::CurrentDateTimeString());
    svd.Compute(interclassCovd * withinclassCovd * interclassCovd);

    //cout <<  "cov*withinclassCov*cov=" << interclassCovd*withinclassCovd*interclassCovd ;

    U = svd.GetV();

    //cout << "U=" << svd.GetU() << "V=" << svd.GetV() << "S=" <<  svd.GetS();
    BIASDOUT(1, "inverting covd "<<FileHandling::CurrentDateTimeString());
    interclassCovd = svd.Invert(interclassCovd);

    //cout <<  "INV(SqrtT(BetweenClassCov))=" << cov ;

    U = interclassCovd * U;
    //cout << "U=" << U ;

    /*
    svd.Compute(withinclassCov*cov);
    U=svd.GetV();
    */
  }

  BIASDOUT(1, "factorizing U "<<FileHandling::CurrentDateTimeString());
  svd.Compute(U);
  U /= svd.GetS()[int(dimensions/2)];

  matrix.newsize(dimensions, reductionSize);
  for (int j=0;j<matrix.num_rows();j++){
    for (int i=0;i<matrix.num_cols();i++){
      matrix[j][i] = LDAType( U[j][i] );
    }
  }
}

void LDA::ComputeMeans(const vector<vector<BIAS::Vector<LDAType> > > &vec){
  const int dimensions = vec[0][0].size();

  classMeans_.resize(vec.size());
  mean_.newsize(dimensions);
  mean_.SetZero();

  allCount_=0;
  for (int i=0;i<(int)vec.size();i++){
    classMeans_[i].newsize(dimensions);
    classMeans_[i].SetZero();
    allCount_+=vec[i].size();
    for (int k=0;k<(int)vec[i].size();k++){
      classMeans_[i]+=vec[i][k];
    }
    mean_+=classMeans_[i];
    classMeans_[i] /= LDAType(vec[i].size());
  }
  // this is the absaolute mean of all points, not the mean of the class means
  mean_ /= LDAType(allCount_);

}

void LDA::GetMean(BIAS::Vector<LDAType> &mean){
  mean=mean_;
}

void LDA::GetClassMeans(std::vector<BIAS::Vector<LDAType> > &means){
  means=classMeans_;
}

void LDA::SetReductionSize(int size){
  reductionSize_=size;
}


void  LDA::ComputeWithinAndInterClassCovs(const vector<BIAS::Vector<LDAType> >&
                                          classmeans,
                                          const vector<BIAS::Matrix<LDAType> >&
                                          covs,
                                          Matrix<LDAType>& withinClassCov,
                                          Matrix<LDAType>& interClassCov){
  BIASASSERT(classmeans.size()==covs.size());
  BIASASSERT(withinClassCov.num_rows()==withinClassCov.num_cols());
  BIASASSERT(interClassCov.num_rows()==interClassCov.num_cols());
  BIASASSERT(withinClassCov.num_rows()==interClassCov.num_cols());
  BIASASSERT(classmeans.size()==covs.size());
  BIASASSERT(withinClassCov.num_rows()==covs[0].num_rows());
  Vector<float> mean(classmeans[0].dim(), 0.0f);
  PCA pca;
  pca.ComputeMean(classmeans, mean);
  pca.ComputeScatter(classmeans, mean, interClassCov);
  for (unsigned int i=0; i<classmeans.size(); i++) {
    withinClassCov += covs[i];
  }
  interClassCov *= 1.0f / float(interClassCov.NormFrobenius());
  withinClassCov*= 1.0f / float(withinClassCov.NormFrobenius());
}


void LDA::ComputeAnonymousReduction(const std::vector<std::vector<
                                    BIAS::Vector<LDAType> > > &vec,
                                    BIAS::Matrix<LDAType> &matrix,
                                    int reductionSize,
                                    const std::vector< std::vector<BIAS::
                                    Matrix<LDAType> > >*covs) {
  BIASASSERT(!vec.empty());
  BIASASSERT(!vec[0].empty());
  const int dim = vec[0][0].Size();
  Matrix<LDAType> WithinClassCov(dim, dim, MatrixZero),
    InterClassCov(dim, dim, MatrixZero),
    scatter(dim, dim, MatrixZero), interscatter(dim, dim, MatrixZero);

  PCA myPCA;
  for (unsigned int i=0; i< vec.size(); i++) {
    // run over all images
    Vector<LDAType> mean;
    myPCA.ComputeMean(vec[i], mean);
    myPCA.ComputeScatter(vec[i], mean, interscatter);
    if (covs!=NULL) {
      scatter.SetZero();
      for (unsigned int c=0; c<(*covs)[i].size(); c++) {
        scatter += (*covs)[i][c];
      }
      //scatter.Normalize();
    }
    InterClassCov += interscatter;
    WithinClassCov += scatter;
  }
  ComputeReductionMatrix(WithinClassCov, InterClassCov, matrix, reductionSize);
}