FastMDS.h

Go to the documentation of this file.

/******************************************************************************
 *
 *  Copyright (C) 2008 Max Hermann
 *
 *  This file is part of the "CEDRIC Event Display" application.
 *
 *****************************************************************************/

#ifndef REDUKT_FastMDS_H
#define REDUKT_FastMDS_H

#include "CMDS.h"
#include "SVD.h"

#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <ext/algorithm> // for random_sample_n()
#include <limits>
#include <iostream>
#include <assert.h>

template<class ValueType>
class FastMDS
{
public:
        typedef ublas::matrix<ValueType> Matrix;
        typedef ublas::vector<ValueType> Vector;

        FastMDS( Matrix& D, int k, int p, int s );

        const Matrix& getPoints() const { return m_points; }

private:
        Matrix m_points;
};


template<class T>
FastMDS<T>::FastMDS( Matrix& D, int k, int p, int s ): m_points(D.size1(),k)
{
        using namespace ublas;
        using namespace ublasTools;
        using std::cout;
        using std::endl;
        using std::flush;
        
        assert( D.size1() == D.size2() );

        int n = D.size1();
        //int p = ;
        //int s = 2*(k+1);
        
        cout << "FastMDS with n=" << n << ", p=" << p << ", dim_low=k=" << k << ", s=" << s << endl;
        
        // compute s*p sample indices
        cout << "Computing sample indices, taking " << s << " samples from each range, total " << s*p << endl;
        std::vector<int> indices(n); for( int i=0; i < n; ++i ) indices[i]=i;
        std::vector<int> sample( s*p );
        std::vector<int>::iterator it = sample.begin();
        for( int i=0; i < p; ++i )
        {
                //cout << "[" << i*(n/p) << "," << (i+1)*(n/p) << "]" << endl;
                
                // take s samples from D_i submatrix range
                it = __gnu_cxx::random_sample_n( indices.begin() +  i   *(n/p), 
                                                 indices.begin() + (i+1)*(n/p),
                                                 it,
                                                                     s   );
        }
        cout << "Sample indices are: " << endl;
        for( int i=0; i < s*p; ++i )
                cout << sample[i] << " ";
        cout << endl;

        // build align matrix M
        Matrix M( s*p, s*p );
        cout << "Build align " << M.size1() << "x" << M.size2() << " matrix M..." << flush;
        for( int row=0; row < s*p; ++row )
                for( int col=0; col < s*p; ++col )
                {
                        M(row,col) = D(sample[row],sample[col]);
                }
        cout << "done" << endl;

        // solve CMDS for M
        cout << "Solve CMDS for align matrix M..." << flush;
        CMDS<T> M_cmds( M );
        M_cmds.project( k );   // M = CMDS(M) projected to k dimensions
        cout << "done (result is a " << M.size1() << "x" << M.size2() << " matrix)" << endl;
        

        for( int i=0; i < p; ++i )
        {
                // sub matrix D_i from proxy
                int first =  i   *(n/p),
                    last  = (i+1)*(n/p);

                cout << "Submatrix D_" << i << " from point " << first << " to " << last << endl;
                Matrix D_i = matrix_range<Matrix>( D,
                                                   range(first,last),
                                                       range(first,last) );

                // solve CMDS for D_i
                cout << "Solve CMDS for D_" << i << "..." << flush;
                CMDS<T> Di_cmds( D_i );                             
                Di_cmds.project( k );                            // D_i = CMDS(D_i) projected to k dimension
                cout << "done (result is a " << D_i.size1() << "x" << D_i.size2() << " matrix)" << endl;
                
                

                // build matrices containing solutions for sampled points in D_i and M
                cout << "  Build " << k << "x" << s << " matrices mdsDs and mdsMs" << endl;
                Matrix mdsDs( k, s );
                Matrix mdsMs( k, s );
                for( int row=0; row < s; ++row )
                {
                        int ofs = i*s + row;
                        //~ cout << "* sample[" << ofs << "] - i*(n/p) = " << flush;
                        int idx = sample[ofs] - i*(n/p);
                        //~ cout << idx << endl;
                        
                        for( int col=0; col < k; ++col )
                        {
                                mdsDs(col,row) = D_i( idx , col );
                                mdsMs(col,row) = M  ( ofs, col );
                        }
                }
                
                // fit affine transformation  mdsMs = A*mdsDs + b
                cout << "  Find affine mapping A_" << i << endl;
                Matrix A; Vector b;
                if( !affine_fit<T>( mdsDs, mdsMs, A, b ) )
                {
                        cout << "BOGUS!" << endl;
                        assert( false );
                }
                
        #if 0
                // map mdsDs into mdsDs for debugging purposes
                cout << "  mdsDs = " << mdsDs << endl;
                cout << "  mdsMs = " << mdsMs << endl;
                Matrix estimate = prod( A, mdsDs );
                for( unsigned int j=0; j < estimate.size2(); ++j )
                        column(estimate,j) += b;
                cout << "  A*mdsDs+b1' = " << estimate << endl;
        #endif
                
        
                

                // affine mapping A_i is stored in A
                
                cout << "Align CMDS(D_" << i << ") with CMDS(M)..." << flush;
                
                // align CMDS(D_i) in common coordinate system
                D_i = trans( D_i );
                D_i = prod( A, D_i );
                cout << "done (result is a " << D_i.size1() << "x" << D_i.size2() << " matrix)" << endl;
                
                cout << "Copy resulting points into solution matrix" << endl;
                D_i = trans(D_i);
                
                // copy points to result matrix
                for( unsigned int row=0; row < D_i.size1(); ++row )
                        for( int col=0; col < k; ++col )  // k must match D_i.size2() (or something is fucked up)
                                m_points(row+first,col) = D_i(row,col);
        }
}