Logo Search packages:      
Sourcecode: openturns version File versions  Download package

NumericalSampleImplementation.cxx

Go to the documentation of this file.
//                                               -*- C++ -*-
/**
 *  @file  NumericalSampleImplementation.cxx
 *  @brief The class NumericalSampleImplementation implements blank free samples
 *
 *  (C) Copyright 2005-2007 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: dutka $
 *  @date:   $LastChangedDate: 2008-10-31 11:52:04 +0100 (ven 31 oct 2008) $
 *  Id:      $Id: NumericalSampleImplementation.cxx 995 2008-10-31 10:52:04Z dutka $
 */
#include <map>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <algorithm>
#include <iterator>
#include <algorithm>
#include <cstdlib>
#include <cstdio>        // std::fopen
#include <cstdio>        // std::errno
#include <cstring>       // std::strerror

#include "NumericalSampleImplementation.hxx"
#include "StorageManager.hxx"
#include "PersistentObjectFactory.hxx"
#include "Log.hxx"
#include "Exception.hxx"
#include "ResourceMap.hxx"
#include "Path.hxx"

#include "csv_parser_state.hxx"
#include "csv_parser.hh"
#include "csv_lexer.h"

int csvparse (OT::Base::Stat::CSVParserState & theState, yyscan_t yyscanner, FILE * theFile, OT::Base::Stat::NumericalSampleImplementation &impl, OT::UnsignedLong & theDimension, const char * separator);


namespace OpenTURNS
{

  namespace Base
  {

    namespace Type
    {
      TEMPLATE_CLASSNAMEINIT(PersistentCollection<NumericalPoint>);

      static Common::Factory<PersistentCollection<NumericalPoint> > RegisteredFactory("PersistentCollection<NumericalPoint>");
    } /* namespace Type */



    namespace Stat
    {

      typedef NumericalSampleImplementation (*BuildMethod) (const FileName & fileName);
      using Common::NotYetImplementedException;
      using Common::OutOfBoundException;
      using Common::Log;
      using Common::ResourceMap;
      using Common::Path;

      static Common::Factory<NumericalSampleImplementation> RegisteredFactory("NumericalSampleImplementation");


      /*
       * This class implements a map that stores a function pointer to a factory that builds
       * a NumericalSampleImplementation according to a file format
       */
00085       class BuildMethodMap : public std::map<NumericalSampleImplementation::ExternalFileFormat, BuildMethod>
      {
      typedef std::map<NumericalSampleImplementation::ExternalFileFormat, BuildMethod> ParentType;
      ParentType & table_;
      
      public:
      BuildMethodMap()
        : std::map<NumericalSampleImplementation::ExternalFileFormat, BuildMethod>(),
          table_(*this)
      {
        table_[NumericalSampleImplementation::CSV] = NumericalSampleImplementation::BuildFromCSVFile;
      }

      const BuildMethod & operator[] (const NumericalSampleImplementation::ExternalFileFormat & format) const
      {
        return table_[format];
      }

      }; /* end class BuildMethodMap */


      static const BuildMethodMap SampleImportationFactoryMap;



      CLASSNAMEINIT(NumericalSampleImplementation);


      /* Constructor from file */
00114       NumericalSampleImplementation NumericalSampleImplementation::GetFromFile(const FileName & fileName,
                                                             const ExternalFileFormat format)
      {
      return SampleImportationFactoryMap[format](fileName);
      }

      /* Factory of NumericalSampleImplementation from CSV file */
00121       NumericalSampleImplementation NumericalSampleImplementation::BuildFromCSVFile(const FileName & fileName) throw(FileNotFoundException, InternalException)
      {
      String csvSeparator = ResourceMap::GetInstance().get("csv-file-separator");

      yyscan_t scanner = 0;
      NumericalSampleImplementation impl(0, 0);
      impl.setName(fileName);

      FILE * theFile = std::fopen(fileName.c_str(),"r");
      if (!theFile) { // theFile can not be found. Errno is set
        throw FileNotFoundException(HERE) << "Can NOT open file '" << fileName
                                  << "'. Reason: " << std::strerror(errno);
      }
      CSVParserState state;
      state.theFileName = fileName;

      csvlex_init(&scanner);
      csvparse(state, scanner, theFile, impl, impl.dimension_, csvSeparator.c_str());
      csvlex_destroy(scanner);
      std::fclose(theFile);

      return impl;
      }

      /* Store a sample in a temporary text file, one realization by line. Returns the file name. */
00146       String NumericalSampleImplementation::storeToTemporaryFile() const
      {
      String dataFileName(Path::BuildTemporaryFileName("RData.txt.XXXXXX"));
        std::ofstream dataFile(dataFileName.c_str());
      // Fill-in the data file
      UnsignedLong size(getSize());
      UnsignedLong dimension(getDimension());
      for (UnsignedLong i = 0; i < size; ++i)
        {
          String separator = "";
          for (UnsignedLong j = 0; j < dimension; ++j, separator = " ")
            {
            dataFile << separator << std::setprecision(20) << operator[](i)[j];
            }
          dataFile << std::endl;
        }
      dataFile.close();
      return dataFileName;
      }

      /* Export a sample as a matrix, one row by realization, in a format suitable to exchange with R */
00167       String NumericalSampleImplementation::streamToRFormat() const
      {
      UnsignedLong size(getSize());
      UnsignedLong dimension(getDimension());
      OSS oss(20);
      oss << "matrix(c(";
      String separator("");
      for (UnsignedLong j = 0; j < dimension; ++j)
        {
          for (UnsignedLong i = 0; i < size; ++i, separator = ",")
            {
            oss << separator << operator[](i)[j];
            }
        }
      oss << "), nrow=" << size << ", ncol=" << dimension << ")";
      return oss;
      }

      /* Default constructor is private */
00186       NumericalSampleImplementation::NumericalSampleImplementation()
      : Type::PersistentCollection<NumericalSampleImplementation::NumericalPoint>(),
        dimension_(),
        p_description_()
      {
        // Nothing to do
      }

      /* Standard constructor */
00195       NumericalSampleImplementation::NumericalSampleImplementation(UnsignedLong size,
                                                   UnsignedLong dim)
      : Type::PersistentCollection<NumericalSampleImplementation::NumericalPoint>(size, NumericalSampleImplementation::NumericalPoint(dim, 0.)),
        dimension_(dim),
        p_description_()
      {
        // Nothing to do
      }

      /* Constructor from a NumericalPoint */
00205       NumericalSampleImplementation::NumericalSampleImplementation(const UnsignedLong size,
                                                   const NumericalPoint & point)
      : Type::PersistentCollection<NumericalSampleImplementation::NumericalPoint>(size, point),
        dimension_(point.getDimension()),
        p_description_()
      {
      // Nothing to do
      }

      /* Partial copy constructor */
00215       NumericalSampleImplementation::NumericalSampleImplementation(const NumericalSampleImplementation & other, iterator first, iterator last)
      : Type::PersistentCollection<NumericalSampleImplementation::NumericalPoint>(first, last),
        dimension_(other.dimension_),
        p_description_(other.p_description_)
      {
        // Nothing to do
      }


      /* Virtual constructor */
00225       NumericalSampleImplementation * NumericalSampleImplementation::clone() const
      {
      return new NumericalSampleImplementation(*this);
      }



      /* Description Accessor */
00233       void NumericalSampleImplementation::setDescription(const Description & description)
      {
      p_description_ = description.getImplementation();
      }


      /* Description Accessor */
      NumericalSampleImplementation::Description NumericalSampleImplementation::getDescription() const
      {
      return p_description_.isNull() ? Description(getDimension()) : *p_description_;
      }

      /* Comparison operator */
00246       Bool NumericalSampleImplementation::operator ==(const NumericalSampleImplementation & other) const
      {
      return true;
      }

      /* String converter */
00252       String NumericalSampleImplementation::str() const
      {
      return OSS() << "class=" << NumericalSampleImplementation::GetClassName()
                 << " name=" << getName()
                 << " size=" << getSize()
                 << " dimension=" << dimension_
                 << " data=" << Type::PersistentCollection<NumericalSampleImplementation::NumericalPoint>::str();
      }



      /* Dimension accessor */
00264       UnsignedLong NumericalSampleImplementation::getDimension() const
      {
      return dimension_;
      }


      /* Appends an element to the collection */
00271       void NumericalSampleImplementation::add(const NumericalPoint & point)
      throw(InvalidArgumentException)
      {
      if ( (dimension_ != 0) && (point.getDimension() != dimension_) )
        throw InvalidArgumentException(HERE) << "Point has invalid dimension ("
                                     << point.getDimension()
                                     << ") expected : "
                                     << getDimension();
      NumericalPoint newPoint = point;
      Type::PersistentCollection<NumericalSampleImplementation::NumericalPoint>::add(newPoint);
      }

      /* 
       * Gives the mean of the sample, based on the formula
       * mean = sum of the elements in the sample / size of the sample
       */
00287       NumericalSampleImplementation::NumericalPoint NumericalSampleImplementation::computeMean() const
      {
      const UnsignedLong size = getSize();
      NumericalPoint mean(getDimension(), 0.0);
      for (UnsignedLong index = 0; index < size; ++index)
        mean += (*this)[index];
      return mean * (1.0 / static_cast<NumericalScalar>(size));
      }

      /* 
       * Gives the covariance matrix of the sample, normalization by 1 / (size - 1) if size > 1
       */
00299       NumericalSampleImplementation::CovarianceMatrix NumericalSampleImplementation::computeCovariance() const
      {
      const UnsignedLong dimension = getDimension();
      const UnsignedLong size = getSize();

      // Special case for a sample of size 1
      if (getSize() == 1) return CovarianceMatrix(dimension);
      NumericalPoint mean(computeMean());
      CovarianceMatrix covariance(dimension);
      for (UnsignedLong index = 0; index < size; ++index) {
        NumericalPoint realization( (*this)[index] );
        for (UnsignedLong i = 0; i < dimension; ++i)
          for (UnsignedLong j = 0; j <= i; ++j)
            covariance(i, j) += (realization[i] - mean[i]) * (realization[j] - mean[j]);
      }
      const NumericalScalar alpha = 1.0 / (static_cast<NumericalScalar>(size) - 1.0);
      for (UnsignedLong i = 0; i < dimension; ++i)
        for (UnsignedLong j = 0; j <= i; ++j)
          covariance(i, j) *= alpha;
      return covariance;
      }

      /*
       * Gives the standard deviation of the sample, i.e. the square-root of the covariance matrix.
       */
00324       NumericalSampleImplementation::SquareMatrix NumericalSampleImplementation::computeStandardDeviation() const
      {
      return computeCovariance().computeCholesky();
      }

      /*
       * Gives the standard deviation of each component of the sample
       */
00332       NumericalSampleImplementation::NumericalPoint NumericalSampleImplementation::computeStandardDeviationPerComponent() const
      {
      const UnsignedLong dimension = getDimension();
      const UnsignedLong size = getSize();

      // Special case for a sample of size 1
      if (size == 1) return NumericalPoint(dimension, 0.0);
      NumericalPoint mean(computeMean());
      NumericalPoint sd(dimension);
      for (UnsignedLong index = 0; index < size; ++index) {
        NumericalPoint centeredRealization((*this)[index] - mean);
        for (UnsignedLong i = 0; i < dimension; ++i)
          sd[i] += centeredRealization[i] * centeredRealization[i];
      }
      sd *= 1.0 / (static_cast<NumericalScalar>(size) - 1.0);
      for (UnsignedLong i = 0; i < dimension; ++i)
        sd[i] = sqrt(sd[i]);
      return sd;
      }



      /*
       * Gives the Pearson correlation matrix of the sample
       */
00357       NumericalSampleImplementation::CorrelationMatrix NumericalSampleImplementation::computePearsonCorrelation() const
      {
      const UnsignedLong dimension = getDimension();
      CorrelationMatrix correlation(dimension);
      if (dimension == 1) return correlation;
      CovarianceMatrix covariance(computeCovariance());
      NumericalPoint sd(dimension);
      for (UnsignedLong i = 0; i < dimension; ++i) {
        sd[i] = sqrt( covariance(i,i) );
        for (UnsignedLong j = 0; j < i; ++j)
          correlation(i, j) = covariance(i, j) / (sd[i] * sd[j]);
      }
      return correlation;
      }

00372       struct Pair
      {
        NumericalScalar value_;
        UnsignedLong index_;
        Pair():value_(0.0), index_() {}
        Pair(const NumericalScalar value, const UnsignedLong index): value_(value), index_(index) {}
        ~Pair() {}
        Bool operator < (const Pair & other) const { return value_ < other.value_; }
      };

      typedef Type::Collection<Pair>                   PairCollection;
      typedef Type::Collection<PairCollection>         PairCollectionCollection;
      typedef Type::Collection<UnsignedLong>           UnsignedLongCollection;
      typedef Type::Collection<UnsignedLongCollection> UnsignedLongCollectionCollection;

      /* Ranked sample */
00388       NumericalSampleImplementation NumericalSampleImplementation::rank() const
      {
      const UnsignedLong dimension = getDimension();
      const UnsignedLong size = getSize();

      NumericalSampleImplementation rankedSample(size, dimension);
      PairCollectionCollection sortedMarginalSamples(dimension, PairCollection(size));

      // Sort and rank all the marginal samples
      for (UnsignedLong i = 0; i < dimension; ++i) {
        for (UnsignedLong j = 0; j < size; ++j) {
          sortedMarginalSamples[i][j].value_ = (*this)[j][i];
          sortedMarginalSamples[i][j].index_ = j;
        }
        // sort
        std::stable_sort(sortedMarginalSamples[i].begin(), sortedMarginalSamples[i].end());
        // rank
        for (UnsignedLong j = 0; j < size; ++j)
          rankedSample[ sortedMarginalSamples[i][j].index_ ][i] = j;
      }
      return rankedSample;
      }

      /* Ranked component */
00412       NumericalSampleImplementation NumericalSampleImplementation::rank(const UnsignedLong index) const
      {
      if (index >= getDimension()) throw OutOfBoundException(HERE) << "Error: the requested index is too large, index=" << index << ", dimension=" << getDimension();
      const UnsignedLong size = getSize();
      NumericalSampleImplementation rankedSample(size, 1);
      PairCollectionCollection sortedMarginalSamples(1, PairCollection(size));

      // Sort and rank the marginal sample number index
      for (UnsignedLong j = 0; j < size; ++j) {
        sortedMarginalSamples[0][j].value_ = (*this)[j][index];
        sortedMarginalSamples[0][j].index_ = j;
      }
      // sort
      std::stable_sort(sortedMarginalSamples[0].begin(), sortedMarginalSamples[0].end());
      // rank
      for (UnsignedLong j = 0; j < size; ++j)
        rankedSample[ sortedMarginalSamples[0][j].index_ ][0] = j;
      return rankedSample;
      }

      /* Sorted sample, component by component */
00433       NumericalSampleImplementation NumericalSampleImplementation::sort() const
      {
      const UnsignedLong dimension = getDimension();
      const UnsignedLong size = getSize();

      NumericalSampleImplementation sortedSample(size, dimension);
      Type::Collection<NumericalScalar> component(size);
      // Sort all the marginal samples
      for (UnsignedLong i = 0; i < dimension; ++i) {
        for (UnsignedLong j = 0; j < size; ++j)
          component[j] = (*this)[j][i];
        // sort
        std::stable_sort(component.begin(), component.end());
        // copy
        for (UnsignedLong j = 0; j < size; ++j)
          sortedSample[j][i] = component[j];
      }
      return sortedSample;
      }

      /* Sorted sample, one component */
00454       NumericalSampleImplementation NumericalSampleImplementation::sort(const UnsignedLong index) const
      {
      if (index >= getDimension()) throw OutOfBoundException(HERE) << "Error: the requested index is too large, index=" << index << ", dimension=" << getDimension();
      const UnsignedLong size = getSize();
      NumericalSampleImplementation sortedSample(size, 1);
      Type::Collection<NumericalScalar> component(size);
      // Sort the requested component
      for (UnsignedLong j = 0; j < size; ++j)
        component[j] = (*this)[j][index];
      // sort
      std::stable_sort(component.begin(), component.end());
      // copy
      for (UnsignedLong j = 0; j < size; ++j)
        sortedSample[j][0] = component[j];
      return sortedSample;
      }

      /* Sorted according a component */
00472       NumericalSampleImplementation NumericalSampleImplementation::sortAccordingAComponent(const UnsignedLong index) const
      {
      NumericalSampleImplementation rankedIndex(rank(index));
      const UnsignedLong size = getSize();
      NumericalSampleImplementation result(size, getDimension());
      for (UnsignedLong i = 0; i < size; ++i)
        result[static_cast<UnsignedLong>( round(rankedIndex[i][0]) ) ] = (*this)[i];
      return result;
      }

      /*
       * Gives the Spearman correlation matrix of the sample
       */
00485       NumericalSampleImplementation::CorrelationMatrix NumericalSampleImplementation::computeSpearmanCorrelation() const
      {
      return rank().computePearsonCorrelation();
      }

      /*
       * Gives the Kendall tau matrix of the sample
       */
00493       NumericalSampleImplementation::CorrelationMatrix NumericalSampleImplementation::computeKendallTau() const
      {
      const UnsignedLong dimension = getDimension();
      const UnsignedLong size = getSize();

      CorrelationMatrix tau(dimension);
      // For the first component of all pairs
      for (UnsignedLong i = 1; i < size; i++)
        {
          NumericalPoint xI(operator[](i));
          // For the second component of all pairs
          for (UnsignedLong j = 0; j < i; j++)
            {
            NumericalPoint xJ(operator[](j));
            // For each first component
            for (UnsignedLong m = 0; m < dimension; m++)
              {
                // For each second component
                for (UnsignedLong n = 0; n < m; n++)
                  {
                  // Contribution is +1 if the two components are strictly ordered the same way, 0 if there is a tie, -1 otherwise.
                  tau(m, n) += (xI[m] < xJ[m] ? -1.0 : xI[m] > xJ[m] ? 1.0 : 0.0) * (xI[n] < xJ[n] ? -1.0 : xI[n] > xJ[n] ? 1.0 : 0.0);
                  }
              }
            }
        }
      // Normalization
      NumericalScalar normalizationFactor = 2.0 / (static_cast<NumericalScalar>(size) * (static_cast<NumericalScalar>(size) + 1.0));
      // For each first component
      for (UnsignedLong m = 1; m < dimension; ++m)
        // For each second component
        for (UnsignedLong n = 0; n < m; ++n)
          tau(m, n) *= normalizationFactor;
      return tau;
      }

      /*
       * Gives the range of the sample (by component)
       */
00532       NumericalSampleImplementation::NumericalPoint NumericalSampleImplementation::computeRangePerComponent() const
      {
      return getMax() - getMin();
      }

      /*
       * Gives the median of the sample (by component)
       */
00540       NumericalSampleImplementation::NumericalPoint NumericalSampleImplementation::computeMedianPerComponent() const
      {
      return computeQuantilePerComponent(0.5);
      }
      
      /*
       * Gives the variance of the sample (by component)
       */
00548       NumericalSampleImplementation::NumericalPoint NumericalSampleImplementation::computeVariancePerComponent() const
      {
      const UnsignedLong dimension = getDimension();
      const UnsignedLong size = getSize();

      // Special case for a sample of size 1
      if (size == 1) return NumericalPoint(dimension, 0.0);
      NumericalPoint mean(computeMean());
      NumericalPoint var(dimension);
      for (UnsignedLong index = 0; index < size; ++index) {
        NumericalPoint centeredRealization((*this)[index] - mean);
        for (UnsignedLong i = 0; i < dimension; ++i)
          var[i] += centeredRealization[i] * centeredRealization[i];
      }
      var *= 1.0 / (static_cast<NumericalScalar>(size) - 1.0);
      return var;
      }        

      /*
       * Gives the skewness of the sample (by component)
       */
00569       NumericalSampleImplementation::NumericalPoint NumericalSampleImplementation::computeSkewnessPerComponent() const
      {
      const UnsignedLong dimension = getDimension();
      const UnsignedLong size = getSize();

      // Special case for a sample of size 1
      if (size == 1) return NumericalPoint(dimension, 0.0);
      NumericalPoint mean(computeMean());
      NumericalPoint skewness(dimension);
      NumericalPoint var(dimension);
      for (UnsignedLong index = 0; index < size; index++)
        {
          NumericalPoint centeredRealization(operator[](index) - mean);
          for (UnsignedLong i = 0; i < dimension; i++)
            {
            // var = sum (Xi - Xmean)^2
            NumericalScalar square(centeredRealization[i] * centeredRealization[i]);
            var[i] += square;
            // skewness = sum (Xi - Xmean)^3
            skewness[i] += square * centeredRealization[i];
            }
        }
      
      const NumericalScalar factor = 1.0 / static_cast<NumericalScalar>(size);
      const NumericalScalar factor1 = 1.0 / (static_cast<NumericalScalar>(size) - 1.0);
      for (UnsignedLong i = 0; i < dimension; ++i) {
        // skewness = 1 / size sum (Xi - Xmean)^3 / (var/(size-1))^3/2
        skewness[i] *= factor * pow(var[i] * factor1, -1.5);
      }
      return skewness;
      }           

      /*
       * Gives the kurtosis of the sample (by component)
       */
00604       NumericalSampleImplementation::NumericalPoint NumericalSampleImplementation::computeKurtosisPerComponent() const
      {
      const UnsignedLong dimension = getDimension();
      const UnsignedLong size = getSize();

      // Special case for a sample of size 1
      if (size == 1) return NumericalPoint(dimension, 0.0);
      NumericalPoint mean(computeMean());
      NumericalPoint kurtosis(dimension);
      NumericalPoint var(dimension);
      for (UnsignedLong index = 0; index < size; index++)
        {
          NumericalPoint centeredRealization(operator[](index) - mean);
          for (UnsignedLong i = 0; i < dimension; i++)
            {
            // var = sum (Xi - Xmean)^2
            NumericalScalar square(centeredRealization[i] * centeredRealization[i]);
            var[i] += square;
            // kurtosis = sum (Xi - Xmean)^4
            kurtosis[i] += square * square;
            }
        }
      
      const NumericalScalar factor = 1.0 / static_cast<NumericalScalar>(size);
      const NumericalScalar factor1 = 1.0 / (static_cast<NumericalScalar>(size) - 1.0);
      for (UnsignedLong i = 0; i < dimension; ++i) {
        // kurtosis = 1 / size sum (Xi - Xmean)^4 / (var/(size-1))^2 - 3
        kurtosis[i] *= factor * pow(var[i] * factor1, -2);
        kurtosis[i] -= 3.0;
      }
      return kurtosis;
      }     

      /*
       * Gives the quantile per component of the sample
       */
00640       NumericalSampleImplementation::NumericalPoint NumericalSampleImplementation::computeQuantilePerComponent(const NumericalScalar prob) const
      {
      // Special case for prob >= 1
      if (prob >= 1.0) return getMax();
      // Special case for prob <= 0.0
      if (prob <= 0.0) return getMin();

      const UnsignedLong dimension = getDimension();
      const UnsignedLong size = getSize();
      const NumericalScalar scalarIndex = prob * static_cast<NumericalScalar>(size);
      const UnsignedLong index    = static_cast<UnsignedLong>( floor( scalarIndex ) );
      const NumericalScalar beta  = scalarIndex - index;
      const NumericalScalar alpha = 1.0 - beta;
      NumericalPoint quantile(dimension);
      NumericalPoint component(size);
      for (UnsignedLong j = 0; j < dimension; ++j) {
        for (UnsignedLong i = 0; i < size; ++i)
          component[i] = operator[](i)[j];
        std::sort(component.begin(), component.end());
        if (index == size - 1)
          quantile[j] = component[index];
        else
          // Interpolation between the two adjacent empirical quantiles
          quantile[j] = alpha * component[index] + beta * component[index + 1];
      }
      return quantile;
      }

      /*
       * Gives the N-dimension quantile of the sample
       */
00671       NumericalSampleImplementation::NumericalPoint NumericalSampleImplementation::computeQuantile(const NumericalScalar prob) const
      {
      if (getDimension() == 1) return computeQuantilePerComponent(prob);
      throw NotYetImplementedException(HERE);
      }

      /*
       * Get the empirical CDF of the sample
       */
00680       NumericalScalar NumericalSampleImplementation::computeEmpiricalCDF(const NumericalPoint & point) const
      {
      const UnsignedLong size = getSize();
      const UnsignedLong dimension = getDimension();
      NumericalScalar cdf = 0.0;
      NumericalScalar p = 1.0 / static_cast<NumericalScalar>(size);
      for (UnsignedLong i = 0; i < size; ++i) {
        NumericalPoint x = (*this)[i];
        UnsignedLong j = 0;
        while ((j < dimension) && (x[j] <= point[j])) ++j;
          if (j == dimension) cdf += p;
      }
      return cdf;
      }

      /* Maximum accessor */
00696       NumericalSampleImplementation::NumericalPoint NumericalSampleImplementation::getMax() const
      {
      const UnsignedLong size = getSize();
      if (size == 0) throw InternalException(HERE) << "Error: impossible to get the maximum of an empty NumericalSample";
      NumericalPoint max = (*this)[0];
      const UnsignedLong dimension = getDimension();
      for(UnsignedLong i = 1; i < size; ++i) {
        NumericalPoint point = (*this)[i];
        for(UnsignedLong j = 0; j < dimension; ++j)
          if (point[j] > max[j]) max[j] = point[j];
      }
      return max;
      }

      /* Minimum accessor */
00711       NumericalSampleImplementation::NumericalPoint NumericalSampleImplementation::getMin() const
      {
      const UnsignedLong size = getSize();
      if (size == 0) throw InternalException(HERE) << "Error: impossible to get the maximum of an empty NumericalSample";
      NumericalPoint min = (*this)[0];
      const UnsignedLong dimension = getDimension();
      for(UnsignedLong i = 1; i < size; ++i) {
        NumericalPoint point = (*this)[i];
        for(UnsignedLong j = 0; j < dimension; ++j)
          if (point[j] < min[j]) min[j] = point[j];
      }
      return min;
      }


      /*
       * Translate realizations in-place
       */
00729       void NumericalSampleImplementation::translate(const NumericalPoint & translation)
      {
      const UnsignedLong size = getSize();
      for (UnsignedLong i = 0; i < size; ++i)
        (*this)[i] += translation;
      }

      /* Get the i-th marginal distribution */
00737       NumericalSampleImplementation NumericalSampleImplementation::getMarginal(const UnsignedLong index) const throw(InvalidArgumentException)
      {
            if (index >= getDimension()) throw InvalidArgumentException(HERE) << "The index of a marginal sample must be in the range [0, dim-1]";
      // Special case for dimension 1
      if (getDimension() == 1) return (*this);
      // General case
      UnsignedLong size(getSize());
      NumericalSampleImplementation marginalSample(size, 1);
      for (UnsignedLong i = 0; i < size; i++)
        {
          // We access directly to the component of the NumericalPoint for performance reason
          marginalSample[i][0] = operator[](i).operator[](index);
        }
      return marginalSample;
      }

      /* Get the distribution of the marginal distribution corresponding to indices dimensions */
00754       NumericalSampleImplementation NumericalSampleImplementation::getMarginal(const Indices & indices) const throw(InvalidArgumentException)
      {
            UnsignedLong dimension(getDimension());
      if (!indices.check(dimension - 1)) throw InvalidArgumentException(HERE) << "The indices of a marginal sample must be in the range [0, dim-1] and  must be different";
      // Special case for dimension 1
      if (dimension == 1) return (*this);
      // General case
            UnsignedLong outputDimension(indices.getSize());
      UnsignedLong size(getSize());
      NumericalSampleImplementation marginalSample(size, outputDimension);
      for (UnsignedLong i = 0; i < size; i++)
        {
          for (UnsignedLong j = 0; j < outputDimension; j++)
            {
            // We access directly to the component of the NumericalPoint for performance reason
            marginalSample[i][j] = operator[](i).operator[](indices[j]);
            }
        }
      return marginalSample;
      }

      /*
       * Scale realizations componentwise in-place
       */
00778       void NumericalSampleImplementation::scale(const NumericalPoint & scaling)
      {
      const UnsignedLong size = getSize();
      const UnsignedLong dimension = getDimension();
      for (UnsignedLong i = 0; i < size; ++i)
        for (UnsignedLong j = 0; j < dimension; ++j)
          (*this)[i][j] *= scaling[j];
      }

      /* Save to CSV file */
00788       void NumericalSampleImplementation::exportToCSVFile(const FileName & filename) const
      {
      String csvSeparator = ResourceMap::GetInstance().get("csv-file-separator");

      std::ofstream csvFile(filename.c_str());
      csvFile.precision(20);

      const UnsignedLong size = getSize();
      const UnsignedLong dimension = getDimension();
      for(UnsignedLong i=0; i < size; ++i, csvFile << '\n') {
        String separator;
        for(UnsignedLong j=0; j < dimension; ++j, separator = csvSeparator)
          csvFile << separator << (*this)[i][j];
      }
      }


      /* Method save() stores the object through the StorageManager */
00806       void NumericalSampleImplementation::save(const StorageManager::Advocate & adv) const
      {
      Type::PersistentCollection<NumericalPoint>::save(adv);
      adv.writeAttribute(StorageManager::DimensionAttribute, dimension_);
      if (!p_description_.isNull())
        adv.writeValue(*p_description_, StorageManager::MemberNameAttribute, "p_description_");
      }

      
      /* Method load() reloads the object from the StorageManager */
00816       void NumericalSampleImplementation::load(const StorageManager::Advocate & adv)
      {
      Type::PersistentCollection<NumericalPoint>::load(adv);
      adv.readAttribute(StorageManager::DimensionAttribute, dimension_);
      Description description;
      adv.readValue(description, StorageManager::MemberNameAttribute, "p_description_");
      if (description.getSize() != 0) setDescription(description);
      }


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

Generated by  Doxygen 1.6.0   Back to index