PCSet.h

Go to the documentation of this file.

/* =====================================================================================
 
                      The UQ Toolkit (UQTk) version 3.1.1
                          Copyright (2021) NTESS
                        https://www.sandia.gov/UQToolkit/
                        https://github.com/sandialabs/UQTk
 
     Copyright 2021 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
     Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government
     retains certain rights in this software.
 
     This file is part of The UQ Toolkit (UQTk)
 
     UQTk is open source software: you can redistribute it and/or modify
     it under the terms of BSD 3-Clause License
 
     UQTk 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
     BSD 3 Clause License for more details.
 
     You should have received a copy of the BSD 3 Clause License
     along with UQTk. If not, see https://choosealicense.com/licenses/bsd-3-clause/.
 
     Questions? Contact the UQTk Developers at <uqtk-developers@software.sandia.gov>
     Sandia National Laboratories, Livermore, CA, USA
===================================================================================== */
 
#ifndef PCSET_H_SEEN
#define PCSET_H_SEEN
 
#include <iostream>
#include <sstream>
#include <string.h>
#include <stdio.h>
#include <map>
#include "Array1D.h"
#include "Array2D.h"
#include "error_handlers.h"
#include "ftndefs.h"
#include "quad.h"
 
/*  CVODE headers  */
#include <cvode/cvode.h>             /* prototypes for CVODE fcts., consts.  */
#include <cvode/cvode_direct.h>        /* access to CVDls interface            */
#include <nvector/nvector_serial.h>  /* serial N_Vector types, fcts., macros */
#include <sunmatrix/sunmatrix_dense.h> /* access to dense SUNMatrix            */
#include <sunlinsol/sunlinsol_dense.h> /* access to dense SUNLinearSolver      */
// #include <cvode/cvode_dense.h>       /* prototype for CVDense                */
#include <sundials/sundials_dense.h> /* definitions DlsMat DENSE_ELEM        */
#include <sundials/sundials_types.h> /* definition of type realtype          */
 
#include <iostream>
#include <string.h>
#include <stdio.h>
#include <sstream>
using namespace std; // needed for python string conversion
 
class PCBasis;
class Quad;
 
typedef enum {TaylorSeries=0, Integration} LogCompMethod;
 
 
class PCSet {
public:
 
  PCSet(const string sp_type, const int order, const int n_dim, const string pc_type,
        const double alpha=0.0, const double betta=1.0);
 
 
  PCSet(const string sp_type, const int order, const int n_dim, const string pc_type, const string pc_seq,
        const double alpha=0.0, const double betta=1.0);
 
 
  PCSet(const string sp_type, const Array1D<int>& maxOrders, const int n_dim, const string pc_type,
        const double alpha=0.0, const double betta=1.0);
 
  PCSet(const string sp_type, const Array2D<int>& customMultiIndex, const string pc_type,
        const double alpha=0.0, const double betta=1.0);
 
  ~PCSet();
 
 
  void dPhi_alpha(Array1D<double>& x, Array1D<int>& alpha, Array1D<double>& grad);
  void dPhi(Array1D<double>& x, Array2D<int>& mindex, Array1D<double>& grad, Array1D<double>& ck);
  void dPhi(Array2D<double>& x, Array2D<int>& mindex, Array2D<double>& grad, Array1D<double>& ck);
 
  void ddPhi_alpha(Array1D<double>& x, Array1D<int>& alpha, Array2D<double>& grad);
  void ddPhi(Array1D<double>& x, Array2D<int>& mindex, Array2D<double>& grad, Array1D<double>& ck);
 
 
  void SetQd1d(Array1D<double>& qdpts1d,Array1D<double>& wghts1d, int nqd);
 
  void SetQuadRule(const string grid_type,const string fs_type,int param);
 
  void SetQuadRule(Quad &quadRule);
 
 
  void PrintMultiIndex() const;
 
  void PrintMultiIndexNormSquared() const;
 
 
 
  string GetPCType() const {return pcType_;}
 
  double GetAlpha() const {return alpha_;}
 
  double GetBeta() const {return beta_;}
 
  void GetMultiIndex(Array2D<int> &mindex) const {mindex=multiIndex_;}
 
  void GetMultiIndex(int *mindex) const;
 
  void GetNormSq(Array1D<double>& normsq) const {normsq=psiSq_;}
 
  int GetNumberPCTerms() const {return nPCTerms_;}
 
  int GetNDim() const {return nDim_;}
 
  int GetOrder() const {return order_;}
 
  int GetNQuadPoints() const {return nQuadPoints_;}
 
  void GetQuadPoints(Array2D<double>& quad) const {quad=quadPoints_;}
 
  void GetQuadPointsWeights(Array2D<double>& quad, Array1D<double>& wghts) const { quad=quadPoints_; wghts=quadWeights_;}
 
  void GetQuadPoints(double* quad) const {for(int i=0;i<nQuadPoints_;i++) for(int j=0;j<nDim_;j++) quad[i*nDim_+j]=quadPoints_(i,j);}
 
  void GetQuadWeights(Array1D<double>& wghts) const {wghts=quadWeights_;}
 
  void GetQuadWeights(double* wghts) const {for(int i=0;i<nQuadPoints_;i++) wghts[i]=quadWeights_(i);}
 
  void GetPsi(Array2D<double>& psi) const {psi=psi_;}
 
  void GetPsi(double* psi) const  {for(int i=0;i<nQuadPoints_;i++) for(int j=0;j<nPCTerms_;j++) psi[i*nPCTerms_+j]=psi_(i,j);}
 
  void GetPsiSq(Array1D<double>& psisq) const {psisq=psiSq_;}
 
  void GetPsiSq(double* psisq) const {for(int i=0;i<nPCTerms_;i++) psisq[i]=psiSq_(i);}
 
  double GetTaylorTolerance() const {return rTolTaylor_;}
 
  void SetTaylorTolerance(const double& rTol) {rTolTaylor_ = rTol;}
 
  int GetTaylorTermsMax() const {return maxTermTaylor_;}
 
  void SetTaylorTermsMax(const int& maxTerm) {maxTermTaylor_ = maxTerm;}
 
  void SetLogCompMethod(const LogCompMethod& logMethod) {logMethod_ = logMethod;}
 
  double GetGMRESDivTolerance() const {return rTolGMRESDiv_;}
 
  void SetGMRESDivTolerance(const double& rTol) {rTolGMRESDiv_ = rTol;}
 
 
 
  void InitMeanStDv(const double& m, const double& s, double* p) const;
 
  void InitMeanStDv(const double& m, const double& s, Array1D<double>& p) const;
 
  void Copy(double* p1, const double* p2) const;
 
  void Copy(Array1D<double>& p1, const Array1D<double>& p2) const;
 
  void Add(const double* p1, const double* p2, double* p3) const;
 
  void Add(const Array1D<double>& p1, const Array1D<double>& p2, Array1D<double>& p3) const;
 
  void AddInPlace(double* p1, const double* p2) const;
 
  void AddInPlace(Array1D<double>& p1, const Array1D<double>& p2) const;
 
  void Multiply(const double* p1, const double& a, double* p2) const;
 
  void Multiply(const Array1D<double>& p1, const double& a, Array1D<double>& p2) const;
 
  void MultiplyInPlace(double* p1, const double& a) const;
 
  void MultiplyInPlace(Array1D<double>& p1, const double& a) const;
 
  void Subtract(const double* p1, const double* p2, double* p3) const;
 
  void Subtract(const Array1D<double>& p1, const Array1D<double>& p2, Array1D<double>& p3) const;
 
  void SubtractInPlace(double* p1, const double* p2) const;
 
  void SubtractInPlace(Array1D<double>& p1, const Array1D<double>& p2) const;
 
  void Prod(const double* p1, const double* p2, double* p3) const;
 
  void Prod(const Array1D<double>& p1, const Array1D<double>& p2, Array1D<double>& p3) const;
 
  void Prod3(const double* p1, const double* p2, const double* p3, double* p4) const;
 
  void Prod3(const Array1D<double>& p1, const Array1D<double>& p2, const Array1D<double>& p3,
            Array1D<double>& p4) const;
 
  void Polyn(const double* polycf, int npoly, const double* p1, double* p2) const;
 
  void Polyn(const Array1D<double>& polycf, const Array1D<double>& p1, Array1D<double>& p2) const;
 
  void PolynMulti(const Array1D<double>& polycf, const Array2D<int>& mindex, const Array2D<double>& p1, Array1D<double>& p2) const;
 
  void Exp(const double* p1, double* p2) const;
 
  void Exp(const Array1D<double>& p1, Array1D<double>& p2) const;
 
  void Log(const double* p1, double* p2) const;
 
  void Log(const Array1D<double>& p1, Array1D<double>& p2) const;
 
  void Log10(const double* p1, double* p2) const;
 
  void Log10(const Array1D<double>& p1, Array1D<double>& p2) const;
 
  void RPow(const double* p1, double* p2, const double& a) const;
 
  void RPow(const Array1D<double>& p1, Array1D<double>& p2, const double& a) const;
 
  void IPow(const double* p1, double* p2, const int& ia) const;
 
  void IPow(const Array1D<double>& p1, Array1D<double>& p2, const int& ia) const;
 
  void Inv(const double* p1, double* p2) const;
 
  void Inv(const Array1D<double>& p1, Array1D<double>& p2) const;
 
  void Div(const double* p1, const double* p2, double* p3) const;
 
  void Div(const Array1D<double>& p1, const Array1D<double>& p2, Array1D<double>& p3) const;
 
  double StDv(const double* p) const;
 
  double StDv(const Array1D<double>& p) const;
 
  double GetModesRMS(const double* p) const;
 
  double GetModesRMS(const Array1D<double>& p) const;
 
  void Derivative(const double* p1, double* p2) const;
 
  void Derivative(const Array1D<double>& p1, Array1D<double>& p2) const;
 
  int GetNumTripleProd() const;
  void GetTripleProd(int *nTriple, int *iProd, int *jProd, double *Cijk) const;
  void GetTripleProd(Array1D<int>& nTriple, Array1D<int>& iProd, Array1D<int>& jProd, Array1D<double>& Cijk) const;
  int  GetNumQuadProd() const;
  void GetQuadProd(int *nQuad, int *iProd, int *jProd, int *kProd, double *Cijkl) const;
  void GetQuadProd(Array1D<int> &nQuad, Array1D<int> &iProd, Array1D<int> &jProd, Array1D<int> &kProd,
                   Array1D<double> &Cijkl) const;
 
 
  void SeedBasisRandNumGen(const int& seed) const;
 
  void DrawSampleSet(const Array1D<double>& p, Array1D<double>& samples);
 
  void DrawSampleSet(const double* p, double* samples, const int& nSamples);
 
  void DrawSampleVar(Array2D<double>& samples) const;
  void DrawSampleVar(double *samples, const int &nS, const int &nD) const;
 
 
  double EvalPC(const Array1D<double>& p, Array1D<double>& randVarSamples);
 
  double EvalPC(const double* p, const double* randVarSamples);
 
  void EvalPCAtCustPoints(Array1D<double>& xch, Array2D<double>& custPoints,Array1D<double>& p);
 
  void EvalBasisAtCustPts(const Array2D<double>& custPoints,Array2D<double>& psi);
 
  void EvalBasisAtCustPts(const double* custPoints, const int npts, double* psi);
  // void EvalBasisAtCustPts(const int npts, const int ndim, const int npc, const double *custPoints, double *psi);
 
 
  void GalerkProjection(const Array1D<double>& fcn, Array1D<double>& ck);
 
  void GalerkProjectionMC(const Array2D<double>& x, const Array1D<double>& fcn, Array1D<double>& ck);
 
 
  int ComputeOrders(Array1D<int>& orders);
 
  int ComputeEffDims(int *effdim);
 
  int ComputeEffDims(Array1D<int> &effdim);
 
  void EncodeMindex(Array1D< Array2D<int> >& sp_mindex);
 
 
  double ComputeMean(const double *coef);
 
  double ComputeMean(Array1D<double>& coef);
 
  double ComputeVarFrac(const double *coef, double *varfrac);
 
  double ComputeVarFrac(Array1D<double>& coef, Array1D<double>& varfrac);
 
  void ComputeMainSens(Array1D<double>& coef, Array1D<double>& mainsens);
 
  void ComputeTotSens(Array1D<double>& coef, Array1D<double>& totsens);
 
  void ComputeJointSens(Array1D<double>& coef, Array2D<double>& jointsens);
 
 
 
  void SetVerbosity(int verbosity) { uqtkverbose_ = verbosity; }
 
  void EvalNormSq(Array1D<double>& normsq);
  void EvalNormSq(double* normsq, const int npc);
 
  void EvalNormSqExact(Array1D<double>& normsq);
 
  bool IsInDomain(double x);
 
 
 
 private:
  PCSet(): order_(0), nDim_(0) {};
 
 
  PCSet(const PCSet &obj):order_(obj.order_), nDim_(obj.nDim_) {};
 
  void ComputeMaxOrdPerDim();
 
  void Initialize(const string ordertype);
 
  //  void InitQuadrature();
 
  void InitISP();
  void InitNISP();
 
  void EvalBasisProd3();
 
  void EvalBasisProd4();
 
  static void GMRESMatrixVectorProdWrapper(int* n, double* x, double* y, int* nelt,
                                    int* ia, int* ja, double* a, int* obj) {
    // Look up *obj in the map that relates integer indices to pointers to PCSet
    OMap_t::iterator it = omap_->find(*obj);
    if(it == omap_->end()) {
      string err_message = (string) "GMRESMatrixVectorProdWrapper():"
                           + " the callback object is not a valid entry in the map";
      throw Tantrum(err_message);
    }
    // Perform callback to the member function of the proper PCSet instance
    it->second->GMRESMatrixVectorProd(x, a, y);
 
    return;
  }
 
  static void GMRESPreCondWrapper(int* n, double* r, double* z, int* nelt,
                                  int* ia, int* ja, double* a, int* obj,
                                  double* rwork, int* iwork) { };
 
  void GMRESMatrixVectorProd(const double* x, const double*a, double* y) const;
 
  void LogTaylor(const double* p1, double* p2) const;
 
  void LogInt(const double* p1, double* p2) const;
 
  static int LogIntRhsWrapper(realtype t, N_Vector y, N_Vector ydot, void *f_data)
  {
    double indxobj = ((double*) f_data)[0] ;
 
    OMap_t::iterator it = omap_->find((int) indxobj);
 
    if (it == omap_->end())
    {
      string err_message = (string) "LogIntRhsWrapper():"
                           + " the callback object is not a valid entry in the map";
      throw Tantrum(err_message);
    }
 
    // Perform callback to the member function of the proper PCSet instance
    it->second->LogIntRhs(t,y,ydot,f_data);
 
    return ( 0 ) ;
 
  }
 
  int LogIntRhs(realtype t, N_Vector y, N_Vector ydot, void *f_data) const;
 
  int uqtkverbose_;
 
  string spType_;
 
  string pcType_;
 
  string pcSeq_;
 
  PCBasis* p_basis_;
 
  int order_;
 
  int maxorddim_;
 
  Array1D<int> maxOrders_;
 
  Array1D<int> maxOrdPerDim_;
 
  const int nDim_;
 
  int nQuadPoints_;
 
  int nPCTerms_;
 
  double rTolTaylor_;
 
  int maxTermTaylor_;
 
  double SMALL_;
 
  double rTolGMRESDiv_;
 
  Array2D<double> psi_;
 
  Array1D<double> psiSq_;
 
  Array2D<double> quadPoints_;
 
  Array1D<double> quadWeights_;
 
  Array2D<int> quadIndices_;
 
  Array2D<int> multiIndex_;
 
  Array1D<Array1D<int> > iProd2_;
 
  Array1D<Array1D<int> > jProd2_;
 
  Array1D<Array1D<double> > psiIJKProd2_;
 
  Array1D<Array1D<int> > iProd3_;
 
  Array1D<Array1D<int> > jProd3_;
 
  Array1D<Array1D<int> > kProd3_;
 
  Array1D<Array1D<double> > psiIJKLProd3_;
 
  LogCompMethod logMethod_ ;
 
  int CVmaxord_;
 
  int CVmaxnumsteps_  ;
 
  double CVinitstep_;
 
  double CVmaxstep_;
 
  double CVrelt_;
 
  double CVabst_  ;
 
  int Check_CVflag(void *flagvalue, const char *funcname, int opt) const;
 
  int my_index_;
 
  int narg_;
 
  double alpha_;
  double beta_;
 
  typedef std::map<int, PCSet*> OMap_t;
  static int next_index_;
  static OMap_t *omap_;
 
};
 
#endif /* !PCSET_H_SEEN */