LCOV - code coverage report
Current view: top level - routines - fa.hpp (source / functions) Hit Total Coverage
Test: clean.info Lines: 25 25 100.0 %
Date: 2013-05-24 Functions: 1 4 25.0 %
Branches: 102 594 17.2 %

           Branch data     Line data    Source code
       1                 :            : /* This software is distributed under BSD 3-clause license (see LICENSE file).
       2                 :            :  *
       3                 :            :  * Copyright (c) 2013 Sergey Lisitsyn, Fernando J. Iglesias Garcia
       4                 :            :  */
       5                 :            : 
       6                 :            : #ifndef TAPKEE_FA_H_
       7                 :            : #define TAPKEE_FA_H_
       8                 :            : 
       9                 :            : namespace tapkee
      10                 :            : {
      11                 :            : namespace tapkee_internal
      12                 :            : {
      13                 :            : 
      14                 :            : template <class RandomAccessIterator, class FeatureVectorCallback>
      15                 :          1 : DenseMatrix project(RandomAccessIterator begin, RandomAccessIterator end, FeatureVectorCallback callback,
      16                 :            :                 IndexType dimension, const IndexType max_iter, const ScalarType epsilon,
      17                 :            :                 const IndexType target_dimension, const DenseVector& mean_vector)
      18                 :            : {
      19 [ +  - ][ +  - ]:          1 :         timed_context context("Data projection");
         [ +  - ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
      20                 :            : 
      21                 :            :         // The number of data points
      22                 :          1 :         const IndexType n = end-begin;
      23                 :            : 
      24                 :            :         // Dense representation of the data points
      25                 :            : 
      26   [ +  -  #  #  :          1 :         DenseVector current_vector(dimension);
                   #  # ]
      27                 :            : 
      28 [ +  - ][ +  - ]:          1 :         DenseMatrix X = DenseMatrix::Zero(dimension,n);
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      29                 :            : 
      30 [ +  + ][ #  # ]:         51 :         for (RandomAccessIterator iter=begin; iter!=end; ++iter)
                 [ #  # ]
      31                 :            :         {
      32 [ +  - ][ #  # ]:         50 :                 callback.vector(*iter,current_vector);
                 [ #  # ]
      33 [ +  - ][ +  - ]:         50 :                 X.col(iter-begin) = current_vector - mean_vector;
         [ +  - ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
      34                 :            :         }
      35                 :            : 
      36                 :            :         // Initialize FA model
      37                 :            : 
      38                 :            :         // Initial variances
      39 [ +  - ][ +  - ]:          1 :         DenseMatrix sig = DenseMatrix::Identity(dimension,dimension);
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      40                 :            :         // Initial linear mapping
      41 [ +  - ][ +  - ]:          1 :         DenseMatrix A = DenseMatrix::Random(dimension, target_dimension).cwiseAbs();
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      42                 :            : 
      43                 :            :         // Main loop
      44                 :          1 :         IndexType iter = 0;
      45 [ +  - ][ +  - ]:          1 :         DenseMatrix invC,M,SC;
         [ +  - ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
      46                 :          1 :         ScalarType ll = 0, newll = 0;
      47 [ +  - ][ #  # ]:         56 :         while (iter < max_iter)
                 [ #  # ]
      48                 :            :         {
      49                 :         56 :                 ++iter;
      50                 :            : 
      51                 :            :                 // Perform E-step
      52                 :            : 
      53                 :            :                 // Compute the inverse of the covariance matrix
      54 [ +  - ][ +  - ]:         56 :                 invC = (A*A.transpose() + sig).inverse();
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      55 [ +  - ][ +  - ]:         56 :                 M = A.transpose()*invC*X;
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      56 [ +  - ][ +  - ]:         56 :                 SC = n*(DenseMatrix::Identity(target_dimension,target_dimension) - A.transpose()*invC*A) + M*M.transpose();
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      57                 :            : 
      58                 :            :                 // Perform M-step
      59 [ +  - ][ +  - ]:         56 :                 A = (X*M.transpose())*SC.inverse();
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
      60 [ +  - ][ +  - ]:         56 :                 sig = DenseMatrix(((X*X.transpose() - A*M*X.transpose()).diagonal()/n).asDiagonal()).array() + epsilon;
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
      61                 :            : 
      62                 :            :                 // Compute log-likelihood of FA model
      63 [ +  - ][ +  - ]:         56 :                 newll = 0.5*(log(invC.determinant()) - (invC*X).cwiseProduct(X).sum()/n);
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      64                 :            : 
      65                 :            :                 // Check for convergence
      66 [ +  + ][ +  + ]:         56 :                 if ((iter > 1) && (fabs(newll - ll) < epsilon))
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      67                 :          1 :                         break;
      68                 :            : 
      69                 :         55 :                 ll = newll;
      70                 :            :         }
      71                 :            : 
      72 [ +  - ][ +  - ]:          1 :         return X.transpose()*A;
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      73                 :            : }
      74                 :            : 
      75                 :            : } /* namespace tapkee */
      76                 :            : } /* namespace tapkee_internal */
      77                 :            : 
      78                 :            : #endif /* TAPKEE_FA_H_ */

Generated by: LCOV version 1.9