|| || || || || ||

---

MatNorm2Dist.h

Go to the documentation of this file.

#ifndef Impala_Core_Matrix_MatNorm2Dist_h
#define Impala_Core_Matrix_MatNorm2Dist_h

#include "Core/Matrix/MatFunc.h"
#include "Core/Matrix/MatSum.h"
#include "Core/Matrix/MatMul.h"
#include "Core/Array/Abs.h"
#include "Core/Array/Add.h"
#include "Core/Array/Mul.h"
#include "Core/Array/MulVal.h"
#include "Core/Array/Sqrt.h"
#include "Basis/Timer.h"
#ifdef CUDA
#include "Link/Cuda/Cuda.h"
#include "Link/Cuda/MatNorm2Dist.h"
#endif
#ifdef SSE_USED
#include "Core/Matrix/MatNorm2DistSSE.h"
#endif
//#include "Core/Array/WriteRaw.h"

namespace Impala
{
namespace Core
{
namespace Matrix
{

/*
def distance(a,b):
    """Computes Euclidean distance matrix

    E = distance(A,B)
    
        A - (DxM) matrix
        B - (DxN) matrix
    
    Returns:
        E - (MxN) Euclidean distances between vectors in A and B
    
    
    Description :
        This fully vectorized (VERY FAST!) script computes the
        Euclidean distance between two vectors by:
    
                     ||A-B|| = sqrt ( ||A||^2 + ||B||^2 - 2*A.B )
    
    Example :
        A = rand(400,100); B = rand(400,200);
        d = distance(A,B);
    
    Credits for original Matlab script:
    % Author   : Roland Bunschoten
    %            University of Amsterdam
    %            Intelligent Autonomous Systems (IAS) group
    %            Kruislaan 403  1098 SJ Amsterdam
    %            tel.(+31)20-5257524
    %            bunschot@wins.uva.nl
    % Last Rev : Oct 29 16:35:48 MET DST 1999
    % Tested   : PC Matlab v5.2 and Solaris Matlab v5.3
    % Thanx    : Nikos Vlassis
    
    % Copyright notice: You are free to modify, extend and distribute
    %    this code granted that the author of the original code is
    %    mentioned as the original author of the code.
    """
    from numpy import multiply, transpose, sum, dot, abs, sqrt
    from numpy.matlib import repmat

    assert a.shape[0] == b.shape[0], 'A and B should be of same dimensionality'

    # aa=sum(a.*a,1); bb=sum(b.*b,1); ab=a'*b;
    #d = sqrt(abs(repmat(aa',[1 size(bb,2)]) + repmat(bb,[size(aa,2) 1]) - 2*ab));

    aa = sum(multiply(a,a), axis=0)      # aa = sum(a.*a,1)
    bb = sum(multiply(b,b), axis=0)      # bb = sum(b.*b,1)
    ab = dot(transpose(a),b)             # ab = a'*b   
    return sqrt(abs(transpose(repmat(aa,bb.shape[0],1)) + repmat(bb,aa.shape[0],1) - 2*ab))
*/

template<class ArrayT>
inline ArrayT*
MatNorm2DistInternal(ArrayT* aT, ArrayT* b)
{
#ifdef SSE_USED
    return MatNorm2DistSSE(aT, b);
#endif
    ILOG_VAR(Core.Matrix.MatNorm2Dist);
    if (MatNrCol(aT) != MatNrRow(b)) {
        ILOG_ERROR("MatNorm2DistInternal operands: dimensionality problem");
    }

    Timer timer;
    
    // aa = sum(multiply(a,a), axis=0)      # aa = sum(a.*a,1)
    ArrayT* tempAT = 0;
    Mul(tempAT, aT, aT);         // elementwise multiplication
    ArrayT* aaT = MatSumAxis1(tempAT);
    delete tempAT;

    ILOG_DEBUG("aaT: " << timer.SplitTimeStr());

    // bb = sum(multiply(b,b), axis=0)      # bb = sum(b.*b,1)
    ArrayT* tempB = 0;
    Mul(tempB, b, b);         // elementwise multiplication
    ArrayT* bb = MatSumAxis0(tempB);
    delete tempB;

    ILOG_DEBUG("bb: " << timer.SplitTimeStr());

    // ab = dot(transpose(a),b)             # ab = a'*b   
    //ArrayT* aT = MatTranspose(a);
    //ILOG_DEBUG("transpose: " << timer.SplitTimeStr());

    ILOG_DEBUG("aT " << MatNrRow(aT) << " " << MatNrCol(aT));
    ILOG_DEBUG("b  " << MatNrRow(b) << " " << MatNrCol(b));

    ArrayT* ab = MatMul(aT, b);

    ILOG_DEBUG("matmul: " << timer.SplitTimeStr());

    // return sqrt(abs(transpose(repmat(aa,bb.shape[0],1)) + repmat(bb,aa.shape[0],1) - 2*ab))
    // #d = sqrt(abs(repmat(aa',[1 size(bb,2)]) + repmat(bb,[size(aa,2) 1]) - 2*ab));

    //ILOG_DEBUG("aaT " << MatNrRow(aaT) << " " << MatNrCol(aaT));
    //ILOG_DEBUG("bb  " << MatNrRow(bb) << " " << MatNrCol(bb));

    #pragma omp parallel for 
    for (int i=0 ; i<MatNrRow(aaT) ; i++)
    {
        Real64 tmp = *MatE(aaT, i, 0);
        for (int j=0 ; j<MatNrCol(bb) ; j++)
        {
           *MatE(ab, i, j) = sqrt(fabs(*MatE(bb, 0, j) + tmp - 2.0 * *MatE(ab, i, j)));
        }
    }
    ILOG_DEBUG("sqrt(abs(aa + bb - 2ab)): " << timer.SplitTimeStr());
    delete bb;
    delete aaT;

    ILOG_DEBUG("delete: " << timer.SplitTimeStr());
//WriteRaw(a, "matrix_a.raw", &Util::Database::GetInstance(), 1);
//WriteRaw(b, "matrix_b.raw", &Util::Database::GetInstance(), 1);
//WriteRaw(repAAT, "matrix_c.raw", &Util::Database::GetInstance(), 1);
    //ILOG_INFO("cpu-matnorm2dist-total: " << timer.SplitTime());
    ILOG_DEBUG(timer.SplitTime() << " (cpu-matnorm2dist-total)");
    return ab;
}


template<class ArrayT>
inline ArrayT*
MatNorm2Dist(ArrayT* a, ArrayT* b)
{
/*
    Computes Euclidean distance matrix

    E = MatNorm2Dist(A,B)
    
        A - (DxM) matrix
        B - (DxN) matrix
    
    Returns:
        E - (MxN) Euclidean distances between vectors in A and B
*/
    ILOG_VAR(Core.Matrix.MatNorm2Dist);
    if (MatNrRow(a) != MatNrRow(b)) {
        ILOG_ERROR("MatNorm2Dist operands: A and B should be of same dimensionality");
    }
    Mat* aT = MatTranspose(a);
    ArrayT* result = MatNorm2DistInternal(aT, b);
    delete aT;
    return result;
}


template<class ArrayT>
inline ArrayT*
MatNorm2DistTransposed(ArrayT* aT, ArrayT* bT)
{
/*
    Computes Euclidean distance matrix

    E = MatNorm2DistTransposed(A,B)
    
        A - (MxD) matrix
        B - (NxD) matrix
    
    Returns:
        E - (MxN) Euclidean distances between vectors in A and B
*/
    ILOG_VAR(Core.Matrix.MatNorm2DistTransposed);
    if (MatNrCol(aT) != MatNrCol(bT)) {
        ILOG_ERROR("MatNorm2DistTransposed operands: A and B should be of same dimensionality");
    }
#ifdef CUDA
    if(Link::Cuda::CudaUsed())
    {
        return Link::Cuda::MatNorm2DistTransposed(aT, bT);
    }
#endif
    Mat* b = MatTranspose(bT);
    Mat* result = MatNorm2DistInternal(aT, b);
    delete b;
    return result;
}

} // namespace Matrix
} // namespace Core
} // namespace Impala

#endif

---