MVAWrapperBase.h

Go to the documentation of this file.

// \version 
//
// \brief Helper functions for MVAReader and MVAWriter wrappers
//
// \author robert.sulej@cern.ch
//
#ifndef ANAB_MVAWRAPPERBASE_H
#define ANAB_MVAWRAPPERBASE_H

#include "canvas/Persistency/Common/Ptr.h"

#include "lardataobj/AnalysisBase/MVAOutput.h"

#include <typeinfo>
#include <functional>
#include <string>
#include <unordered_map>
#include <cmath>

namespace anab {

class FVectorWrapperBase {
public:

protected:

    std::string getProductName(std::type_info const & ti) const;
    size_t getProductHash(std::type_info const & ti) const { return ti.hash_code(); }

};

class MVAWrapperBase {
public:

protected:

    // all mva outputs in the feature vecor sum up to p=1

    template <class T, size_t N>
    std::array<float, N> pAccumulate(
        std::vector< art::Ptr<T> > const & items,
        std::vector< FeatureVector<N> > const & outs) const;

    template <class T, size_t N>
    std::array<float, N> pAccumulate(
        std::vector< art::Ptr<T> > const & items, std::vector<float> const & weights,
        std::vector< FeatureVector<N> > const & outs) const;

    template <class T, size_t N>
    std::array<float, N> pAccumulate(
        std::vector< art::Ptr<T> > const & items, std::function<float (T const &)> fweight,
        std::vector< FeatureVector<N> > const & outs) const;

    template <class T, size_t N>
    std::array<float, N> pAccumulate(
        std::vector< art::Ptr<T> > const & items, std::function<float (art::Ptr<T> const &)> fweight,
        std::vector< FeatureVector<N> > const & outs) const;

    // outputs in the feature vecor sum up to p=1 in groups:
    // - members of a group are tagged in the mask with the same non-negative number
    // - entries with negative tag in the mask are ignored

    template <class T, size_t N>
    std::array<float, N> pAccumulate(
        std::vector< art::Ptr<T> > const & items,
        std::vector< FeatureVector<N> > const & outs,
        std::array<char, N> const & mask) const;
};

} // namespace anab

//----------------------------------------------------------------------------
// MVAReader functions.
//
template <class T, size_t N>
std::array<float, N> anab::MVAWrapperBase::pAccumulate(
    std::vector< art::Ptr<T> > const & items,
    std::vector< anab::FeatureVector<N> > const & outs) const
{
    std::array<double, N> acc;
    acc.fill(0);

        float pmin = 1.0e-6, pmax = 1.0 - pmin;
        float log_pmin = std::log(pmin), log_pmax = std::log(pmax);

        for (auto const & ptr : items)
        {
                auto const & vout = outs[ptr.key()];
                for (size_t i = 0; i < vout.size(); ++i)
                {
                    float v;
                        if (vout[i] < pmin) v = log_pmin;
                        else if (vout[i] > pmax) v = log_pmax;
                        else v = std::log(vout[i]);

                        acc[i] += v;
                }
        }

        if (!items.empty())
        {
                double totp = 0.0;
                for (size_t i = 0; i < N; ++i)
                {
                        acc[i] = exp(acc[i] / items.size());
                        totp += acc[i];
                }
                for (size_t i = 0; i < N; ++i)
                {
                        acc[i] /= totp;
                }
        }
        else std::fill(acc.begin(), acc.end(), 1.0 / N);


    std::array<float, N> result;
    for (size_t i = 0; i < N; ++i) result[i] = acc[i];
    return result;
}
//----------------------------------------------------------------------------

template <class T, size_t N>
std::array<float, N> anab::MVAWrapperBase::pAccumulate(
    std::vector< art::Ptr<T> > const & items, std::vector<float> const & weights,
    std::vector< anab::FeatureVector<N> > const & outs) const
{
    std::array<double, N> acc;
    acc.fill(0);

        float pmin = 1.0e-6, pmax = 1.0 - pmin;
        float log_pmin = std::log(pmin), log_pmax = std::log(pmax);
        double totw = 0.0;

        for (size_t k = 0; k < items.size(); ++k)
        {
            auto const & ptr = items[k];
                float w = weights[k];

                if (w == 0) continue;

                auto const & vout = outs[ptr.key()];
                for (size_t i = 0; i < vout.size(); ++i)
                {
                    float v;
                        if (vout[i] < pmin) v = log_pmin;
                        else if (vout[i] > pmax) v = log_pmax;
                        else v = std::log(vout[i]);

                        acc[i] += w * v;
                }
                totw += w;
        }

        if (!items.empty())
        {
                double totp = 0.0;
                for (size_t i = 0; i < N; ++i)
                {
                        acc[i] = exp(acc[i] / totw);
                        totp += acc[i];
                }
                for (size_t i = 0; i < N; ++i)
                {
                        acc[i] /= totp;
                }
        }
        else std::fill(acc.begin(), acc.end(), 1.0 / N);


    std::array<float, N> result;
    for (size_t i = 0; i < N; ++i) result[i] = acc[i];
    return result;
}
//----------------------------------------------------------------------------

template <class T, size_t N>
std::array<float, N> anab::MVAWrapperBase::pAccumulate(
    std::vector< art::Ptr<T> > const & items, std::function<float (T const &)> fweight,
    std::vector< anab::FeatureVector<N> > const & outs) const
{
    std::array<double, N> acc;
    acc.fill(0);

        float pmin = 1.0e-6, pmax = 1.0 - pmin;
        float log_pmin = std::log(pmin), log_pmax = std::log(pmax);
        double totw = 0.0;

        for (size_t k = 0; k < items.size(); ++k)
        {
            auto const & ptr = items[k];
                float w = fweight(*ptr);

                if (w == 0) continue;

                auto const & vout = outs[ptr.key()];
                for (size_t i = 0; i < vout.size(); ++i)
                {
                    float v;
                        if (vout[i] < pmin) v = log_pmin;
                        else if (vout[i] > pmax) v = log_pmax;
                        else v = std::log(vout[i]);

                        acc[i] += w * v;
                }
                totw += w;
        }

        if (!items.empty())
        {
                double totp = 0.0;
                for (size_t i = 0; i < N; ++i)
                {
                        acc[i] = exp(acc[i] / totw);
                        totp += acc[i];
                }
                for (size_t i = 0; i < N; ++i)
                {
                        acc[i] /= totp;
                }
        }
        else std::fill(acc.begin(), acc.end(), 1.0 / N);


    std::array<float, N> result;
    for (size_t i = 0; i < N; ++i) result[i] = acc[i];
    return result;
}
//----------------------------------------------------------------------------

template <class T, size_t N>
std::array<float, N> anab::MVAWrapperBase::pAccumulate(
    std::vector< art::Ptr<T> > const & items, std::function<float (art::Ptr<T> const &)> fweight,
    std::vector< anab::FeatureVector<N> > const & outs) const
{
    std::array<double, N> acc;
    acc.fill(0);

        float pmin = 1.0e-6, pmax = 1.0 - pmin;
        float log_pmin = std::log(pmin), log_pmax = std::log(pmax);
        double totw = 0.0;

        for (size_t k = 0; k < items.size(); ++k)
        {
            auto const & ptr = items[k];
                float w = fweight(ptr);

                if (w == 0) continue;

                auto const & vout = outs[ptr.key()];
                for (size_t i = 0; i < vout.size(); ++i)
                {
                    float v;
                        if (vout[i] < pmin) v = log_pmin;
                        else if (vout[i] > pmax) v = log_pmax;
                        else v = std::log(vout[i]);

                        acc[i] += w * v;
                }
                totw += w;
        }

        if (!items.empty())
        {
                double totp = 0.0;
                for (size_t i = 0; i < N; ++i)
                {
                        acc[i] = exp(acc[i] / totw);
                        totp += acc[i];
                }
                for (size_t i = 0; i < N; ++i)
                {
                        acc[i] /= totp;
                }
        }
        else std::fill(acc.begin(), acc.end(), 1.0 / N);


    std::array<float, N> result;
    for (size_t i = 0; i < N; ++i) result[i] = acc[i];
    return result;
}
//----------------------------------------------------------------------------

//----------------------------------------------------------------------------
// functions with the mask tagging groups og labels
//----------------------------------------------------------------------------

template <class T, size_t N>
std::array<float, N> anab::MVAWrapperBase::pAccumulate(
    std::vector< art::Ptr<T> > const & items,
    std::vector< anab::FeatureVector<N> > const & outs,
    std::array<char, N> const & mask) const
{
    size_t n_groups = 0;
    std::unordered_map<char, size_t> label2group;
    std::vector<size_t> nb_entries;
    std::array<int, N> groupidx;
    for (size_t i = 0; i < N; ++i)
    {
        int idx = -1;
        if (mask[i] >= 0)
        {
            auto search = label2group.find(mask[i]);
            if (search == label2group.end())
            {
                idx = n_groups;
                label2group[mask[i]] = idx;
                nb_entries.push_back(0);
                ++n_groups;
            }
            else
            {
                idx = search->second;
                nb_entries[idx]++;
            }
        }
        groupidx[i] = idx;
    }

    std::array<double, N> acc;
    acc.fill(0);

        float pmin = 1.0e-6, pmax = 1.0 - pmin;
        float log_pmin = std::log(pmin), log_pmax = std::log(pmax);

        for (auto const & ptr : items)
        {
                auto const & vout = outs[ptr.key()];
                for (size_t i = 0; i < vout.size(); ++i)
                {
                    if (groupidx[i] < 0) continue;

                    float v;
                        if (vout[i] < pmin) v = log_pmin;
                        else if (vout[i] > pmax) v = log_pmax;
                        else v = std::log(vout[i]);

                        acc[i] += v;
                }
        }

        if (!items.empty())
        {
                std::vector<double> totp(n_groups, 0.0);
                for (size_t i = 0; i < N; ++i)
                {
                    if (groupidx[i] >= 0)
            {
                        acc[i] = exp(acc[i] / items.size());
                        totp[groupidx[i]] += acc[i];
                }
                }
                for (size_t i = 0; i < N; ++i)
                {
                        if (groupidx[i] >= 0) { acc[i] /= totp[groupidx[i]]; }
                }
        }
        else
        {
            for (size_t i = 0; i < N; ++i)
            {
                 if (groupidx[i] >= 0) { acc[i] = 1 / nb_entries[groupidx[i]]; }
            }
    }

    std::array<float, N> result;
    for (size_t i = 0; i < N; ++i) result[i] = acc[i];
    return result;
}
//----------------------------------------------------------------------------

#endif //ANAB_MVAWRAPPERBASE