MAP.cpp

Go to the documentation of this file.

/*  Copyright (C) 2008 National Institute For Space Research (INPE) - Brazil.

    This file is part of the TerraLib - a Framework for building GIS enabled applications.

    TerraLib 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.

    TerraLib 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 TerraLib. See COPYING. If not, write to
    TerraLib Team at <terralib-team@terralib.org>.
 */

/*!
  \file terralib/classification/MAP.cpp

  \brief MAP (Maximum a Posteriori) strategy for classification.
*/

// TerraLib
#include "MAP.h"

// STL
#include <limits>

namespace te
{
  namespace cl
  {
    MAP::Parameters::Parameters()
    {
      reset();
    }

    MAP::Parameters::~Parameters() = default;

    const MAP::Parameters& MAP::Parameters::operator=(
      const MAP::Parameters& rhs)
    {
      reset();

      m_prioriProbs = rhs.m_prioriProbs;
      m_prioriCalcSampleStep = rhs.m_prioriCalcSampleStep;

      return *this;
    }

    void MAP::Parameters::reset() throw(Exception)
    {
      m_prioriProbs.clear();
      m_prioriCalcSampleStep = 5;
    }

    te::common::AbstractParameters* MAP::Parameters::clone() const
    {
      return new MAP::Parameters(*this);
    }

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

    MAP::MAP()
    {
      reset();
    }

    MAP::~MAP() = default;

    bool MAP::initialize(const MAP::Parameters& params) throw(Exception)
    {
      reset();

      if( !params.m_prioriProbs.empty() )
      {
        double sum = 0.0;

        for( std::size_t prioriProbsIdx = 0 ; prioriProbsIdx < params.m_prioriProbs.size() ;
          ++prioriProbsIdx )
        {
          if(
              ( params.m_prioriProbs[ prioriProbsIdx ] < 0.0 )
              ||
              ( params.m_prioriProbs[ prioriProbsIdx ] > 1.0 )
            )
          {
            return false;
          }

          sum += params.m_prioriProbs[ prioriProbsIdx ];
        }

        if( sum != 1.0 )
        {
          return false;
        }
      }

      if( params.m_prioriCalcSampleStep < 1 )
      {
        return false;
      }

      m_parameters = params;

      m_isInitialized = true;

      return true;
    }

    void MAP::reset()
    {
      m_isInitialized = false;
      m_parameters.reset();
      m_classesMeans.clear();
      m_classesCovarianceMatrixes.clear();
      m_classesCovarianceInvMatrixes.clear();
      m_classLabels.clear();
      m_classesOptizedMAPDiscriminantTerm.clear();
    }

    bool MAP::train( const InputAdaptor< double >& samples,
      const std::vector<unsigned int>& attributesIndices,
      const std::vector<unsigned int>& sampleLabels,
      const bool enableProgressInterface) throw(Exception)
    {
      if( ! m_isInitialized ) return false;
      if( samples.getElementsCount() == 0 ) return false;
      if( sampleLabels.size() != samples.getElementsCount() ) return false;
      if( attributesIndices.empty() ) return false;

      const unsigned int attributesIndicesSize = static_cast<unsigned int>(attributesIndices.size());

      std::unique_ptr< te::common::TaskProgress > progressPtr;
      if( enableProgressInterface )
      {
        progressPtr.reset( new te::common::TaskProgress );
        progressPtr->setTotalSteps( 3 );
        progressPtr->setMessage( "Trainning" );
      }

      // grouping samples by class

      std::vector< std::vector< std::vector< double > > > samplesByClass;

      {
        m_classLabels.clear();

        const unsigned int samplesCount = samples.getElementsCount();
        std::map< unsigned int, unsigned int > labels2ClassIndexMap;
        std::map< unsigned int, unsigned int >::iterator labels2ClassIndexMapIt;
        unsigned int attributesIndicesIdx = 0;

        for( unsigned int sampleIdx = 0 ; sampleIdx < samplesCount ; ++sampleIdx )
        {
          const unsigned int& sampleLabel = sampleLabels[ sampleIdx ];

          labels2ClassIndexMapIt = labels2ClassIndexMap.find( sampleLabel );

          if( labels2ClassIndexMapIt  == labels2ClassIndexMap.end() )
          {
            labels2ClassIndexMap[sampleLabel] = static_cast<unsigned int>(m_classLabels.size());
            m_classLabels.push_back( sampleLabel );

            samplesByClass.push_back( std::vector< std::vector< double > >() );
            samplesByClass.back().push_back( std::vector< double >() );

            std::vector< double >& sample = samplesByClass.back().back();
            sample.resize( attributesIndicesSize, 0.0 );

            for( attributesIndicesIdx = 0; attributesIndicesIdx < attributesIndicesSize;
              ++attributesIndicesIdx )
            {
              samples.getFeature( sampleIdx, attributesIndicesIdx, sample[ attributesIndicesIdx ] );
            }
          }
          else
          {
            samplesByClass[ labels2ClassIndexMapIt->second ].push_back( std::vector< double >() );
            std::vector< double >& sample = samplesByClass[ labels2ClassIndexMapIt->second ].back();
            sample.resize( attributesIndicesSize, 0.0 );

            for( attributesIndicesIdx = 0; attributesIndicesIdx < attributesIndicesSize;
              ++attributesIndicesIdx )
            {
              samples.getFeature( sampleIdx, attributesIndicesIdx, sample[ attributesIndicesIdx ] );
            }
          }
        }
      }

      if( enableProgressInterface )
      {
        progressPtr->pulse();
        if( ! progressPtr->isActive() ) return false;
      }

      // Checking priori probabilities

      if( ! m_parameters.m_prioriProbs.empty() )
      {
        if( m_parameters.m_prioriProbs.size() != samplesByClass.size() )
        {
          return false;
        }
      }

      // classes means calcule

      {
        const std::size_t samplesByClassSize = samplesByClass.size();

        m_classesMeans.resize( samplesByClassSize );
        std::vector< double > dymmyMeansVec( attributesIndicesSize, 0.0 );
        std::fill( m_classesMeans.begin(), m_classesMeans.end(), dymmyMeansVec );

        unsigned int attributeIdx = 0;
        double mean = 0;
        std::size_t classSamplesIdx = 0;
        std::size_t classSamplesSize = 0;

        for( std::size_t samplesByClassIdx = 0 ; samplesByClassIdx < samplesByClassSize ;
          ++samplesByClassIdx )
        {
          const std::vector< std::vector< double > >& classSamples = samplesByClass[
            samplesByClassIdx ];
          classSamplesSize = classSamples.size();

          for( attributeIdx = 0; attributeIdx < attributesIndicesSize; ++attributeIdx )
          {
            mean = 0;

            for( classSamplesIdx = 0 ; classSamplesIdx < classSamplesSize ; ++classSamplesIdx )
            {
              mean += classSamples[ classSamplesIdx ][ attributeIdx ];
            }

            mean /= (double)classSamplesSize;

            m_classesMeans[ samplesByClassIdx ][ attributeIdx ] = mean;
          }
        }
      }

      if( enableProgressInterface )
      {
        progressPtr->pulse();
        if( ! progressPtr->isActive() ) return false;
      }

      // covariance matrixes calcule

      {
        const std::size_t samplesByClassSize = samplesByClass.size();
        const boost::numeric::ublas::matrix< double > dummyCovarianceMatrix(
          attributesIndicesSize, attributesIndicesSize );

        m_classesCovarianceMatrixes.clear();
        m_classesCovarianceMatrixes.resize( samplesByClassSize, dummyCovarianceMatrix );

        m_classesCovarianceInvMatrixes.clear();
        m_classesCovarianceInvMatrixes.resize( samplesByClassSize, dummyCovarianceMatrix );

        m_classesOptizedMAPDiscriminantTerm.clear();
        m_classesOptizedMAPDiscriminantTerm.resize( samplesByClassSize );

        unsigned int attributeIdx1 = 0;
        unsigned int attributeIdx2 = 0;
        std::size_t classSamplesIdx = 0;
        std::size_t classSamplesSize = 0;
        double mean1 = 0;
        double mean2 = 0;
        double covariance = 0;

        for( std::size_t samplesByClassIdx = 0 ; samplesByClassIdx < samplesByClassSize ;
          ++samplesByClassIdx )
        {
          const std::vector< std::vector< double > >& classSamples = samplesByClass[
            samplesByClassIdx ];
          classSamplesSize = classSamples.size();

          for( attributeIdx1 = 0; attributeIdx1 < attributesIndicesSize; ++attributeIdx1 )
          {
            mean1 = m_classesMeans[ samplesByClassIdx ][ attributeIdx1 ];

            for( attributeIdx2 = 0; attributeIdx2 < attributesIndicesSize; ++attributeIdx2 )
            {
              mean2 = m_classesMeans[ samplesByClassIdx ][ attributeIdx2 ];
              covariance = 0.0;

              for( classSamplesIdx = 0 ; classSamplesIdx < classSamplesSize ; ++classSamplesIdx )
              {
                covariance +=
                  (
                    ( classSamples[ classSamplesIdx ][ attributeIdx1 ] - mean1 )
                    *
                    ( classSamples[ classSamplesIdx ][ attributeIdx2 ] - mean2 )
                  );
              }

              covariance /= (double)( classSamplesSize - 1 );

              m_classesCovarianceMatrixes[ samplesByClassIdx ]( attributeIdx1,
                attributeIdx2 ) = covariance;
            }
          }

          if( ! te::common::GetInverseMatrix( m_classesCovarianceMatrixes[ samplesByClassIdx ],
            m_classesCovarianceInvMatrixes[ samplesByClassIdx ] ) )
          {
            return false;
          }

          double classCovarianceMatrixDet = 0;
          if( ! te::common::GetDeterminant( m_classesCovarianceMatrixes[ samplesByClassIdx ],
            classCovarianceMatrixDet ) )
          {
            return false;
          }

          if( classCovarianceMatrixDet > 0.0 )
          {
            m_classesOptizedMAPDiscriminantTerm.push_back( -0.5 * std::log(
              classCovarianceMatrixDet ) );
          }
          else
          {
            m_classesOptizedMAPDiscriminantTerm.push_back( 0.0 );
          }
        }
      }

      if( enableProgressInterface )
      {
        progressPtr->pulse();
        if( ! progressPtr->isActive() ) return false;
      }

      return true;
    }

    bool MAP::classify(const InputAdaptor<double>& input,
                       const std::vector<unsigned int>& attributesIndices,
                       const std::vector<double>& inputNoDataValues,
                       OutputAdaptor<unsigned int>& output,
                       const unsigned int outputIndex,
                       const double outputNoDataValue,
                       const bool /*enableProgressInterface*/) throw(Exception)
    {
      if( !m_isInitialized ) return false;
      if( m_classesCovarianceMatrixes.empty() ) return false;
      if( input.getElementsCount() != output.getElementsCount() ) return false;
      if( attributesIndices.size() != m_classesMeans[ 0 ].size() ) return false;
      if( outputIndex >= output.getFeaturesCount() ) return false;
      if( inputNoDataValues.size() != attributesIndices.size() ) return false;

      // globals

      const unsigned int attributesIndicesSize = static_cast<unsigned int>(attributesIndices.size());
      const unsigned int inputCount = input.getElementsCount();
      const unsigned int inputFeaturesCount = input.getFeaturesCount();
      const unsigned int classesNmb = static_cast<unsigned int>(m_classesMeans.size());
      const unsigned int featuresCnt = static_cast<unsigned int>(m_classesMeans[0].size());

      // Dealing with the logarithm of priori probabilities

      std::vector< double > logPrioriProbs;

      if( m_parameters.m_prioriProbs.empty() )
      {
        if( ! getPrioriProbabilities( input,attributesIndices, logPrioriProbs ) )
        {
          return false;
        }

        for( unsigned int pIdx = 0 ; pIdx < logPrioriProbs.size() ; ++pIdx )
        {
          logPrioriProbs[ pIdx ] = ( logPrioriProbs[ pIdx ] > 0.0 ) ?
            std::log( logPrioriProbs[ pIdx ] ) : 0.0;
        }
      }
      else
      {
        for( unsigned int pIdx = 0 ; pIdx < m_parameters.m_prioriProbs.size() ;
          ++pIdx )
        {
          logPrioriProbs.push_back( std::log( m_parameters.m_prioriProbs[ pIdx ] ) );
        }
      }

      // classifying

      boost::numeric::ublas::matrix< double > sample( attributesIndicesSize, 1 );
      boost::numeric::ublas::matrix< double > sampleMinusMean( attributesIndicesSize, 1 );
      boost::numeric::ublas::matrix< double > sampleMinusMeanT( 1, attributesIndicesSize );
      boost::numeric::ublas::matrix< double > auxMatrix;
      boost::numeric::ublas::matrix< double > mahalanobisDistanceMatrix;
      unsigned int attributesIndicesIdx = 0;
      unsigned int featureIdx = 0;
      unsigned int classIdx = 0;
      double closestClassdiscriminantFunctionValue = 0;
      double discriminantFunctionValue = 0;
      unsigned int closestClassIdx = 0;
      bool isValidInput = false;

      for( unsigned int inputIdx = 0 ; inputIdx < inputCount ; ++inputIdx )
      {
        // load sample

        isValidInput = true;

        for( attributesIndicesIdx = 0 ; attributesIndicesIdx < attributesIndicesSize ;
          ++attributesIndicesIdx )
        {
          featureIdx = attributesIndices[ attributesIndicesIdx ];
          if( featureIdx >= inputFeaturesCount ) return false;

          input.getFeature( inputIdx, featureIdx, sample( attributesIndicesIdx, 0 ) );

          if( sample( attributesIndicesIdx, 0 ) == inputNoDataValues[ attributesIndicesIdx ] )
          {
            isValidInput = false;
            break;
          }
        }

        if( isValidInput )
        {
          // finding the clesest class

          closestClassdiscriminantFunctionValue = -1.0 * std::numeric_limits< double >::max();

          for( classIdx = 0 ; classIdx < classesNmb ; ++classIdx )
          {
            for( featureIdx = 0 ; featureIdx < featuresCnt ; ++featureIdx )
            {
              sampleMinusMean( featureIdx, 0 ) = sampleMinusMeanT( 0, featureIdx ) =
                ( sample( featureIdx, 0 ) - m_classesMeans[ classIdx ][ featureIdx ] );
            }

            //calculate the mahalanobis distance: (x-mean)T * CovMatrizInvert * (x-mean)

            auxMatrix = boost::numeric::ublas::prod( sampleMinusMeanT,
              m_classesCovarianceInvMatrixes[ classIdx ] );

            mahalanobisDistanceMatrix = boost::numeric::ublas::prod( auxMatrix, sampleMinusMean );
            assert( mahalanobisDistanceMatrix.size1() == 1 );
            assert( mahalanobisDistanceMatrix.size2() == 1 );

            discriminantFunctionValue = logPrioriProbs[ classIdx ]
              + m_classesOptizedMAPDiscriminantTerm[ classIdx ]
              - ( 0.5 * mahalanobisDistanceMatrix( 0, 0 ) );

            if( discriminantFunctionValue > closestClassdiscriminantFunctionValue )
            {
              closestClassdiscriminantFunctionValue = discriminantFunctionValue;
              closestClassIdx = classIdx;
            }
          }

          output.setFeature( inputIdx, 0, m_classLabels[ closestClassIdx ] );
        }
        else
        {
          output.setFeature(inputIdx, 0, static_cast<unsigned int>(outputNoDataValue));
        }
      }

      return true;
    }


    bool MAP::getPrioriProbabilities( const InputAdaptor< double >& input,
      const std::vector<unsigned int>& attributesIndices,
      std::vector< double >& prioriProbs ) const
    {
      const unsigned int inputCount = input.getElementsCount();
      const unsigned int inputFeaturesCount = input.getFeaturesCount();
      const unsigned int attributesIndicesSize = static_cast<unsigned int>(attributesIndices.size());
      const unsigned int classesNmb = static_cast<unsigned int>(m_classesMeans.size());
      const unsigned int featuresCnt = static_cast<unsigned int>(m_classesMeans[0].size());
      const double initialPrioriProbLog = std::log( 1.0 / ( (double)classesNmb ) );

      prioriProbs.resize( classesNmb );
      std::fill( prioriProbs.begin(), prioriProbs.end(), 0.0 );

      boost::numeric::ublas::matrix< double > sample( attributesIndicesSize, 1 );
      boost::numeric::ublas::matrix< double > sampleMinusMean( attributesIndicesSize, 1 );
      boost::numeric::ublas::matrix< double > sampleMinusMeanT( 1, attributesIndicesSize );
      boost::numeric::ublas::matrix< double > auxMatrix;
      boost::numeric::ublas::matrix< double > mahalanobisDistanceMatrix;
      unsigned int attributesIndicesIdx = 0;
      unsigned int featureIdx = 0;
      unsigned int classIdx = 0;
      double closestClassdiscriminantFunctionValue = 0;
      double discriminantFunctionValue = 0;
      unsigned int closestClassIdx = 0;
      unsigned int processedSamplesNmb = 0;

      for( unsigned int inputIdx = 0 ; inputIdx < inputCount ; inputIdx +=
        m_parameters.m_prioriCalcSampleStep )
      {
        // load sample

        for( attributesIndicesIdx = 0 ; attributesIndicesIdx < attributesIndicesSize ;
          ++attributesIndicesIdx )
        {
          featureIdx = attributesIndices[ attributesIndicesIdx ];
          if( featureIdx >= inputFeaturesCount ) return false;

          input.getFeature( inputIdx, featureIdx, sample( attributesIndicesIdx, 0 ) );
        }

        // finding the clesest class

        closestClassdiscriminantFunctionValue = -1.0 * std::numeric_limits< double >::max();

        for( classIdx = 0 ; classIdx < classesNmb ; ++classIdx )
        {
          for( featureIdx = 0 ; featureIdx < featuresCnt ; ++featureIdx )
          {
            sampleMinusMean( featureIdx, 0 ) = sampleMinusMeanT( 0, featureIdx ) =
              ( sample( featureIdx, 0 ) - m_classesMeans[ classIdx ][ featureIdx ] );
          }

          //calculate the mahalanobis distance: (x-mean)T * CovMatrizInvert * (x-mean)

          auxMatrix = boost::numeric::ublas::prod( sampleMinusMeanT,
            m_classesCovarianceInvMatrixes[ classIdx ] );

          mahalanobisDistanceMatrix = boost::numeric::ublas::prod( auxMatrix, sampleMinusMean );
          assert( mahalanobisDistanceMatrix.size1() == 1 );
          assert( mahalanobisDistanceMatrix.size2() == 1 );

          discriminantFunctionValue = initialPrioriProbLog
            + m_classesOptizedMAPDiscriminantTerm[ classIdx ]
            - ( 0.5 * mahalanobisDistanceMatrix( 0, 0 ) );

          if( discriminantFunctionValue > closestClassdiscriminantFunctionValue )
          {
            closestClassdiscriminantFunctionValue = discriminantFunctionValue;
            closestClassIdx = classIdx;
          }
        }

        prioriProbs[ closestClassIdx ] += 1.0;
        ++processedSamplesNmb;
      }

      for( classIdx = 0 ; classIdx < classesNmb ; ++classIdx )
      {
        prioriProbs[ classIdx ] /= ((double)processedSamplesNmb);
      }

      return true;
    }


  }  // end namespace cl
}    // end namespace te