~ubuntu-branches/ubuntu/maverick/openturns/maverick

//                                               -*- C++ -*-

/**

 *  @file  ImportanceSampling.cxx

 *  @brief ImportanceSampling is an implementation of the importance sampling Monte Carlo simulation method

 *

 *  (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: 2009-05-28 14:47:53 +0200 (jeu. 28 mai 2009) $

 *  Id:      $Id: ImportanceSampling.cxx 1262 2009-05-28 12:47:53Z dutka $

 */

#include "ImportanceSampling.hxx"

#include "ComparisonOperatorImplementation.hxx"

 

namespace OpenTURNS

{

 

  namespace Uncertainty 

  {

 

    namespace Algorithm

    {

 

      /*

       * @class ImportanceSampling

       */

 

      CLASSNAMEINIT(ImportanceSampling);

 

      /* Constructor with parameters */

      ImportanceSampling::ImportanceSampling(const Simulation::Event & event, const Distribution & importanceDistribution) throw(InvalidArgumentException):

        Simulation(event),

        importanceDistribution_(importanceDistribution)

      {

        // Check if the importance distribution dimension is compatible with the event

        if (importanceDistribution.getDimension() != event.getImplementation()->getAntecedent()->getDimension()) throw InvalidArgumentException(HERE) << "The importance distribution must have the same dimension as the event";

      }

 

      /* Virtual constructor */

      ImportanceSampling * ImportanceSampling::clone() const

      {

        return new ImportanceSampling(*this);

      }

 

      /* Compute the block sample */

      ImportanceSampling::NumericalSample ImportanceSampling::computeBlockSample()

      {

        UnsignedLong blockSize(getBlockSize());

        // First, compute a sample of the importance distribution

        NumericalSample inputSample(importanceDistribution_.getNumericalSample(blockSize));

        // Then, evaluate the function on this sample

        NumericalSample blockSample(getEvent().getImplementation()->getFunction()(inputSample));

        // Then, modify in place this sample to take into account the change in the input distribution

        for (UnsignedLong i = 0; i < blockSize; i++)

          {

            inputStrategy_.store(inputSample[i]);

            outputStrategy_.store(blockSample[i]);

            if (getEvent().getOperator()(blockSample[i][0], getEvent().getThreshold()))

              {

                // If the event occured, the value is p_initial(x[i]) / p_importance(x[i])

                // Having access to p_initial is a long trip...

                blockSample[i][0] = getEvent().getImplementation()->getAntecedent()->getDistribution().computePDF(inputSample[i]) / importanceDistribution_.computePDF(inputSample[i]);

              }

            else

              {

                blockSample[i][0] = 0.0;

              }

          }

        return blockSample;

      }

 

      /* Importance distribution accessor */

      ImportanceSampling::Distribution ImportanceSampling::getImportanceDistribution() const

      {

        return importanceDistribution_;

      }

 

      /* String converter */

      String ImportanceSampling::__repr__() const

      {

        OSS oss;

        oss << "class=" << ImportanceSampling::GetClassName()

            << " derived from " << Simulation::__repr__();

        return oss;

      }

 

    } /* namespace Algorithm */

  } /* namespace Uncertainty */

} /* namespace OpenTURNS */