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

StrongMaximumTest.cxx

Go to the documentation of this file.
//                                               -*- C++ -*-
/**
 *  @file  StrongMaximumTest.cxx
 *  @brief StrongMaxTest implements an algorithm to check if a given design point
 *
 *  (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-09-22 11:34:11 +0200 (lun 22 sep 2008) $
 *  Id:      $Id: StrongMaximumTest.cxx 941 2008-09-22 09:34:11Z dutka $
 */
#include <cmath>

#include "StrongMaximumTest.hxx"
#include "Normal.hxx"
#include "IdentityMatrix.hxx"
#include "NearestPointChecker.hxx"
#include "DistributionImplementation.hxx"
#include "ComparisonOperatorImplementation.hxx"

namespace OpenTURNS
{

  namespace Uncertainty
  {

    namespace Algorithm
    {

      /* Initialize static attributes of Analytical::StrongMaximumTest class */
      const NumericalScalar StrongMaximumTest::DefaultDeltaPrecision;
      const NumericalScalar StrongMaximumTest::Epsilon;

      typedef Model::DistributionImplementation                                 DistributionImplementation;
00049       typedef DistributionImplementation::Implementation                        Implementation;
      typedef DistributionImplementation::InverseIsoProbabilisticTransformation InverseIsoProbabilisticTransformation;
      typedef Distribution::Normal                                              Normal;
      typedef Base::Optimisation::NearestPointChecker                           NearestPointChecker;

      CLASSNAMEINIT(StrongMaximumTest);

      /*
       * @class StrongMaximumTest
       * StrongMaximumTest allows to validate a design point
       */

      /* Standard constructor */
      StrongMaximumTest::StrongMaximumTest(const StandardEvent & event,
                                           const NumericalPoint & standardSpaceDesignPoint,
                                 const NumericalScalar importanceLevel,
                                 const NumericalScalar accuracyLevel,
                                 const NumericalScalar confidenceLevel)
      throw(InvalidArgumentException):
      PersistentObject(),
      event_(event),
      standardSpaceDesignPoint_(standardSpaceDesignPoint),
      nearDesignPointVerifyingEventPoints_(NumericalSample(0, standardSpaceDesignPoint)),
      nearDesignPointVerifyingEventValues_(NumericalSample(0, NumericalPoint(1))),
      farDesignPointVerifyingEventPoints_(NumericalSample(0, standardSpaceDesignPoint)),
00074       farDesignPointVerifyingEventValues_(NumericalSample(0, NumericalPoint(1))),
      nearDesignPointViolatingEventPoints_(NumericalSample(0, standardSpaceDesignPoint)),
      nearDesignPointViolatingEventValues_(NumericalSample(0, NumericalPoint(1))),
      farDesignPointViolatingEventPoints_(NumericalSample(0, standardSpaceDesignPoint)),
      farDesignPointViolatingEventValues_(NumericalSample(0, NumericalPoint(1)))
      {
      setImportanceLevel(importanceLevel);
      setAccuracyLevel(accuracyLevel);
      setConfidenceLevel(confidenceLevel);
      initializeParametersGivenConfidenceLevel();
      }

      /* Standard constructor */
      StrongMaximumTest::StrongMaximumTest(const StandardEvent & event,
                                           const NumericalPoint & standardSpaceDesignPoint,
                                 const NumericalScalar importanceLevel,
                                 const NumericalScalar accuracyLevel,
                                 const UnsignedLong pointNumber)
      throw(InvalidArgumentException):
      PersistentObject(),
      event_(event),
      standardSpaceDesignPoint_(standardSpaceDesignPoint),
      nearDesignPointVerifyingEventPoints_(NumericalSample(0, standardSpaceDesignPoint)),
      nearDesignPointVerifyingEventValues_(NumericalSample(0, NumericalPoint(1))),
      farDesignPointVerifyingEventPoints_(NumericalSample(0, standardSpaceDesignPoint)),
      farDesignPointVerifyingEventValues_(NumericalSample(0, NumericalPoint(1))),
00100       nearDesignPointViolatingEventPoints_(NumericalSample(0, standardSpaceDesignPoint)),
      nearDesignPointViolatingEventValues_(NumericalSample(0, NumericalPoint(1))),
      farDesignPointViolatingEventPoints_(NumericalSample(0, standardSpaceDesignPoint)),
      farDesignPointViolatingEventValues_(NumericalSample(0, NumericalPoint(1)))
      {
      if (standardSpaceDesignPoint.norm() < Epsilon) throw InvalidArgumentException(HERE) << "Error: the given standard space design point is too close to the origin for the strong maximum test to work, norm=" << standardSpaceDesignPoint.norm();
00106       setImportanceLevel(importanceLevel);
      setAccuracyLevel(accuracyLevel);
      setPointNumber(pointNumber);
      initializeParametersGivenPointNumber();
      }

00112       /* Virtual constructor */
      StrongMaximumTest * StrongMaximumTest::clone() const
      {
      return new StrongMaximumTest(*this);
      }

      /* Destructor */
      StrongMaximumTest::~StrongMaximumTest()
      {
      // Nothing to do
00122       }

      /* standardSpaceDesignPoint accessor */
      void StrongMaximumTest::setStandardSpaceDesignPoint(const NumericalPoint & standardSpaceDesignPoint)
      {
      if(standardSpaceDesignPoint.norm() <= 0.0) 
00128         {
          throw InvalidRangeException(HERE) << "DesignPoint with norm <= 0.0";
        }
      standardSpaceDesignPoint_ = standardSpaceDesignPoint;
      }

00134       /* standardSpaceDesignPoint accessor */
      StrongMaximumTest::NumericalPoint StrongMaximumTest::getStandardSpaceDesignPoint() const 
      {
      return standardSpaceDesignPoint_;
      }

00140       /* deltaEpsilon accessor */
      NumericalScalar StrongMaximumTest::getDeltaEpsilon() const
      {
      return deltaEpsilon_;
      }

00146       /* Event accessor */
      StrongMaximumTest::StandardEvent StrongMaximumTest::getEvent() const
      {
      return event_;
      }

      /* ImportanceLevel  accessor */
      NumericalScalar StrongMaximumTest::getImportanceLevel() const
      {
      return importanceLevel_;
      }
00157 
      /* ImportanceLevel accessor */
      void StrongMaximumTest::setImportanceLevel(const NumericalScalar importanceLevel)
      throw(InvalidRangeException)
      {
      if((importanceLevel >= 1.0) or (importanceLevel <= 0.0))
00163         {
          throw InvalidRangeException(HERE) <<"importanceLevel is not within 0.0 and 1.0";
        }
      importanceLevel_ = importanceLevel ;
      }

      /* AccuracyLevel  accessor */
      NumericalScalar StrongMaximumTest::getAccuracyLevel() const
      {
      return accuracyLevel_;
      }
00174 
      /* AccuracyLevel accessor */
      void StrongMaximumTest::setAccuracyLevel(const NumericalScalar accuracyLevel)
      throw(InvalidRangeException)
      {
      if (accuracyLevel <= 1.0) 
00180         {
          throw InvalidRangeException(HERE) <<"accuracyLevel is not > 1.0";
        }
      accuracyLevel_ = accuracyLevel ;
      }

      /* ConfidenceLevel accessor */
      NumericalScalar StrongMaximumTest::getConfidenceLevel() const 
      {
      return confidenceLevel_;
      }
00191 
      /* ConfidenceLevel accessor */
      void StrongMaximumTest::setConfidenceLevel(const NumericalScalar confidenceLevel)
      throw(InvalidRangeException)
      {
      if((confidenceLevel > 1.0) or (confidenceLevel < 0.0))
00197         {
          throw InvalidRangeException(HERE) << "confidenceLevel is not within 0.0 and 1.0";
        }
      confidenceLevel_ = confidenceLevel ;
      }

      /* DesignPointVicinity accessor */
      NumericalScalar StrongMaximumTest::getDesignPointVicinity() const
      {
      return designPointVicinity_;
      }
00208 
      /* StrongMaxTestDesignPointVicinity accessor */
      void  StrongMaximumTest::setDesignPointVicinity(const NumericalScalar designPointVicinity)
      throw(InvalidRangeException)
      {
      if((designPointVicinity > 1.0) or (designPointVicinity < 0.0))
00214         {
          throw InvalidRangeException(HERE) <<"designPointVicinity is not within 0.0 and 1.0";
        }
      designPointVicinity_ =  designPointVicinity;
      }

      /* StrongMaxTestPointNumber  accessor */
      UnsignedLong  StrongMaximumTest::getPointNumber() const
      {
      return pointNumber_;
      }
00225 
      /* StrongMaxTestPointNumber accessor */
      void StrongMaximumTest::setPointNumber(const UnsignedLong pointNumber)
      throw(InvalidRangeException)
      {
      if(pointNumber == 0)
00231         {
          throw InvalidRangeException(HERE) << "pointNumber is equal to 0";
        }
      pointNumber_ = pointNumber ;
      }

00237       /* NearDesignPointVerifyingEventPoints accessor */
      StrongMaximumTest::NumericalSample StrongMaximumTest::getNearDesignPointVerifyingEventPoints() const
      {
      return nearDesignPointVerifyingEventPoints_;
      }

00243       /* NearDesignPointVerifyingEventValues accessor */
      StrongMaximumTest::NumericalSample StrongMaximumTest::getNearDesignPointVerifyingEventValues() const
      {
      return nearDesignPointVerifyingEventValues_;
      }

00249       /* FarDesignPointVerifyingEventPoints accessor */
      StrongMaximumTest::NumericalSample StrongMaximumTest::getFarDesignPointVerifyingEventPoints() const
      {
      return farDesignPointVerifyingEventPoints_;
      }

00255       /* FarDesignPointVerifyingEventValues accessor */
      StrongMaximumTest::NumericalSample StrongMaximumTest::getFarDesignPointVerifyingEventValues() const
      {
      return farDesignPointVerifyingEventValues_;
      }

00261       /* NearDesignPointViolatingEventPoints accessor */
      StrongMaximumTest::NumericalSample StrongMaximumTest::getNearDesignPointViolatingEventPoints() const
      {
      return nearDesignPointViolatingEventPoints_;
      }

00267       /* NearDesignPointViolatingEventValues accessor */
      StrongMaximumTest::NumericalSample StrongMaximumTest::getNearDesignPointViolatingEventValues() const
      {
      return nearDesignPointViolatingEventValues_;
      }

      /* FarDesignPointViolatingEventPoints accessor */
00274       StrongMaximumTest::NumericalSample StrongMaximumTest::getFarDesignPointViolatingEventPoints() const
      {
      return farDesignPointViolatingEventPoints_;
      }

      /* FarDesignPointViolatingEventValues accessor */
      StrongMaximumTest::NumericalSample StrongMaximumTest::getFarDesignPointViolatingEventValues() const
      {
      return farDesignPointViolatingEventValues_;
      }


      /* String converter */
      String StrongMaximumTest::str() const
      {
      OSS oss;
      oss << "class=" << StrongMaximumTest::GetClassName()
          << " event=" << event_
          << " standardSpaceDesignPoint=" << standardSpaceDesignPoint_
          << " importanceLevel=" << importanceLevel_
          << " accuracyLevel=" << accuracyLevel_
          << " confidenceLevel=" << confidenceLevel_
          << " designPointVicinity=" << designPointVicinity_
          << " pointNumber=" << pointNumber_
00298           << " deltaEpsilon=" << deltaEpsilon_
          << " nearDesignPointVerifyingEventPoints=" << nearDesignPointVerifyingEventPoints_
          << " nearDesignPointVerifyingEventValues=" << nearDesignPointVerifyingEventValues_
          << " farDesignPointVerifyingEventPoints=" << farDesignPointVerifyingEventPoints_
          << " farDesignPointVerifyingEventValues=" << farDesignPointVerifyingEventValues_
          << " nearDesignPointViolatingEventPoints=" << nearDesignPointViolatingEventPoints_
          << " nearDesignPointViolatingEventValues=" << nearDesignPointViolatingEventValues_
          << " farDesignPointViolatingEventPoints=" << farDesignPointViolatingEventPoints_
          << " farDesignPointViolatingEventValues=" << farDesignPointViolatingEventValues_;
      return oss;
      }

      /* Initialize Strong Max Test Parameters : method 1 */
      void StrongMaximumTest::initializeParametersGivenConfidenceLevel()
      throw(InvalidRangeException)
      {
      /* evaluate the intermediate parameter delta_epsilon (see documentation) */
        deltaEpsilon_ = computeDeltaEpsilon();

      /* evaluate the HyperSphereSurfaceRatio (see documentation) */
      NumericalScalar p(computeHyperSphereSurfaceRatio());
      // put eps1 and eps2 instead of 1.0 and 0.0
      if((p >= 1.0) or (p <= 0.0))
00321         {
          throw InvalidRangeException(HERE) << "hyperSphereSurfaceRatio is not strictly within 0.0 and 1.0";
        }

      /* evaluate and affect the pointNumber_ */
      setPointNumber(UnsignedLong(round(log(1.0 - confidenceLevel_) / log(1.0 - p))));


      /* evaluate and affect the designPointVicinity_ */
      setDesignPointVicinity(1.0 / (1.0 + accuracyLevel_ * deltaEpsilon_));
      }

      /* Initialize Strong Max Test Parameters : method 2 */
      void StrongMaximumTest::initializeParametersGivenPointNumber()
      throw(InvalidRangeException)
      {
      /* evaluate the intermediate parameter delta_epsilon (see documentation) */
      deltaEpsilon_ = computeDeltaEpsilon();

      /* evaluate the HyperSphereSurfaceRatio (see documentation) */
      NumericalScalar p(computeHyperSphereSurfaceRatio());

      // put eps1 and eps2 instead of 1.0 and 0.0
00344       if((p >= 1.0) or (p <= 0.0))
        {
          throw InvalidRangeException(HERE) << "hyperSphereSurfaceRatio is not strictly within 0.0 and 1.0";
        }
      /* evaluate and affect the confidenceLevel */

      setConfidenceLevel(1.0 - pow(1.0 - p, pointNumber_));

      /* evaluate and affect the designPointVicinity_ */
      setDesignPointVicinity(1.0 / (1.0 + accuracyLevel_ * deltaEpsilon_));
      }

      /*  the function that evaluates the HyperSphereSurfaceRatio (see documentation) */
      NumericalScalar StrongMaximumTest::computeHyperSphereSurfaceRatio()
      {
      UnsignedLong dimension(standardSpaceDesignPoint_.getDimension());
      NumericalScalar a( acos((1.0 + deltaEpsilon_) / ( 1.0 + accuracyLevel_ * deltaEpsilon_) ) );
      NumericalScalar sinA(sin(a));
      NumericalScalar squareSinA(sinA * sinA);
      NumericalScalar sum(0.0);
      /* even dimension  */
      if (dimension % 2 == 0)  
        {
          UnsignedLong indexMax(dimension / 2 - 1);
          NumericalScalar u(sinA);
          for (UnsignedLong index = 0; index< indexMax; index++)
            {
            sum += u;
            u *= (1.0 - 1.0 / (2.0 * index + 3.0)) * squareSinA;
            }
          /* M_1_PI = 1/PI cf cmath */
          return M_1_PI * (a - cos(a) * sum);
        }
      else
        /* odd dimension  */  
00379         {
          UnsignedLong indexMax((dimension - 1) / 2);
          NumericalScalar u(1.0);
          for (UnsignedLong index = 0; index< indexMax; index++)
            {
            sum += u;
            u *= (1.0 - 1.0 / (2.0 * index + 2.0)) * squareSinA;
            }
          return 0.5 * (1.0 - cos(a) * sum);
        }
      }

      /*  the function that evaluates the  intermediate parameter delta_epsilon (see documentation) */
      NumericalScalar StrongMaximumTest::computeDeltaEpsilon()
      {
      /* evaluate the reliability index */
      NumericalScalar betaSquare(standardSpaceDesignPoint_.norm2());

      /* get the input distribution in the standard space */
        Implementation p_inputStandardDistribution(event_.getImplementation()->getAntecedent()->getDistribution().getImplementation());

      /* evaluate the generator at beta square */
      NumericalScalar pdfMin(importanceLevel_ * p_inputStandardDistribution->computeDensityGenerator(betaSquare));

      /* research the interval [deltaMin deltaMax] including the solution */
      NumericalScalar deltaMax(1.0);

      while ( p_inputStandardDistribution->computeDensityGenerator(betaSquare * pow(1.0 + deltaMax, 2)) > pdfMin )
        {
          deltaMax++;
        }
      NumericalScalar deltaMin(deltaMax - 1.0);

      /* we proceed to the dichotomie on [deltaMin deltaMax] */
      NumericalScalar deltaMiddle(0.0);
      while ( (deltaMax - deltaMin) > DefaultDeltaPrecision )
        {
          /* we evaluate the middle of  [deltaMin deltaMax] */
          deltaMiddle = 0.5 * (deltaMax + deltaMin);
00418           if(  p_inputStandardDistribution->computeDensityGenerator(betaSquare * pow(1.0 + deltaMiddle, 2)) > pdfMin )
            {
            deltaMin = deltaMiddle;
            }
          else 
            {
00424             deltaMax = deltaMiddle;
            }
        }
      return 0.5 * (deltaMax + deltaMin);
      }

      /* the function that evaluate if a point is in the vicinity of the design point */
      Bool StrongMaximumTest::isInTheVicinityOfTheDesignPoint(const NumericalPoint & numericalPoint)
      { 
      return (NumericalPoint::dot(numericalPoint, standardSpaceDesignPoint_) > numericalPoint.norm() * standardSpaceDesignPoint_.norm() * designPointVicinity_);
      }

      /* The function that runs the Strong Max Test */
      void StrongMaximumTest::run()
      throw(InvalidArgumentException, InternalException)
      {
      /* prepare test parameters */

      /* radius of the inner sphere */
      NumericalScalar beta(standardSpaceDesignPoint_.norm());
      /* radius of the sphere to be sampled */
      NumericalScalar radius = beta * (1.0 + accuracyLevel_ * deltaEpsilon_);
      /* sample of the sphere */
      NumericalSample sample(sampleSphere(radius, standardSpaceDesignPoint_.getDimension(), pointNumber_));
      /* create a nearestPointChecker, in charge of the evaluation of the level function over the sample and to classify the points according to the operator and the threshold */
      NearestPointChecker nearestPointChecker(event_.getImplementation()->getFunction(), event_.getOperator(), event_.getThreshold(), sample);
      /* access to the inverse isoprobabilistic transformation */
      InverseIsoProbabilisticTransformation inverseIsoProbabilisticTransformation(event_.getImplementation()->getAntecedent()->getDistribution().getInverseIsoProbabilisticTransformation());
      /* run test */
      try{
        nearestPointChecker.run();
      }
      catch(InvalidArgumentException & ex)
        {
          throw InvalidArgumentException(HERE) << ex.str();
        }
      catch(InternalException & ex)
        {
          throw  InternalException(HERE) << ex.str();
        }

      /* get nearestPointChecker result */
      NearestPointChecker::Result nearestPointCheckerResult(nearestPointChecker.getResult());
      /* split the two samples according to the vicinity of the design point
       * everything is done in place, using the attributs of the class in order
       * to limit the memory usage */

      nearDesignPointVerifyingEventPoints_ = nearestPointCheckerResult.getVerifyingConstraintPoints();
      nearDesignPointVerifyingEventValues_ = nearestPointCheckerResult.getVerifyingConstraintValues();

      UnsignedLong sampleSize(nearDesignPointVerifyingEventPoints_.getSize());
      /* If there is something to classify */
      if (sampleSize > 0)
        {
          UnsignedLong leftCounter(0);
          UnsignedLong rightCounter(sampleSize - 1);

          /* we sort among the nearDesignPointVerifyingEventPoints_ (ie which realise the event) the ones which are in the vicinity of the design point */
          while (leftCounter < rightCounter)
            {
            if (isInTheVicinityOfTheDesignPoint(nearDesignPointVerifyingEventPoints_[leftCounter]))
              {
                /* we leave at the beginning of the sample all the points (and the corresponding values) in the vicinity of the design point */
                leftCounter ++;
              }
            else
              {
                /* we put at the end of the sample  all the points (and the corresponding values) not in the vicinity of the design point */
                NumericalPoint point(nearDesignPointVerifyingEventPoints_[leftCounter]);
                NumericalPoint value(nearDesignPointVerifyingEventValues_[leftCounter]);
                nearDesignPointVerifyingEventPoints_[leftCounter] = nearDesignPointVerifyingEventPoints_[rightCounter];
                nearDesignPointVerifyingEventValues_[leftCounter] = nearDesignPointVerifyingEventValues_[rightCounter];
                nearDesignPointVerifyingEventPoints_[rightCounter] = point;
                nearDesignPointVerifyingEventValues_[rightCounter] = value;
                rightCounter --;
              }
            }
          /* At the end, we still have to check the point at the position leftCounter but without updating rightCounter, which might be already equals to 0 and then might try to become < 0 */
          if (isInTheVicinityOfTheDesignPoint(nearDesignPointVerifyingEventPoints_[leftCounter])) leftCounter++;
          /* substitute physical points to standard points */
          nearDesignPointVerifyingEventPoints_ = inverseIsoProbabilisticTransformation(nearDesignPointVerifyingEventPoints_);

          /* we build the final two NumericalSamples (points, values) for each group */
          /* we split the sortedVerifyingConstraintPoints and the sortedVerifyingConstraintValues in 2 NumericalSamples each : one with the left side (points in the vicinity of the design point) and the other with the right side (points not in the vicinity of the design point) */
          if (leftCounter < sampleSize)
            {
            farDesignPointVerifyingEventPoints_ = nearDesignPointVerifyingEventPoints_.split(leftCounter);
            farDesignPointVerifyingEventValues_ = nearDesignPointVerifyingEventValues_.split(leftCounter);
            }
        }
      /* we do the same thing for points which violate the constraints (ie which don't realise the event) */
      farDesignPointViolatingEventPoints_ = nearestPointCheckerResult.getViolatingConstraintPoints();
      farDesignPointViolatingEventValues_ = nearestPointCheckerResult.getViolatingConstraintValues();

      sampleSize = farDesignPointViolatingEventPoints_.getSize();
      /* If there is something to classify */
      if (sampleSize > 0)
        {
          UnsignedLong leftCounter(0);
          UnsignedLong rightCounter(sampleSize - 1);

          /* we sort among the nearDesignPointViolatingEventPoints_ (ie which realise the event) the ones which are in the vicinity of the design point */
          while (leftCounter < rightCounter)
            {
            if (isInTheVicinityOfTheDesignPoint(farDesignPointViolatingEventPoints_[leftCounter]))
              {
                /* we put at the end of the sample  all the points (and the corresponding values) not in the vicinity of the design point */
                NumericalPoint point(farDesignPointViolatingEventPoints_[leftCounter]);
                NumericalPoint value(farDesignPointViolatingEventValues_[leftCounter]);
                farDesignPointViolatingEventPoints_[leftCounter] = farDesignPointViolatingEventPoints_[rightCounter];
                farDesignPointViolatingEventValues_[leftCounter] = farDesignPointViolatingEventValues_[rightCounter];
                farDesignPointViolatingEventPoints_[rightCounter] = point;
                farDesignPointViolatingEventValues_[rightCounter] = value;
                rightCounter --;
              }
            else
              {
                /* we leave at the beginning of the sample all the points (and the corresponding values) in the vicinity of the design point */
                leftCounter ++;
              }
            }
          /* At the end, we still have to check the point at the position leftCounter but without updating rightCounter, which should be already equals to 0 and then could try to become < 0 */
          if (!isInTheVicinityOfTheDesignPoint(farDesignPointViolatingEventPoints_[leftCounter])) leftCounter++;
          /* substitute physical points to standard points */
00548           farDesignPointViolatingEventPoints_ = inverseIsoProbabilisticTransformation(farDesignPointViolatingEventPoints_);

          /* we build the final two NumericalSamples (points, values) for each group */
          /* we split the sortedViolatingConstraintPoints and the sortedViolatingConstraintValues in 2 NumericalSamples each : one with the left side (points in the vicinity of the design point) and the other with the right side (points not in the vicinity of the design point) */
          if (leftCounter < sampleSize)
            {
            nearDesignPointViolatingEventPoints_ = farDesignPointViolatingEventPoints_.split(leftCounter);
            nearDesignPointViolatingEventValues_ = farDesignPointViolatingEventValues_.split(leftCounter);
            }
        }
      }

      /* the function that samples the sphere (radius) with N points */
      StrongMaximumTest::NumericalSample StrongMaximumTest::sampleSphere(const NumericalScalar radius,
                                                                         const UnsignedLong dimension,
                                                                         const UnsignedLong pointNumber) const
      {
      // First, generate a sample of a standard normal distribution of the proper size and dimension
      Normal standardNormal(dimension);
      NumericalSample sample(standardNormal.getNumericalSample(pointNumber));
      // Then, normalize the points to have length radius
      for (UnsignedLong i = 0; i < pointNumber; i++)
        {
          NumericalScalar norm(sample[i].norm());
          // If the point is the origin, we reject it
          while (norm == 0.0)
            {
            sample[i] = standardNormal.getRealization();
            norm = sample[i].norm();
            }
          sample[i] *= (radius / norm);
        }
      sample.setName("StrongMaximumTest sample");
      // The normalize sample follow the uniform distribution over the hypersphere
      return sample;
      }

    } /* namespace Algorithm */
  } /* namespace Uncertainty */
} /* namespace OpenTURNS */

Generated by  Doxygen 1.6.0   Back to index