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 "Helpers/XmlHelper.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