libopenturns-dev 0.15-2 / usr / include / openturns / NumericalSampleImplementation.hxx

//                                               -*- C++ -*-
/**
 *  @file  NumericalSampleImplementation.hxx
 *  @brief The class NumericalSampleImplementation implements blank free samples
 *
 *  (C) Copyright 2005-2011 EDF-EADS-Phimeca
 *
 *  This library 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 2.1 of the License.
 *
 *  This library 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 this library; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
 *
 *  @author: $LastChangedBy: schueller $
 *  @date:   $LastChangedDate: 2011-05-24 19:30:41 +0200 (Tue, 24 May 2011) $
 *  Id:      $Id: NumericalSampleImplementation.hxx 1910 2011-05-24 17:30:41Z schueller $
 */
#ifndef OPENTURNS_NUMERICALSAMPLEIMPLEMENTATION_HXX
#define OPENTURNS_NUMERICALSAMPLEIMPLEMENTATION_HXX

#include <stdint.h> // for uint64_t

#include "NumericalPoint.hxx"
#include "Description.hxx"
#include "Indices.hxx"
#include "PersistentCollection.hxx"
#include "SquareMatrix.hxx"
#include "CovarianceMatrix.hxx"
#include "CorrelationMatrix.hxx"
#include "Exception.hxx"
#include "Collection.hxx"


namespace OpenTURNS
{

  namespace Base
  {

    namespace Stat
    {

#ifndef SWIG
      class NumericalSampleImplementation;
      class NSI_const_point;

      /****************************************/


      class NSI_point
      {
        NumericalSampleImplementation * p_nsi_;
        UnsignedLong index_;
        UnsignedLong dimension_;

        friend class NSI_const_point;
        friend std::ostream & operator << (std::ostream & os, const NSI_point & point);

      public:
        typedef       NumericalScalar *       iterator;
        typedef const NumericalScalar * const_iterator;

      public:
        NSI_point(NumericalSampleImplementation * p_nsi,
                  const UnsignedLong index);

        NSI_point & operator = (const NSI_point & rhs);
        NSI_point & operator = (const NSI_const_point & rhs);
        NSI_point & operator = (const Type::NumericalPoint & rhs);

        NumericalScalar & operator [] (const UnsignedLong i);
        const NumericalScalar & operator [] (const UnsignedLong i) const;
        NumericalScalar & at (const UnsignedLong i);
        const NumericalScalar & at (const UnsignedLong i) const;

        inline Type::Collection<NumericalScalar> getCollection() const { return Type::Collection<NumericalScalar>( begin(), end() ); }
        inline operator Type::NumericalPoint () const { return getCollection(); }

        inline UnsignedLong getDimension() const { return dimension_; }

        inline iterator begin() { return &operator[](0); }
        inline iterator end() { return &operator[](dimension_); }
        inline const_iterator begin() const { return &operator[](0); }
        inline const_iterator end() const { return &operator[](dimension_); }

        NSI_point & operator += (const NSI_point & other);
        NSI_point & operator -= (const NSI_point & other);
        NSI_point & operator += (const Type::NumericalPoint & other);
        NSI_point & operator -= (const Type::NumericalPoint & other);
        NSI_point & operator *= (const NumericalScalar val);
      };

      bool operator == (const NSI_point & lhs, const NSI_point & rhs);
      bool operator != (const NSI_point & lhs, const NSI_point & rhs);
      bool operator <  (const NSI_point & lhs, const NSI_point & rhs);
      bool operator >  (const NSI_point & lhs, const NSI_point & rhs);
      bool operator <= (const NSI_point & lhs, const NSI_point & rhs);
      bool operator >= (const NSI_point & lhs, const NSI_point & rhs);

      inline
      std::ostream & operator <<(std::ostream & os, const NSI_point & point)
      {
        return os << Type::NumericalPoint( point );
      }

      inline
      Common::OStream & operator << (Common::OStream & OS, const NSI_point & point)
      {
        OS.getStream() << Type::NumericalPoint( point ).__repr__();
        return OS;
      }


      class NSI_const_point
      {
        const NumericalSampleImplementation * p_nsi_;
        UnsignedLong index_;
        UnsignedLong dimension_;

        friend std::ostream & operator << (std::ostream & os, const NSI_const_point & point);

      public:
        typedef       NumericalScalar *       iterator;
        typedef const NumericalScalar * const_iterator;

      public:
        NSI_const_point(const NumericalSampleImplementation * p_nsi, UnsignedLong index);
        NSI_const_point(const NSI_point & point);

        const NumericalScalar & operator [] (const UnsignedLong i) const;
        const NumericalScalar & at (const UnsignedLong i) const;

        inline Type::Collection<NumericalScalar> getCollection() const { return Type::Collection<NumericalScalar>( begin(), end() ); }
        inline operator Type::NumericalPoint () const { return getCollection(); }

        inline UnsignedLong getDimension() const { return dimension_; }

        inline const_iterator begin() const { return &operator[](0); }
        inline const_iterator end() const { return &operator[](dimension_); }
      };

      bool operator == (const NSI_const_point & lhs, const NSI_const_point & rhs);
      bool operator != (const NSI_const_point & lhs, const NSI_const_point & rhs);
      bool operator <  (const NSI_const_point & lhs, const NSI_const_point & rhs);
      bool operator >  (const NSI_const_point & lhs, const NSI_const_point & rhs);
      bool operator <= (const NSI_const_point & lhs, const NSI_const_point & rhs);
      bool operator >= (const NSI_const_point & lhs, const NSI_const_point & rhs);

      inline
      std::ostream & operator <<(std::ostream & os, const NSI_const_point & point)
      {
        return os << Type::NumericalPoint( point );
      }

      inline
      Common::OStream & operator << (Common::OStream & OS, const NSI_const_point & point)
      {
        OS.getStream() << Type::NumericalPoint( point ).__repr__();
        return OS;
      }

      inline
      Type::NumericalPoint operator * (const NSI_const_point & point,
                                       const NumericalScalar val)
      {
        Type::NumericalPoint res(point.getDimension(), 0.);
        for(UnsignedLong i=0; i<point.getDimension(); ++i) res[i] = point[i] * val;
        return res;
      }

      inline
      Type::NumericalPoint operator * (const NumericalScalar val,
                                       const NSI_const_point & point)
      {
        return point * val;
      }




      /****************************************/


      class NSI_iterator
      {
        NumericalSampleImplementation * p_nsi_;
        UnsignedLong current_;

      public:
        typedef size_t difference_type;
        typedef std::random_access_iterator_tag iterator_category;
        typedef NSI_point value_type;
        typedef NSI_point * pointer;
        typedef NSI_point   reference;

        friend class NSI_const_iterator;
        friend bool operator == (const NSI_iterator & lhs, const NSI_iterator & rhs);
        friend bool operator <  (const NSI_iterator & lhs, const NSI_iterator & rhs);
        friend bool operator >  (const NSI_iterator & lhs, const NSI_iterator & rhs);
        friend difference_type operator - (const NSI_iterator & lhs, const NSI_iterator & rhs);

      public:
        NSI_iterator() : p_nsi_(0), current_(0) {}
        NSI_iterator(NumericalSampleImplementation & nsi, const UnsignedLong index) : p_nsi_(&nsi), current_(index) {}

        inline NSI_iterator & operator ++ () { ++current_; return *this; }
        inline NSI_iterator & operator -- () { --current_; return *this; }
        inline NSI_iterator   operator ++ (int) { NSI_iterator old(*this); ++current_; return old; }
        inline NSI_iterator   operator -- (int) { NSI_iterator old(*this); --current_; return old; }

        inline NSI_iterator & operator += (difference_type n) { current_ += n; return *this; }
        inline NSI_iterator & operator -= (difference_type n) { current_ -= n; return *this; }
        inline NSI_iterator   operator +  (difference_type n) const { NSI_iterator old(*this); old.current_ += n; return old; }
        inline NSI_iterator   operator -  (difference_type n) const { NSI_iterator old(*this); old.current_ -= n; return old; }

        inline NSI_point operator *  () const { return NSI_point(p_nsi_, current_); }
      };

      inline
      bool operator == (const NSI_iterator & lhs, const NSI_iterator & rhs)
      {
        return (lhs.p_nsi_ == rhs.p_nsi_) && (lhs.current_ == rhs.current_);
      }

      inline
      bool operator != (const NSI_iterator & lhs, const NSI_iterator & rhs)
      {
        return ! (lhs == rhs);
      }

      inline
      bool operator < (const NSI_iterator & lhs, const NSI_iterator & rhs)
      {
        return (lhs.p_nsi_ == rhs.p_nsi_) && (lhs.current_ < rhs.current_);
      }

      inline
      bool operator > (const NSI_iterator & lhs, const NSI_iterator & rhs)
      {
        return (lhs.p_nsi_ == rhs.p_nsi_) && (lhs.current_ > rhs.current_);
      }

      inline
      NSI_iterator::difference_type operator - (const NSI_iterator & lhs, const NSI_iterator & rhs)
      {
        return (lhs.current_ - rhs.current_);
      }




      class NSI_const_iterator
      {
        const NumericalSampleImplementation * p_nsi_;
        UnsignedLong current_;

      public:
        typedef size_t difference_type;
        typedef std::random_access_iterator_tag iterator_category;
        typedef NSI_const_point value_type;
        typedef NSI_const_point * pointer;
        typedef NSI_const_point   reference;

        friend bool operator == (const NSI_const_iterator & lhs, const NSI_const_iterator & rhs);
        friend bool operator <  (const NSI_const_iterator & lhs, const NSI_const_iterator & rhs);
        friend bool operator >  (const NSI_const_iterator & lhs, const NSI_const_iterator & rhs);
        friend difference_type operator - (const NSI_const_iterator & lhs, const NSI_const_iterator & rhs);

      public:
        NSI_const_iterator() : p_nsi_(0), current_(0) {}
        NSI_const_iterator(const NumericalSampleImplementation & nsi, const UnsignedLong index) : p_nsi_(&nsi), current_(index) {}
        NSI_const_iterator(const NSI_iterator & iter) : p_nsi_(iter.p_nsi_), current_(iter.current_) {}

        inline NSI_const_iterator & operator ++ () { ++current_; return *this; }
        inline NSI_const_iterator & operator -- () { --current_; return *this; }
        inline NSI_const_iterator   operator ++ (int) { NSI_const_iterator old(*this); ++current_; return old; }
        inline NSI_const_iterator   operator -- (int) { NSI_const_iterator old(*this); --current_; return old; }

        inline NSI_const_iterator & operator += (difference_type n) { current_ += n; return *this; }
        inline NSI_const_iterator & operator -= (difference_type n) { current_ -= n; return *this; }
        inline NSI_const_iterator   operator +  (difference_type n) const { NSI_const_iterator old(*this); old.current_ += n; return old; }
        inline NSI_const_iterator   operator -  (difference_type n) const { NSI_const_iterator old(*this); old.current_ -= n; return old; }

        inline NSI_const_point operator *  () const { return NSI_const_point(p_nsi_, current_); }
      };

      inline
      bool operator == (const NSI_const_iterator & lhs, const NSI_const_iterator & rhs)
      {
        return (lhs.p_nsi_ == rhs.p_nsi_) && (lhs.current_ == rhs.current_);
      }

      inline
      bool operator != (const NSI_const_iterator & lhs, const NSI_const_iterator & rhs)
      {
        return ! (lhs == rhs);
      }

      inline
      bool operator < (const NSI_const_iterator & lhs, const NSI_const_iterator & rhs)
      {
        return (lhs.p_nsi_ == rhs.p_nsi_) && (lhs.current_ < rhs.current_);
      }

      inline
      bool operator > (const NSI_const_iterator & lhs, const NSI_const_iterator & rhs)
      {
        return (lhs.p_nsi_ == rhs.p_nsi_) && (lhs.current_ > rhs.current_);
      }

      inline
      NSI_const_iterator::difference_type operator - (const NSI_const_iterator & lhs, const NSI_const_iterator & rhs)
      {
        return (lhs.current_ - rhs.current_);
      }

#endif


      /****************************************/



      /**
       * @class NumericalSampleImplementation
       */

      class NumericalSampleImplementation
        : public Common::PersistentObject
      {
        CLASSNAME;
        friend class NSI_point;
        friend class NSI_const_point;
        friend class NSI_iterator;
        friend class NSI_const_iterator;
        friend class Common::Factory<NumericalSampleImplementation>;
        friend class BuildMethodMap;
        friend class ExportMethodMap;
#ifndef SWIG
        friend Bool operator ==(const NumericalSampleImplementation & lhs, const NumericalSampleImplementation & rhs);
#endif

      public:

        /* Some typedefs for easy reading */
        typedef NSI_iterator               iterator;
        typedef NSI_const_iterator         const_iterator;

        typedef Type::NumericalPoint             NumericalPoint;
        typedef Type::Description                Description;
        typedef Type::Indices                    Indices;
        typedef Type::Collection<UnsignedLong>   UnsignedLongCollection;
        typedef Common::FileNotFoundException    FileNotFoundException;
        typedef Common::InvalidArgumentException InvalidArgumentException;
        typedef Common::InternalException        InternalException;
        typedef Stat::CovarianceMatrix           CovarianceMatrix;
        typedef Stat::CorrelationMatrix          CorrelationMatrix;
        typedef Type::SquareMatrix               SquareMatrix;
        typedef Common::StorageManager           StorageManager;
        typedef Common::PersistentObject         PersistentObject;

        enum ExternalFileFormat { CSV };

      public:

        /** Constructor from file */
        static NumericalSampleImplementation GetFromFile(const FileName & fileName, const ExternalFileFormat format = CSV);

        /** Factory of NumericalSampleImplementation from CSV file */
        static NumericalSampleImplementation BuildFromCSVFile(const FileName & fileName) /* throw(FileNotFoundException, InternalException) */;

        /** Store a sample in a temporary text file, one realization by line. Returns the file name. */
        String storeToTemporaryFile() const;

        /** Export a sample as a matrix, one row by realization, in a format suitable to exchange with R */
        String streamToRFormat() const;

      protected:

        /**
         * Default constructor is protected
         */
        NumericalSampleImplementation();

      public:

        /** Standard constructor */
        NumericalSampleImplementation(const UnsignedLong size,
                                      const UnsignedLong dim);

        /** Constructor from a NumericalPoint */
        NumericalSampleImplementation(const UnsignedLong size,
                                      const NumericalPoint & point);

#ifndef SWIG
        /** Constructor from a collection of NumericalPoint */
        NumericalSampleImplementation(const Type::Collection<NumericalPoint> & coll);

        /** Partial copy constructor */
        NumericalSampleImplementation(const NumericalSampleImplementation & other,
                                      iterator first,
                                      iterator last);

        inline iterator begin() { return iterator(*this, 0); }
        inline iterator end() { return iterator(*this, size_); }
        inline const_iterator begin() const { return const_iterator(*this, 0); }
        inline const_iterator end() const { return const_iterator(*this, size_); }

        void erase(const UnsignedLong first, const UnsignedLong last);
        void erase(iterator first, iterator last);
        void clear();

        inline NSI_point operator [] (const UnsignedLong index) { return NSI_point(this, index); }
        inline NSI_const_point operator [] (const UnsignedLong index) const { return NSI_const_point(this, index); }
        inline NumericalScalar & atLinearIndex (const UnsignedLong index) { return data_[index]; }
        inline const NumericalScalar & atLinearIndex (const UnsignedLong index) const { return data_[index]; }

        void swap_points(const UnsignedLong a, const UnsignedLong b);
        void swap_range_points(const UnsignedLong fa, const UnsignedLong ta, const UnsignedLong fb);
#endif

        /** Virtual constructor */
        NumericalSampleImplementation * clone() const;

        /** Description Accessor */
        void setDescription(const Description & description);
        Description getDescription() const;

        /**
         * String converter
         * This method shows human readable information on the
         * internal state of an NumericalSampleImplementation. It is used when streaming
         * the NumericalSampleImplementation or for user information.
         */
        virtual String __repr__() const;
        virtual String __str__(const String & offset = "") const;

        inline Bool __eq__(const NumericalSampleImplementation & rhs) const
        { return (*this == rhs); }

        /* Method contains() is for Python */
        Bool contains(const NumericalPoint & val) const;

        /** Size accessor */
        inline UnsignedLong getSize() const { return size_; }

        /** Dimension accessor */
        inline UnsignedLong getDimension() const { return dimension_; }

        /** Appends an element to the collection */
        NumericalSampleImplementation & add(const NumericalPoint & point) /* throw(InvalidArgumentException) */;

        /** Appends another sample to the collection */
        NumericalSampleImplementation & add(const NumericalSampleImplementation & sample) /* throw(InvalidArgumentException) */;

        /**
         * Gives the mean of the sample, based on the formula
         * mean = sum of the elements in the sample / size of the sample
         */
        virtual NumericalPoint computeMean() const;

        /**
         * Gives the covariance matrix of the sample, normalization by 1 / (size - 1) if size > 1
         */
        virtual CovarianceMatrix computeCovariance() const;

        /**
         * Gives the standard deviation of the sample, i.e. the square-root of the covariance matrix.
         */
        SquareMatrix computeStandardDeviation() const;

        /**
         * Gives the standard deviation of each component of the sample
         */
        virtual NumericalPoint computeStandardDeviationPerComponent() const;

        /**
         * Gives the Pearson correlation matrix of the sample
         */
        virtual CorrelationMatrix computePearsonCorrelation() const;

        /**
         * Gives the Spearman correlation matrix of the sample
         */
        virtual CorrelationMatrix computeSpearmanCorrelation() const;

        /**
         * Gives the Kendall correlation matrix of the sample
         */
        virtual CorrelationMatrix computeKendallTau() const;

        /**
         * Gives the range of the sample (by component)
         */
        NumericalPoint computeRangePerComponent() const;

        /**
         * Gives the median of the sample (by component)
         */
        NumericalPoint computeMedianPerComponent() const;

        /**
         * Gives the variance of the sample (by component)
         */
        virtual NumericalPoint computeVariancePerComponent() const;

        /**
         * Gives the skewness of the sample (by component)
         */
        virtual NumericalPoint computeSkewnessPerComponent() const;

        /**
         * Gives the kurtosis of the sample (by component)
         */
        virtual NumericalPoint computeKurtosisPerComponent() const;

        /**
         * Gives the centered moment of order k of the sample (by component)
         */
        NumericalPoint computeCenteredMomentPerComponent(const UnsignedLong k) const;

        /**
         * Gives the quantile per component of the sample
         */
        virtual NumericalPoint computeQuantilePerComponent(const NumericalScalar prob) const;

        /**
         * Gives the N-dimension quantile of the sample
         */
        NumericalPoint computeQuantile(const NumericalScalar prob) const;

        /**
         * Get the empirical CDF of the sample
         */
        NumericalScalar computeEmpiricalCDF(const NumericalPoint & point,
                                            const Bool tail = false) const;

        /** Maximum accessor */
        virtual NumericalPoint getMax() const;

        /** Minimum accessor */
        virtual NumericalPoint getMin() const;

        /** Ranked sample */
        NumericalSampleImplementation rank() const;

        /** Ranked component */
        NumericalSampleImplementation rank(const UnsignedLong index) const;

        /** Sorted sample */
        NumericalSampleImplementation sort() const;

        /** Sorted component */
        NumericalSampleImplementation sort(const UnsignedLong index) const;

        /** Sorted component */
        NumericalSampleImplementation sortAccordingToAComponent(const UnsignedLong index) const;

        /** Get the i-th marginal sample */
        NumericalSampleImplementation getMarginal(const UnsignedLong index) const /* throw(InvalidArgumentException) */;

        /** Get the marginal sample corresponding to indices dimensions */
        NumericalSampleImplementation getMarginal(const Indices & indices) const /* throw(InvalidArgumentException) */;

        /**
         * Translate realizations in-place
         */
        void translate(const NumericalPoint & translation);

        /**
         * Scale realizations componentwise in-place
         */
        void scale(const NumericalPoint & scaling);

        /** Save to CSV file */
        void exportToCSVFile(const FileName & filename,
                             const Bool withDescription = false) const;

        /** Method save() stores the object through the StorageManager */
        virtual void save(StorageManager::Advocate & adv) const;

        /** Method load() reloads the object from the StorageManager */
        virtual void load(StorageManager::Advocate & adv);


      private:

        /** Merge sort counting the swaps, intermediate for the fast implementation of computeKendallTau() */
        uint64_t mergeSort(const UnsignedLong offset,
                           const UnsignedLong length,
                           const UnsignedLong activeDimension,
                           UnsignedLongCollection & ordering,
                           UnsignedLongCollection & buffer) const;

        /** Lexical ordering between two dimensions: ranking with respect to i then to j if ties */
        UnsignedLongCollection lexicalRank(const UnsignedLong i,
                                           const UnsignedLong j) const;

        /** The size of the sample */
        UnsignedLong size_;

        /** The dimension of the sample */
        UnsignedLong dimension_;

        /** The collection of unwrapped points */
        Type::PersistentCollection<NumericalScalar> data_;

        /** The description of all components */
        Description::Implementation p_description_;

      }; /* class NumericalSampleImplementation */


      /* Comparison function */
      Bool operator ==(const NumericalSampleImplementation & lhs, const NumericalSampleImplementation & rhs);


    } /* namespace Stat */
  } /* namespace Base */
} /* namespace OpenTURNS */

#endif /* OPENTURNS_NUMERICALSAMPLEIMPLEMENTATION_HXX */