JaakkolaHeuristic.h

Go to the documentation of this file.

//===========================================================================
/*!
 * 
 *
 * \brief       Jaakkola's heuristic and related quantities for Gaussian kernel selection
 * 
 * 
 *
 * \author      T. Glasmachers, O. Krause, C. Igel
 * \date        2010
 *
 *
 * \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/>.
 *
 */
//===========================================================================


#ifndef SHARK_ALGORITHMS_JAAKKOLAHEURISTIC_H
#define SHARK_ALGORITHMS_JAAKKOLAHEURISTIC_H


#include <shark/Data/Dataset.h>
#include <shark/Core/Traits/ProxyReferenceTraits.h>

#include <boost/range/adaptor/filtered.hpp>
#include <algorithm>

namespace shark{


/// \brief Jaakkola's heuristic and related quantities for Gaussian kernel selection
///
/// \par
/// Jaakkola's heuristic method for setting the width parameter of the
/// Gaussian radial basis function kernel is to pick a quantile (usually
/// the median) of the distribution of Euclidean distances between points
/// having different labels. The present implementation computes the kernel
/// width \f$ \sigma \f$ and the bandwidth
///    \f[ \gamma = \frac{1}{2 \sigma^2} \f]
/// based on the median or on any other quantile of the empirical
/// distribution.
///
/// By default, only the distance to the closest point with different
/// label is considered. This behavior can be turned off by an option
/// of the constructor. This is faster andin accordance with the
/// original paper.
class JaakkolaHeuristic
{
public:
    /// Constructor
    /// \param dataset           vector-valued input data
    /// \param nearestFalseNeighbor  if true, only the nearest neighboring point with different label is considered (default true)
    template<class InputType>
    JaakkolaHeuristic(LabeledData<InputType,unsigned int> const& dataset, bool nearestFalseNeighbor = true)
    {
        typedef typename LabeledData<InputType,unsigned int>::const_element_range Elements;
        typedef typename ConstProxyReference<InputType const>::type Element;
        Elements elements = dataset.elements();
        if(!nearestFalseNeighbor) {
            for(typename Elements::iterator it = elements.begin(); it != elements.end(); ++it){
                Element x = it->input;
                typename Elements::iterator itIn = it;
                itIn++;
                for (; itIn != elements.end(); itIn++) {
                    if (itIn->label == it->label) continue;
                    Element y = itIn->input;
                    double dist = distanceSqr(x,y);
                    m_stat.push_back(dist);
                }
            }

        } else {
            std::size_t classes = numberOfClasses(dataset);
            std::size_t dim = inputDimension(dataset);
            m_stat.resize(dataset.numberOfElements());
            std::fill(m_stat.begin(),m_stat.end(), std::numeric_limits<double>::max());
            std::size_t blockStart = 0;
            for(std::size_t c = 0; c != classes; ++c){
                
                typename Elements::iterator leftIt = elements.begin();
                typename Elements::iterator end = elements.end();
                while(leftIt != end){
                    //todo: use a filter on the iterator
                    //create the next batch containing only elements of class c as left argument to distanceSqr
                    typename Batch<InputType>::type leftBatch(512, dim);
                    std::size_t leftElements = 0;
                    while(leftElements < 512 && leftIt != end){
                        if(leftIt->label == c){
                            row(leftBatch,leftElements) = leftIt->input;
                            ++leftElements;
                        }
                        ++leftIt;
                    }
                    //now go through all elements and again create batches, this time of all elements which are not of class c
                    typename Elements::iterator rightIt = elements.begin();
                    while(rightIt != end){
                        typename Batch<InputType>::type rightBatch(512, dim);
                        std::size_t rightElements = 0;
                        while(rightElements < 512 && rightIt != end){
                            if(rightIt->label != c){
                                row(rightBatch,rightElements) = rightIt->input;
                                ++rightElements;
                            }
                            ++rightIt;
                        }

                        //now compute distances and update shortest distance
                        RealMatrix distances = distanceSqr(leftBatch,rightBatch);
                        for(std::size_t i = 0; i != leftElements;++i){
                            m_stat[blockStart+i]=std::min(min(subrange(row(distances,i),0,rightElements)),m_stat[blockStart+i]);
                        }
                    }
                    blockStart+= leftElements;
                }
            }
            
            //~ for(typename Elements::iterator it = elements.begin(); it != elements.end(); ++it){
                //~ double minDistSqr = std::numeric_limits<double>::max();//0;
                //~ Element x = it->input;
                //~ for (typename Elements::iterator itIn = elements.begin(); itIn != elements.end(); itIn++) {
                    //~ if (itIn->label == it->label) continue;
                    //~ Element y = itIn->input;
                    //~ double dist = distanceSqr(x,y);
                    //~ //if( (minDistSqr == 0) || (dist < minDistSqr))  minDistSqr = dist;
                    //~ if(dist < minDistSqr)  minDistSqr = dist;
                //~ }
                //~ m_stat.push_back(minDistSqr);
            //~ }
            
        }
        std::sort(m_stat.begin(), m_stat.end());
    }
        
    /// Compute the given quantile (usually median)
    /// of the empirical distribution of Euclidean distances
    /// of data pairs with different labels.
    double sigma(double quantile = 0.5)
    {
        std::size_t ic = m_stat.size();
        SHARK_ASSERT(ic > 0);

        std::sort(m_stat.begin(), m_stat.end());

        if (quantile < 0.0)
        {
            // TODO: find minimum
            return std::sqrt(m_stat[0]);
        }
        if (quantile >= 1.0)
        {
            // TODO: find maximum
            return std::sqrt(m_stat[ic-1]);
        }
        else
        {
            // TODO: partial sort!
            double t = quantile * (ic - 1);
            std::size_t i = (std::size_t)floor(t);
            double rest = t - i;
            return ((1.0 - rest) * std::sqrt(m_stat[i]) + rest * std::sqrt(m_stat[i+1]));
        }
    }

    /// Compute the given quantile (usually the median)
    /// of the empirical distribution of Euclidean distances
    /// of data pairs with different labels converted into
    /// a value usable as the gamma parameter of the GaussianRbfKernel.
    double gamma(double quantile = 0.5)
    {
        double s = sigma(quantile);
        return 0.5 / (s * s);
    }


private:
    /// all pairwise distances
    std::vector<double> m_stat;
};

}
#endif