helper_funcs.hpp

Go to the documentation of this file.

#ifndef RBC_HELPERFUNCS_HPP
#define RBC_HELPERFUNCS_HPP

#include <cassert>
#include <algorithm>
#include <functional>
#include <RBC/data_types.hpp>

#if defined(__APPLE__) || defined(__MACOSX)
#include <OpenCL/cl.hpp>
#else
#include <CL/cl.hpp>
#endif


namespace RBC
{

    bool setProfilingFlag (int argc, char **argv);


    template <typename T>
    uint64_t nextPow2 (T num)
    {
        assert (num >= 0);

        uint64_t pow;
        for (pow = 1; pow < (uint64_t) num; pow <<= 1) ;

        return pow;
    }


    template <typename T>
    void printBuffer (const char *title, T *ptr, uint32_t width, uint32_t height)
    {
        std::cout << title << std::endl;

        for (int row = 0; row < height; ++row)
        {
            for (int col = 0; col < width; ++col)
            {
                std::cout << std::setw (3 * sizeof (T)) << +ptr[row * width + col] << " ";
            }
            std::cout << std::endl;
        }

        std::cout << std::endl;
    }


    template <typename T>
    void printBufferF (const char *title, T *ptr, uint32_t width, uint32_t height, uint32_t prec)
    {
        std::ios::fmtflags f (std::cout.flags ());
        std::cout << title << std::endl;
        std::cout << std::fixed << std::setprecision (prec);

        for (int row = 0; row < height; ++row)
        {
            for (int col = 0; col < width; ++col)
            {
                std::cout << std::setw (5 + prec) << ptr[row * width + col] << " ";
            }
            std::cout << std::endl;
        }

        std::cout << std::endl;
        std::cout.flags (f);
    }


    template <typename T>
    void cpuReduce (T *in, T *out, uint32_t cols, uint32_t rows, std::function<bool (T, T)> func)
    {
        for (uint r = 0; r < rows; ++r)
        {
            T rec = in[r * cols];
            for (uint c = 1; c < cols; ++c)
            {
                T tmp = in[r * cols + c];
                if (func (tmp, rec)) rec = tmp;
            }
            out[r] = rec;
        }
    }


    template <typename T>
    void cpuInScan (T *in, T *out, uint32_t width, uint32_t height)
    {
        // Initialize the first element of each row
        for (uint32_t row = 0; row < height; ++row)
            out[row * width] = in[row * width];
        // Perform the scan
        for (uint32_t row = 0; row < height; ++row)
            for (uint32_t col = 1; col < width; ++col)
                out[row * width + col] = out[row * width + col - 1] + in[row * width + col];
    }


    template <typename T>
    void cpuExScan (T *in, T *out, uint32_t width, uint32_t height)
    {
        // Initialize the first element of each row
        for (uint32_t row = 0; row < height; ++row)
            out[row * width] = 0;
        // Perform the scan
        for (uint32_t row = 0; row < height; ++row)
            for (uint32_t col = 1; col < width; ++col)
                out[row * width + col] = out[row * width + col - 1] + in[row * width + col - 1];
    }


    template <typename T>
    void cpuRBCComputeDists (T *X, T *R, T *D, uint32_t nx, uint32_t nr, uint32_t d)
    {
        for (uint x = 0; x < nx; ++x)
        {
            for (uint r = 0; r < nr; ++r)
            {
                T dist = 0.f;
                for (uint j = 0; j < d; ++j)
                    dist += std::pow (X[x * d + j] - R[r * d + j], 2);

                D[x * nr + r] = dist;
            }
        }
    }


    template <typename T>
    void cpuRBCComputeDists8 (T *X, T *R, T *D, uint32_t nx, uint32_t nr, uint32_t d, T a)
    {
        for (uint x = 0; x < nx; ++x)
        {
            for (uint r = 0; r < nr; ++r)
            {
                T dist = 0.f;
                for (uint j = 0; j < 4; ++j)
                    dist += 1.f / (1.f + a) * std::pow (X[x * d + j] - R[r * d + j], 2);
                for (uint j = 4; j < 8; ++j)
                    dist += a * std::pow (X[x * d + j] - R[r * d + j], 2);

                D[x * nr + r] = dist;
            }
        }
    }


    template <typename T>
    void cpuRBCMinDists (T *in, rbc_dist_id *out, uint32_t *N, uint32_t *Rnk, 
                         uint32_t cols, uint32_t rows, bool accCounters)
    {
        if (accCounters)
            for (uint c = 0; c < cols; ++c)
                N[c] = 0;

        rbc_dist_id rbcMin;

        for (uint r = 0; r < rows; ++r)
        {
            rbcMin.dist = in[r * cols];
            rbcMin.id = 0;

            for (uint c = 1; c < cols; ++c)
            {
                T tmp = in[r * cols + c];
                if (tmp < rbcMin.dist)
                {
                    rbcMin.dist = tmp;
                    rbcMin.id = c;
                }
            }

            out[r] = rbcMin;

            if (accCounters)
                Rnk[r] = N[rbcMin.id]++;
        }
    }


    template <typename T>
    void cpuRBCPermute (T *X, rbc_dist_id *ID, T *Xp, rbc_dist_id *IDp, uint32_t *O, uint32_t *Rnk, 
                        uint32_t nx, uint32_t nr, uint32_t d, bool permID)
    {
        for (uint32_t x = 0; x < nx; ++x)
        {
            uint id = ID[x].id;
            uint offset = O[id];
            uint rank = Rnk[x];

            for (uint32_t j = 0; j < d; ++j)
                Xp[(offset + rank) * d + j] = X[x * d + j];

            if (permID == true)
                IDp[offset + rank] = ID[x];
        }
    }


    template <typename T>
    T euclideanMetric (T *p1, T *p2, uint32_t d)
    {
        T dist = 0.f;

        for (int i = 0; i < d; ++i)
            dist += std::pow (p1[i] - p2[i], 2);

        return std::sqrt (dist);
    }
    

    template <typename T>
    T euclideanMetric8Squared (T *p1, T *p2, float a)
    {
        T dist = 0.f;

        for (int i = 0; i < 4; ++i)
            dist += 1.f / (1.f + a) * std::pow (p1[i] - p2[i], 2);
        for (int i = 4; i < 8; ++i)
            dist += a * std::pow (p1[i] - p2[i], 2);

        return dist;
    }


    template <typename T>
    void cpuRBCSearch (T *Qp, rbc_dist_id *RID, T *Xp, cl_uint *O, cl_uint *N, rbc_dist_id *NNID, T *NN, 
                       uint32_t nq, uint32_t nr, uint32_t nx, uint32_t d)
    {
        cl_uint max_n = std::accumulate (N, N + nr, 0, 
            [](cl_uint a, cl_uint b) -> cl_uint { return std::max (a, b); });

        cl_float *D = new cl_float[nq * max_n];
        
        // Compute distances Q - X[L]
        for (uint q = 0; q < nq; ++q)
        {
            cl_uint rID = RID[q].id;
            cl_uint o = O[rID];
            cl_uint n = N[rID];

            for (uint x = 0; x < max_n; ++x)
            {
                float dist = 0.f;

                if (x < n)
                    for (uint j = 0; j < d; ++j)
                        dist += std::pow (Qp[q * d + j] - Xp[(o + x) * d + j], 2);
                else
                    dist = std::numeric_limits<float>::infinity ();

                D[q * max_n + x] = dist;
            }
        }
        
        cpuRBCMinDists (D, NNID, nullptr, nullptr, max_n, nq, false);

        for (uint32_t q = 0; q < nq; ++q)
        {
            uint nnID = O[RID[q].id] + NNID[q].id;

            for (uint32_t j = 0; j < d; ++j)
                NN[q * d + j] = Xp[nnID * d + j];
        }

        delete[] D;
    }


    template <typename T>
    void cpuRBCSearch8 (T *Qp, rbc_dist_id *RID, T *Xp, cl_uint *O, cl_uint *N, rbc_dist_id *NNID, T *NN, 
                        uint32_t nq, uint32_t nr, uint32_t nx, T a)
    {
        const unsigned int d = 8;
        cl_uint max_n = std::accumulate (N, N + nr, 0, 
            [](cl_uint a, cl_uint b) -> cl_uint { return std::max (a, b); });

        cl_float *D = new cl_float[nq * max_n];
        
        // Compute distances Q - X[L]
        for (uint q = 0; q < nq; ++q)
        {
            cl_uint rID = RID[q].id;
            cl_uint o = O[rID];
            cl_uint n = N[rID];

            for (uint x = 0; x < max_n; ++x)
            {
                float dist = 0.f;

                if (x < n)
                {
                    for (uint j = 0; j < 4; ++j)
                        dist += 1.f / (1.f + a) * std::pow (Qp[q * d + j] - Xp[(o + x) * d + j], 2);
                    for (uint j = 4; j < d; ++j)
                        dist += a * std::pow (Qp[q * d + j] - Xp[(o + x) * d + j], 2);
                }
                else
                    dist = std::numeric_limits<float>::infinity ();

                D[q * max_n + x] = dist;
            }
        }
        
        cpuRBCMinDists (D, NNID, nullptr, nullptr, max_n, nq, false);

        for (uint32_t q = 0; q < nq; ++q)
        {
            uint nnID = O[RID[q].id] + NNID[q].id;

            for (uint32_t j = 0; j < d; ++j)
                NN[q * d + j] = Xp[nnID * d + j];
        }

        delete[] D;
    }


    template <typename T>
    void cpuNNSearch (T *Q, T *X, T *NN, uint32_t nq, uint32_t nx, uint32_t d)
    {
        // Compute distance array
        std::vector<T> D (nq * nx);

        for (uint32_t q = 0; q < nq; ++q)
            for (uint32_t x = 0; x < nx; ++x)
                D[q * nx + x] = euclideanMetric (&Q[q * d], &X[x * d], d);

        // Find NN IDs
        std::vector<uint32_t> minID (nq);

        for (uint32_t q = 0; q < nq; ++q)
        {
            T minDist = D[q * nx];
            uint32_t id = 0;

            for (uint32_t x = 1; x < nx; ++x)
            {
                T tmp = D[q * nx + x];
                if (tmp < minDist)
                {
                    minDist = tmp;
                    id = x;
                }
            }

            minID[q] = id;
        }

        // Collect NNs
        for (uint32_t q = 0; q < nq; ++q)
            for (uint32_t j = 0; j < d; ++j)
                NN[q * d + j] = X[minID[q] * d + j];     
    }


    template <typename T>
    T meanError (T *Q, T *NN, uint32_t n, uint32_t d)
    {
        std::vector<T> D (n);
        unsigned int i = -1;

        std::generate (D.begin (), D.end (), 
            [&]() { i++; return euclideanMetric (&Q[i * d], &NN[i * d], d); });

        T sumD = std::accumulate (D.begin (), D.end (), 0.f);

        return sumD / (T) n;
    }

}

#endif  // RBC_HELPERFUNCS_HPP