StatisticalInverseProblem.h

Go to the documentation of this file.

//-----------------------------------------------------------------------bl-
//--------------------------------------------------------------------------
//
// QUESO - a library to support the Quantification of Uncertainty
// for Estimation, Simulation and Optimization
//
// Copyright (C) 2008-2017 The PECOS Development Team
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the Version 2.1 GNU Lesser General
// Public License as published by the Free Software Foundation.
//
// 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. 51 Franklin Street, Fifth Floor,
// Boston, MA  02110-1301  USA
//
//-----------------------------------------------------------------------el-

#ifndef UQ_SIP_H
#define UQ_SIP_H

#include <queso/StatisticalInverseProblemOptions.h>
#include <queso/MetropolisHastingsSG.h>
#include <queso/MLSampling.h>
#include <queso/InstantiateIntersection.h>
#include <queso/VectorRealizer.h>
#include <queso/SequentialVectorRealizer.h>
#include <queso/VectorRV.h>
#include <queso/ScalarFunction.h>
#include <queso/GPMSA.h>
#include <queso/ScopedPtr.h>

namespace QUESO {

class GslVector;
class GslMatrix;

/* OLD STUFF: The constructor of the 'scalar function' asks for the user to specify which
 * value the routine is actually computing, so that the 'scalar function' class can properly
 * implements both class operations 'actualValue()' and 'minus2LnValue()'*/

template <class P_V = GslVector, class P_M = GslMatrix>
class StatisticalInverseProblem
{
public:


  StatisticalInverseProblem(const char*                               prefix,
                                   const SipOptionsValues*            alternativeOptionsValues, // dakota
                                   const BaseVectorRV      <P_V,P_M>& priorRv,
                                   const BaseScalarFunction<P_V,P_M>& likelihoodFunction,
                                         GenericVectorRV   <P_V,P_M>& postRv);


  StatisticalInverseProblem(const char * prefix,
                            const SipOptionsValues * alternativeOptionsValues,
                            const GPMSAFactory<P_V,P_M> & gpmsaFactory,
                            GenericVectorRV<P_V,P_M> & postRv);

  ~StatisticalInverseProblem();


  bool                             computeSolutionFlag             () const;


  void solveWithBayesMetropolisHastings(const MhOptionsValues* alternativeOptionsValues, // dakota
                                        const P_V&                    initialValues,
                                        const P_M*                    initialProposalCovMatrix);


  void solveWithBayesMetropolisHastings(const MhOptionsValues * alternativeOptionsValues,
                                        const P_V & initialValues);


  void solveWithBayesMetropolisHastings(const MhOptionsValues * alternativeOptionsValues);


  void solveWithBayesMetropolisHastings();


  void seedWithMAPEstimator();

  void                             solveWithBayesMLSampling        ();

  const MetropolisHastingsSG<P_V, P_M> & sequenceGenerator() const;

  const BaseVectorRV   <P_V,P_M>& priorRv                   () const;


  const GenericVectorRV<P_V,P_M>& postRv                    () const;


  const BaseVectorSequence<P_V,P_M>& chain() const;


  const ScalarSequence<double>& logLikelihoodValues() const;


  const ScalarSequence<double>& logTargetValues() const;

  double                           logEvidence                     () const;

  double                           meanLogLikelihood               () const;

  //\TODO Related to ML.
  double                           eig                             () const;



  void print(std::ostream& os) const;
  friend std::ostream& operator<<(std::ostream& os,
      const StatisticalInverseProblem<P_V,P_M>& obj) {
    obj.print(os);
    return os;
  }
private:
  const BaseEnvironment&                 m_env;

  const BaseVectorRV        <P_V,P_M>&   m_priorRv;
  const BaseScalarFunction  <P_V,P_M>&   m_likelihoodFunction;
        GenericVectorRV     <P_V,P_M>&   m_postRv;

  typename ScopedPtr<VectorSet           <P_V,P_M> >::Type m_solutionDomain;
  typename ScopedPtr<BaseJointPdf        <P_V,P_M> >::Type m_solutionPdf;
  typename ScopedPtr<BaseVectorMdf       <P_V,P_M> >::Type m_subSolutionMdf;
  typename ScopedPtr<BaseVectorCdf       <P_V,P_M> >::Type m_subSolutionCdf;
  typename ScopedPtr<BaseVectorRealizer  <P_V,P_M> >::Type m_solutionRealizer;

  typename ScopedPtr<MetropolisHastingsSG<P_V,P_M> >::Type m_mhSeqGenerator;
  typename ScopedPtr<MLSampling          <P_V,P_M> >::Type m_mlSampler;
  typename ScopedPtr<BaseVectorSequence  <P_V,P_M> >::Type m_chain;
  ScopedPtr<ScalarSequence<double> >::Type m_logLikelihoodValues;
  ScopedPtr<ScalarSequence<double> >::Type m_logTargetValues;

  ScopedPtr<const SipOptionsValues>::Type m_optionsObj;

  bool m_seedWithMAPEstimator;

#ifdef UQ_ALSO_COMPUTE_MDFS_WITHOUT_KDE
  typename ScopedPtr<ArrayOfOneDGrids    <P_V,P_M> > m_subMdfGrids;
  typename ScopedPtr<ArrayOfOneDTables   <P_V,P_M> > m_subMdfValues;
#endif
};

}  // End namespace QUESO

#endif // UQ_SIP_H