MissingFeaturesKernelExpansion.h

Go to the documentation of this file.

//===========================================================================
/*!
 * 
 *
 * \brief       A kernel expansion with support of missing features
 * 
 * 
 *
 * \author      B. Li
 * \date        2012
 *
 *
 * \par Copyright 1995-2017 Shark Development Team
 * 
 * <BR><HR>
 * This file is part of Shark.
 * <http://shark-ml.org/>
 * 
 * Shark is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published 
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * Shark is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with Shark.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
//===========================================================================
#include <shark/Data/Dataset.h>
#include <shark/Data/DataView.h>
#include <shark/Models/Kernels/EvalSkipMissingFeatures.h>
#include <shark/Models/Kernels/KernelExpansion.h>

namespace shark {

/// Kernel expansion with missing features support
template<class InputType>
class MissingFeaturesKernelExpansion : public KernelExpansion<InputType>
{
private:
    typedef KernelExpansion<InputType> Base;
public:
    typedef typename Base::KernelType KernelType;
    typedef typename Base::BatchInputType BatchInputType;
    typedef typename Base::BatchOutputType BatchOutputType;
    /// Constructors from the base class
    ///@{
    MissingFeaturesKernelExpansion(){}


    MissingFeaturesKernelExpansion(KernelType* kernel)
    : Base(kernel)
    {}

    MissingFeaturesKernelExpansion(KernelType* kernel, Data<InputType> const& basis, bool offset)
    : Base(kernel, basis, offset, 1u)
    {}
    ///@}

    /// \brief From INameable: return the class name.
    std::string name() const
    { return "MissingFeaturesKernelExpansion"; }

    boost::shared_ptr<State> createState()const{
        return boost::shared_ptr<State>(new EmptyState());
    }

    /// Override eval(...) in the base class
    virtual void eval(BatchInputType const& patterns, BatchOutputType& outputs)const{
        SHARK_ASSERT(Base::mep_kernel);
        SIZE_CHECK(Base::m_alpha.size1() > 0u);
        
        //Todo: i am too lazy to us iterated loops in this function.
        //so i am using a DataView to have O(1) random access lookup. but this is not needed!
        DataView<Data<InputType> const > indexedBasis(Base::m_basis);
        
        ensure_size(outputs,batchSize(patterns),Base::outputShape().numElements());
        if (Base::hasOffset())
                noalias(outputs) = repeat(Base::m_b,batchSize(patterns));
            else
                outputs.clear();
        
        for(std::size_t p = 0; p != batchSize(patterns); ++p){


            // Calculate scaling coefficient for the 'pattern'
            const double patternNorm = computeNorm(column(Base::m_alpha, 0), m_scalingCoefficients, row(patterns,p));
            const double patternSc = patternNorm / m_classifierNorm;

            // Do normal classification except that we use kernel which supports inputs with Missing features
            //TODO: evaluate k for all i and replace the += with a matrix-vector operation. 
            //better: do this for all p and i and go matrix-matrix-multiplication
            for (std::size_t i = 0; i != indexedBasis.size(); ++i){
                const double k = evalSkipMissingFeatures(
                    *Base::mep_kernel,
                    indexedBasis[i],
                    row(patterns,p)) / m_scalingCoefficients[i] / patternSc;
                noalias(row(outputs,p)) += k * row(Base::m_alpha, i);
                
            }
        }
    }
    void eval(BatchInputType const& patterns, BatchOutputType& outputs, State & state)const{
        eval(patterns, outputs);
    }

    /// Calculate norm of classifier, i.e., ||w||
    ///
    /// formula:
    /// \f$ \sum_{i,j=1}^{n}\alpha_i\frac{y_i}{s_i}K\left(x_i,x_j)\right)\frac{y_j}{s_j}\alpha_j \f$
    /// where \f$ s_i \f$ is scaling coefficient, and \f$ K \f$ is kernel function,
    /// \f$ K\left(x_i,x_j)\right) \f$ is taken only over features that are valid for both \f$ x_i \f$ and \f$ x_j \f$
    template<class InputTypeT>
    double computeNorm(
        const RealVector& alpha,
        const RealVector& scalingCoefficient,
        InputTypeT const& missingness
    ) const{
        SHARK_ASSERT(Base::mep_kernel);
        SIZE_CHECK(alpha.size() == scalingCoefficient.size());
        SIZE_CHECK(Base::m_basis.numberOfElements() == alpha.size());

        // Calculate ||w||^2
        double norm_sqr = 0.0;
        
        //Todo: i am too lazy to use iterated loops in this function.
        //so i am using a DataView to have O(1) random access lookup. but this is not needed!
        DataView<Data<InputType> const > indexedBasis(Base::m_basis);

        for (std::size_t i = 0; i < alpha.size(); ++i){
            for (std::size_t j = 0; j < alpha.size(); ++j){
                const double evalResult = evalSkipMissingFeatures(
                    *Base::mep_kernel,
                    indexedBasis[i],
                    indexedBasis[j],
                    missingness);
                // Note that in Shark solver, we do axis flip by substituting \alpha with y \times \alpha
                norm_sqr += evalResult * alpha(i) * alpha(j) / scalingCoefficient(i) / scalingCoefficient(j);
            }
        }

        // Return ||w||
        return std::sqrt(norm_sqr);
    }
    
    double computeNorm(
        const RealVector& alpha,
        const RealVector& scalingCoefficient
    ) const{
        SHARK_ASSERT(Base::mep_kernel);
        SIZE_CHECK(alpha.size() == scalingCoefficient.size());
        SIZE_CHECK(Base::m_basis.numberOfElements() == alpha.size());
        
        //Todo: i am too lazy to us iterated loops in this function.
        //so i am using a DataView to have O(1) random access lookup. but this is not needed!
        DataView<Data<InputType> const > indexedBasis(Base::m_basis);

        // Calculate ||w||^2
        double norm_sqr = 0.0;
        
        for (std::size_t i = 0; i < alpha.size(); ++i){
            for (std::size_t j = 0; j < alpha.size(); ++j){
                const double evalResult = evalSkipMissingFeatures(
                    *Base::mep_kernel,
                    indexedBasis[i],
                    indexedBasis[j]);
                // Note that in Shark solver, we do axis flip by substituting \alpha with y \times \alpha
                norm_sqr += evalResult * alpha(i) * alpha(j) / scalingCoefficient(i) / scalingCoefficient(j);
            }
        }

        // Return ||w||
        return std::sqrt(norm_sqr);
    }

    void setScalingCoefficients(const RealVector& scalingCoefficients)
    {
#if DEBUG
        for(double v: scalingCoefficients)
        {
            SHARK_ASSERT(v > 0.0);
        }
#endif
        m_scalingCoefficients = scalingCoefficients;
    }

    void setClassifierNorm(double classifierNorm)
    {
        SHARK_ASSERT(classifierNorm > 0.0);
        m_classifierNorm = classifierNorm;
    }

protected:
    /// The scaling coefficients
    RealVector m_scalingCoefficients;

    /// The norm of classifier(w)
    double m_classifierNorm;
};

} // namespace shark {