GPMSA.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-2015 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_GPMSA_HELPER_H
#define UQ_GPMSA_HELPER_H

#include <vector>

#include <queso/Environment.h>
#include <queso/ScalarFunction.h>
#include <queso/VectorSpace.h>
#include <queso/VectorRV.h>
#include <queso/ConcatenatedVectorRV.h>
#include <queso/GammaVectorRV.h>
#include <queso/BetaVectorRV.h>
#include <queso/UniformVectorRV.h>
#include <queso/GPMSAOptions.h>
#include <queso/ScopedPtr.h>

namespace QUESO {

class GslVector;
class GslMatrix;

template <class V = GslVector, class M = GslMatrix>
class GPMSAEmulator : public BaseScalarFunction<V, M>
{
public:
  GPMSAEmulator(const VectorSet<V, M> & domain,
                const VectorSpace<V, M> & m_scenarioSpace,
                const VectorSpace<V, M> & m_parameterSpace,
                const VectorSpace<V, M> & m_simulationOutputSpace,
                const VectorSpace<V, M> & m_experimentOutputSpace,
                const unsigned int m_numSimulations,
                const unsigned int m_numExperiments,
                const std::vector<V *> & m_simulationScenarios,
                const std::vector<V *> & m_simulationParameters,
                const std::vector<V *> & m_simulationOutputs,
                const std::vector<V *> & m_experimentScenarios,
                const std::vector<V *> & m_experimentOutputs,
                const std::vector<V>   & m_discrepancyBases,
                const std::vector<M>   & m_observationErrorMatrices,
                const M & m_experimentErrors,
                const ConcatenatedVectorRV<V, M> & m_totalPrior,
                const V & residual_in,
                const M & BT_Wy_B_inv_in,
                const M & KT_K_inv_in);

  virtual ~GPMSAEmulator();

  virtual double lnValue(const V & domainVector,
                         const V * domainDirection,
                         V * gradVector,
                         M * hessianMatrix,
                         V * hessianEffect) const;

  virtual double actualValue(const V & domainVector,
                             const V * domainDirection,
                             V * gradVector,
                             M * hessianMatrix,
                             V * hessianEffect) const;

  const VectorSpace<V, M> & m_scenarioSpace;
  const VectorSpace<V, M> & m_parameterSpace;
  const VectorSpace<V, M> & m_simulationOutputSpace;
  const VectorSpace<V, M> & m_experimentOutputSpace;

  const unsigned int m_numSimulations;
  const unsigned int m_numExperiments;

  const std::vector<V *> & m_simulationScenarios;
  const std::vector<V *> & m_simulationParameters;
  const std::vector<V *> & m_simulationOutputs;
  const std::vector<V *> & m_experimentScenarios;
  const std::vector<V *> & m_experimentOutputs;

        std::vector<V>     m_discrepancyBases;

  const std::vector<M> & m_observationErrorMatrices;

  //
  // Intermediate calculations we can cache
  //
  unsigned int num_svd_terms;

  // Total observation error covriance matrix
  const M & m_experimentErrors;

  const ConcatenatedVectorRV<V, M> & m_totalPrior;

  //
  // Intermediate calculations cached by factory
  //
  const V & residual;

  const M & BT_Wy_B_inv;

  const M & KT_K_inv;
};

template <class V = GslVector, class M = GslMatrix>
class GPMSAFactory
{
public:
  GPMSAFactory(const BaseEnvironment & env,
               GPMSAOptions * opts,
               const BaseVectorRV<V, M> & parameterPrior,
               const VectorSpace<V, M> & scenarioSpace,
               const VectorSpace<V, M> & parameterSpace,
               const VectorSpace<V, M> & simulationOutputSpace,
               const VectorSpace<V, M> & experimentOutputSpace,
               unsigned int numSimulations,
               unsigned int numExperiments);

  ~GPMSAFactory();



  unsigned int numSimulations() const;

  unsigned int numExperiments() const;

  const VectorSpace<V, M> & scenarioSpace() const;

  const VectorSpace<V, M> & parameterSpace() const;

  const VectorSpace<V, M> & simulationOutputSpace() const;

  const VectorSpace<V, M> & experimentOutputSpace() const;


  const V & simulationScenario(unsigned int simulationId) const;


  const std::vector<V *> & simulationScenarios() const;


  const V & simulationParameter(unsigned int simulationId) const;


  const std::vector<V *> & simulationParameters() const;


  const V & simulationOutput(unsigned int simulationId) const;


  const std::vector<V *> & simulationOutputs() const;


  const V & experimentScenario(unsigned int experimentId) const;


  const std::vector<V *> & experimentScenarios() const;


  const V & experimentOutput(unsigned int experimentId) const;


  const std::vector<V *> & experimentOutputs() const;

  const M & experimentErrors() const;

  const BaseEnvironment & env() const;

  const GPMSAEmulator<V, M> & getGPMSAEmulator() const;



  void addSimulation(V & simulationScenario,
                     V & simulationParameter,
                     V & simulationOutput);


  void addSimulations(const std::vector<V *> & simulationScenarios,
                      const std::vector<V *> & simulationParameters,
                      const std::vector<V *> & simulationOutputs);


  void addExperiments(const std::vector<V *> & experimentScenarios,
                      const std::vector<V *> & experimentOutputs,
                      const M * experimentErrors);


  void setDiscrepancyBases(const std::vector<V *> & discrepancyBases);

  M & getObservationErrorCovariance(unsigned int simulationNumber);

  const M & getObservationErrorCovariance(unsigned int simulationNumber) const;

  const ConcatenatedVectorRV<V, M> & prior() const;

  void print(std::ostream& os) const;
  friend std::ostream & operator<<(std::ostream& os,
                                   const GPMSAFactory<V, M> & obj)
  {
    obj.print(os);
    return os;
  }

  const BaseEnvironment & m_env;

  const BaseVectorRV<V, M> & m_parameterPrior;

  const VectorSpace<V, M> & m_scenarioSpace;
  const VectorSpace<V, M> & m_parameterSpace;
  const VectorSpace<V, M> & m_simulationOutputSpace;
  const VectorSpace<V, M> & m_experimentOutputSpace;

  unsigned int m_numSimulations;
  unsigned int m_numExperiments;

  std::vector<V *> m_simulationScenarios;
  std::vector<V *> m_simulationParameters;
  std::vector<V *> m_simulationOutputs;
  std::vector<V *> m_experimentScenarios;
  std::vector<V *> m_experimentOutputs;

  // We will be recentering data around the simulation output mean
  typename ScopedPtr<V>::Type simulationOutputMeans;

  std::vector<V> m_discrepancyBases;

  std::vector<M> m_observationErrorMatrices;

  // Total observation error covriance matrix
  const M * m_experimentErrors;

  // Counter for the number of adds that happen
  unsigned int m_numSimulationAdds;
  unsigned int m_numExperimentAdds;

  // The space in which the emulator (simulator) lives
  // const VectorSpace<V, M> & m_emulatorSpace;

  // The emulator state
  // const V & m_emulator;

  // Build the emulator once all data has been added
  void setUpEmulator();

  // All the GP priors information for a scalar GP follows:
  void setUpHyperpriors();

  // Domains for all the hyperpriors
  typename ScopedPtr<VectorSpace<V, M> >::Type oneDSpace;

  // Emulator mean
  typename ScopedPtr<V>::Type emulatorMeanMin;
  typename ScopedPtr<V>::Type emulatorMeanMax;
  typename ScopedPtr<BoxSubset<V, M> >::Type emulatorMeanDomain;

  // Emulator precision
  typename ScopedPtr<VectorSpace<V, M> >::Type emulatorPrecisionSpace;
  typename ScopedPtr<V>::Type emulatorPrecisionMin;
  typename ScopedPtr<V>::Type emulatorPrecisionMax;
  typename ScopedPtr<BoxSubset<V, M> >::Type emulatorPrecisionDomain;

  // Emulator correlation strength
  typename ScopedPtr<VectorSpace<V, M> >::Type emulatorCorrelationSpace;
  typename ScopedPtr<V>::Type emulatorCorrelationMin;
  typename ScopedPtr<V>::Type emulatorCorrelationMax;
  typename ScopedPtr<BoxSubset<V, M> >::Type emulatorCorrelationDomain;

  // Observational precision
  typename ScopedPtr<VectorSpace<V, M> >::Type observationalPrecisionSpace;
  typename ScopedPtr<V>::Type observationalPrecisionMin;
  typename ScopedPtr<V>::Type observationalPrecisionMax;
  typename ScopedPtr<BoxSubset<V, M> >::Type observationalPrecisionDomain;

  // Discrepancy precision
  typename ScopedPtr<V>::Type discrepancyPrecisionMin;
  typename ScopedPtr<V>::Type discrepancyPrecisionMax;
  typename ScopedPtr<BoxSubset<V, M> >::Type discrepancyPrecisionDomain;

  // Discrepancy correlation strength
  typename ScopedPtr<VectorSpace<V, M> >::Type discrepancyCorrelationSpace;
  typename ScopedPtr<V>::Type discrepancyCorrelationMin;
  typename ScopedPtr<V>::Type discrepancyCorrelationMax;
  typename ScopedPtr<BoxSubset<V, M> >::Type discrepancyCorrelationDomain;

  // Emulator data precision
  typename ScopedPtr<V>::Type emulatorDataPrecisionMin;
  typename ScopedPtr<V>::Type emulatorDataPrecisionMax;
  typename ScopedPtr<BoxSubset<V, M> >::Type emulatorDataPrecisionDomain;

  // Now form full prior
  typename ScopedPtr<VectorSpace<V, M> >::Type totalSpace;
  typename ScopedPtr<V>::Type totalMins;
  typename ScopedPtr<V>::Type totalMaxs;

  typename ScopedPtr<BoxSubset<V, M> >::Type totalDomain;

  std::vector<const BaseVectorRV<V, M> *> priors;

  // The hyperpriors
  typename ScopedPtr<UniformVectorRV<V, M> >::Type m_emulatorMean;  // scalar
  typename ScopedPtr<GammaVectorRV<V, M> >::Type m_emulatorPrecision;  // (dim num_svd_terms) gamma(a, b) shape-rate
  typename ScopedPtr<GammaVectorRV<V, M> >::Type m_observationalPrecision;  // scalar gamma(a, b) shape-rate
  typename ScopedPtr<BetaVectorRV<V, M> >::Type m_emulatorCorrelationStrength;  // (dim scenariosspace + dim parameterspace)
  typename ScopedPtr<GammaVectorRV<V, M> >::Type m_discrepancyPrecision;  // (scalar) shape-rate
  typename ScopedPtr<BetaVectorRV<V, M> >::Type m_discrepancyCorrelationStrength;  // (dim scenariospace)
  typename ScopedPtr<GammaVectorRV<V, M> >::Type m_emulatorDataPrecision;  // (scalar) shape-rate
  typename ScopedPtr<ConcatenatedVectorRV<V, M> >::Type m_totalPrior;  // prior for joint parameters and hyperparameters

  typename ScopedPtr<V>::Type m_emulatorPrecisionShapeVec;
  typename ScopedPtr<V>::Type m_emulatorPrecisionScaleVec;
  typename ScopedPtr<V>::Type m_observationalPrecisionShapeVec;
  typename ScopedPtr<V>::Type m_observationalPrecisionScaleVec;
  typename ScopedPtr<V>::Type m_emulatorCorrelationStrengthAlphaVec;
  typename ScopedPtr<V>::Type m_emulatorCorrelationStrengthBetaVec;
  typename ScopedPtr<V>::Type m_discrepancyPrecisionShapeVec;
  typename ScopedPtr<V>::Type m_discrepancyPrecisionScaleVec;
  typename ScopedPtr<V>::Type m_discrepancyCorrelationStrengthAlphaVec;
  typename ScopedPtr<V>::Type m_discrepancyCorrelationStrengthBetaVec;
  typename ScopedPtr<V>::Type m_emulatorDataPrecisionShapeVec;
  typename ScopedPtr<V>::Type m_emulatorDataPrecisionScaleVec;

  // The gaussian process object to build
  typename ScopedPtr<GPMSAEmulator<V, M> >::Type gpmsaEmulator;
  bool m_constructedGP;

  // Block diagonal matrix; sacrificing efficiency for clarity
  typename ScopedPtr<M>::Type m_observationErrorMatrix;

  //
  // Intermediate calculations we can cache
  //
  typename ScopedPtr<M>::Type m_TruncatedSVD_simulationOutputs;

  std::vector<M> m_discrepancyMatrices;

  typename ScopedPtr<M>::Type m_BMatrix;

  // Matrix of svd basis vectors
  typename ScopedPtr<M>::Type K;

  typename ScopedPtr<V>::Type residual;

  // Cached calculation of (B^T*W_y*B)^-1
  typename ScopedPtr<M>::Type BT_Wy_B_inv;

  // Cached calculation of (K^T*K)^-1
  typename ScopedPtr<M>::Type KT_K_inv;


private:
  bool allocated_m_opts;
  GPMSAOptions * m_opts;

};

}  // End namespace QUESO

#endif // UQ_GPMSA_HELPER_H