PtPartitioningIndex.cc

#include <cassert>
#include <algorithm>
#include <limits.h>

#include "PtPartitioningIndex.h"

namespace stab {

    unsigned PtPartitioningIndex::max_element( const unsigned * clu, unsigned size ) const
    {
    unsigned max = 0U;
    for ( unsigned i=0; i != size; i++ ) {
      const unsigned tmp = clu[i];
      if ( tmp > max && tmp < UINT_MAX ) {
        max = tmp; 
      }
    }

    return max;

    }

    double PtPartitioningIndex::index(const PseudoJet* inputParticles,
                                      const unsigned* clustering1,
                                      const unsigned* clustering2,
                                      const unsigned len) const
    {
        assert(inputParticles);
        assert(clustering1);
        assert(clustering2);
        assert(len);

        //unsigned n1 = *std::max_element(clustering1, clustering1+len)+1U;
        //unsigned n2 = *std::max_element(clustering2, clustering2+len)+1U;
        unsigned n1 = this->max_element(clustering1,len)+1U;
    unsigned n2 = this->max_element(clustering2,len)+1U;

        if (n1 > n2)
        {
            std::swap(clustering1, clustering2);
            std::swap(n1, n2);
        }

        const unsigned matlen = n1*n2;
        assert(matlen);

        if (storage_.size() < matlen)
            storage_.resize(matlen);
        double* data = &storage_[0];

        for (unsigned i=0; i<matlen; ++i)
            data[i] = 0.0;

        long double ptsum = 0.0L;
        long double ptsumsq = 0.0L;
        for (unsigned i=0; i<len; ++i)
        {
        if ( clustering1[i] == UINT_MAX || clustering2[i] == UINT_MAX ) continue;
            const double pt = inputParticles[i].pt();
            ptsum += pt;
            ptsumsq += pt*pt;
            data[clustering1[i]*n2 + clustering2[i]] += pt;
        }
        assert(ptsum > 0.0L);

        return processInterClusterMatrix(data, n1, n2, ptsum, ptsumsq);
    }

    double BestPtMatchIndex::processInterClusterMatrix(
        const double* mat, const unsigned nrows,
        const unsigned ncols, const double etSum,
        double /* ptSumsq */) const
    {
        const unsigned matlen = nrows*ncols;
        if (storage_.size() < matlen)
            storage_.resize(matlen);
        ElementSort* data = &storage_[0];

        const unsigned nflags = nrows + ncols;
        if (intstorage_.size() < nflags)
            intstorage_.resize(nflags);
        int *rowTaken = &intstorage_[0];
        int *colTaken = rowTaken + nrows;
        for (unsigned i=0; i<nflags; ++i)
            rowTaken[i] = 0;

        ElementSort* d = data;
        for (unsigned i=0; i<nrows; ++i)
            for (unsigned j=0; j<ncols; ++j, ++d)
            {
                d->ptsum = *mat++;
                d->row = i;
                d->col = j;
            }

        std::sort(data, data+matlen);
        long double usedsum = 0.0L;
        ElementSort* const dmax = data + matlen;
        unsigned k = 0;
        const unsigned kmax = nrows < ncols ? nrows : ncols;
        for (d = data; d!=dmax && k<kmax; ++d)
            if (!rowTaken[d->row] && !colTaken[d->col])
            {
                rowTaken[d->row] = 1;
                colTaken[d->col] = 1;
                usedsum += d->ptsum;
                ++k;
            }

        return static_cast<double>(usedsum)/etSum;
    }

    double TransverseRandIndex2::processInterClusterMatrix(
        const double* mat, const unsigned nrows,
        const unsigned ncols, const double etSum,
        const double etSumsq) const
    {
        const unsigned nsums = nrows + ncols;
        if (storage_.size() < nsums)
            storage_.resize(nsums);
        long double* rowSums = &storage_[0];
        long double* colSums = rowSums + nrows;
        for (unsigned i=0; i<ncols; ++i)
            colSums[i] = 0.0L;

        long double sumsq = 0.0L;
        for (unsigned i=0; i<nrows; ++i)
        {
            long double rsum = 0.0L;
            const double* m = mat + i*ncols;
            for (unsigned j=0; j<ncols; ++j)
            {
                rsum += m[j];
                colSums[j] += m[j];
                sumsq += m[j]*m[j];
            }
            rowSums[i] = rsum;
        }

        long double crsq = 0.0L;
        for (unsigned i=0; i<nsums; ++i)
            crsq += rowSums[i]*rowSums[i];

        const double norm = 0.5*(etSum*etSum - etSumsq);
        assert(norm > 0.0);

        return 1.0 - (0.5*crsq - sumsq)/norm;
    }
}