LArSupportVectorMachine.h

Go to the documentation of this file.

#ifndef LAR_SUPPORT_VECTOR_MACHINE_H
#define LAR_SUPPORT_VECTOR_MACHINE_H 1

#include "larpandoracontent/LArHelpers/LArMvaHelper.h"

#include "larpandoracontent/LArObjects/LArMvaInterface.h"

#include "Pandora/StatusCodes.h"

#include <functional>
#include <map>
#include <vector>

//------------------------------------------------------------------------------------------------------------------------------------------

namespace lar_content
{

class SupportVectorMachine : public MvaInterface
{
public:
    typedef std::function<double(const LArMvaHelper::MvaFeatureVector &, const LArMvaHelper::MvaFeatureVector &, const double)> KernelFunction;

    enum KernelType
    {
        USER_DEFINED = 0,
        LINEAR       = 1,
        QUADRATIC    = 2,
        CUBIC        = 3,
        GAUSSIAN_RBF = 4
    };

    SupportVectorMachine();

    pandora::StatusCode Initialize(const std::string &parameterLocation, const std::string &svmName);

    bool Classify(const LArMvaHelper::MvaFeatureVector &features) const;

    double CalculateClassificationScore(const LArMvaHelper::MvaFeatureVector &features) const;

    double CalculateProbability(const LArMvaHelper::MvaFeatureVector &features) const;

    bool IsInitialized() const;

    unsigned int GetNFeatures() const;

    void SetKernelFunction(KernelFunction kernelFunction);

private:
    class SupportVectorInfo
    {
    public:
        SupportVectorInfo(const double yAlpha, LArMvaHelper::MvaFeatureVector supportVector);

        double                          m_yAlpha;        
        LArMvaHelper::MvaFeatureVector  m_supportVector; 
    };

    class FeatureInfo
    {
    public:
        FeatureInfo(const double muValue, const double sigmaValue);

        FeatureInfo();

        double StandardizeParameter(const double parameter) const;

        double    m_muValue;       
        double    m_sigmaValue;    
    };

    typedef std::vector<SupportVectorInfo> SVInfoList;
    typedef std::vector<FeatureInfo>       FeatureInfoVector;

    typedef std::map<KernelType, KernelFunction> KernelMap;

    bool              m_isInitialized;       

    bool              m_enableProbability;   
    double            m_probAParameter;      
    double            m_probBParameter;      

    bool              m_standardizeFeatures; 
    unsigned int      m_nFeatures;           
    double            m_bias;                
    double            m_scaleFactor;         

    SVInfoList        m_svInfoList;          
    FeatureInfoVector m_featureInfoList;     

    KernelType        m_kernelType;          
    KernelFunction    m_kernelFunction;      
    KernelMap         m_kernelMap;           

    void ReadXmlFile(const std::string &svmFileName, const std::string &svmName);

    pandora::StatusCode ReadComponent(pandora::TiXmlElement *pCurrentXmlElement);

    pandora::StatusCode ReadMachine(const pandora::TiXmlHandle &currentHandle);

    pandora::StatusCode ReadFeatures(const pandora::TiXmlHandle &currentHandle);

    pandora::StatusCode ReadSupportVector(const pandora::TiXmlHandle &currentHandle);

    double CalculateClassificationScoreImpl(const LArMvaHelper::MvaFeatureVector &features) const;

    static double QuadraticKernel(const LArMvaHelper::MvaFeatureVector &supportVector, const LArMvaHelper::MvaFeatureVector &features, const double scaleFactor = 1.);

    static double CubicKernel(const LArMvaHelper::MvaFeatureVector &supportVector, const LArMvaHelper::MvaFeatureVector &features, const double scaleFactor = 1.);

    static double LinearKernel(const LArMvaHelper::MvaFeatureVector &supportVector, const LArMvaHelper::MvaFeatureVector &features, const double scaleFactor = 1.);

    static double GaussianRbfKernel(const LArMvaHelper::MvaFeatureVector &supportVector, const LArMvaHelper::MvaFeatureVector &features, const double scaleFactor = 1.);
};

//------------------------------------------------------------------------------------------------------------------------------------------

inline bool SupportVectorMachine::Classify(const LArMvaHelper::MvaFeatureVector &features) const
{
    return (this->CalculateClassificationScoreImpl(features) > 0.);
}

//------------------------------------------------------------------------------------------------------------------------------------------

inline double SupportVectorMachine::CalculateClassificationScore(const LArMvaHelper::MvaFeatureVector &features) const
{
    return this->CalculateClassificationScoreImpl(features);
}

//------------------------------------------------------------------------------------------------------------------------------------------

inline double SupportVectorMachine::CalculateProbability(const LArMvaHelper::MvaFeatureVector &features) const
{
    if (!m_enableProbability)
    {
        std::cout << "LArSupportVectorMachine: cannot calculate probabilities for this SVM" << std::endl;
        throw pandora::STATUS_CODE_NOT_INITIALIZED;
    }

    // Use the logistic function to map the linearly-transformed score on the interval (-inf,inf) to a probability on [0,1] - the two free
    // parameters in the linear transformation are trained such that the logistic map produces an accurate probability
    const double scaledScore = m_probAParameter * this->CalculateClassificationScoreImpl(features) + m_probBParameter;

    return 1. / (1. + std::exp(scaledScore));
}

//------------------------------------------------------------------------------------------------------------------------------------------

inline bool SupportVectorMachine::IsInitialized() const
{
    return m_isInitialized;
}

//------------------------------------------------------------------------------------------------------------------------------------------

inline unsigned int SupportVectorMachine::GetNFeatures() const
{
    return m_nFeatures;
}

//------------------------------------------------------------------------------------------------------------------------------------------

inline void SupportVectorMachine::SetKernelFunction(KernelFunction kernelFunction)
{
    m_kernelFunction = std::move(kernelFunction);
}

//------------------------------------------------------------------------------------------------------------------------------------------

inline double SupportVectorMachine::LinearKernel(const LArMvaHelper::MvaFeatureVector &supportVector, const LArMvaHelper::MvaFeatureVector &features, const double scaleFactor)
{
    const double denominator(scaleFactor * scaleFactor);
    if (denominator < std::numeric_limits<double>::epsilon())
        throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);

    double total(0.);
    for (unsigned int i = 0; i < features.size(); ++i)
        total += supportVector.at(i).Get() * features.at(i).Get();

    return total / denominator;
}

//------------------------------------------------------------------------------------------------------------------------------------------

inline double SupportVectorMachine::QuadraticKernel(const LArMvaHelper::MvaFeatureVector &supportVector, const LArMvaHelper::MvaFeatureVector &features, const double scaleFactor)
{
    const double denominator(scaleFactor * scaleFactor);
    if (denominator < std::numeric_limits<double>::epsilon())
        throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);

    double total(0.);
    for (unsigned int i = 0; i < features.size(); ++i)
        total += supportVector.at(i).Get() * features.at(i).Get();

    total = total / denominator + 1.;
    return total * total;
}

//------------------------------------------------------------------------------------------------------------------------------------------

inline double SupportVectorMachine::CubicKernel(const LArMvaHelper::MvaFeatureVector &supportVector, const LArMvaHelper::MvaFeatureVector &features, const double scaleFactor)
{
    const double denominator(scaleFactor * scaleFactor);
    if (denominator < std::numeric_limits<double>::epsilon())
        throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);

    double total(0.);
    for (unsigned int i = 0; i < features.size(); ++i)
        total += supportVector.at(i).Get() * features.at(i).Get();

    total = total / denominator + 1.;
    return total * total * total;
}

//------------------------------------------------------------------------------------------------------------------------------------------

inline double SupportVectorMachine::GaussianRbfKernel(const LArMvaHelper::MvaFeatureVector &supportVector, const LArMvaHelper::MvaFeatureVector &features, const double scaleFactor)
{
    double total(0.);
    for (unsigned int i = 0; i < features.size(); ++i)
        total += (supportVector.at(i).Get() - features.at(i).Get()) * (supportVector.at(i).Get() - features.at(i).Get());

    return std::exp(-scaleFactor * total);
}

//------------------------------------------------------------------------------------------------------------------------------------------

inline SupportVectorMachine::SupportVectorInfo::SupportVectorInfo(const double yAlpha, LArMvaHelper::MvaFeatureVector supportVector) :
    m_yAlpha(yAlpha),
    m_supportVector(std::move(supportVector))
{
}

//------------------------------------------------------------------------------------------------------------------------------------------

inline SupportVectorMachine::FeatureInfo::FeatureInfo(const double muValue, const double sigmaValue) :
    m_muValue(muValue),
    m_sigmaValue(sigmaValue)
{
}

//------------------------------------------------------------------------------------------------------------------------------------------

inline SupportVectorMachine::FeatureInfo::FeatureInfo() :
    m_muValue(0.),
    m_sigmaValue(0.)
{
}

//------------------------------------------------------------------------------------------------------------------------------------------

inline double SupportVectorMachine::FeatureInfo::StandardizeParameter(const double parameter) const
{
    if (m_sigmaValue < std::numeric_limits<double>::epsilon())
        throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);

    return (parameter - m_muValue) / m_sigmaValue;
}

} // namespace lar_content

#endif // #ifndef LAR_SUPPORT_VECTOR_MACHINE_H