1/*! 2 * \file include/TFEL/Math/Kriging/FactorizedKriging.ixx 3 * \brief 4 * \author Thomas Helfer 5 * \brief 12 avr 2009 6 * \copyright Copyright (C) 2006-2018 CEA/DEN, EDF R&D. All rights 7 * reserved. 8 * This project is publicly released under either the GNU GPL Licence 9 * or the CECILL-A licence. A copy of thoses licences are delivered 10 * with the sources of TFEL. CEA or EDF may also distribute this 11 * project under specific licensing conditions. 12 */ 13 14#ifndef LIB_TFEL_MATH_FACTORIZEDKRIGINGIXX 15#define LIB_TFEL_MATH_FACTORIZEDKRIGINGIXX 16 17#include<algorithm> 18 19#include"TFEL/Math/matrix.hxx" 20#ifdef HAVE_ATLAS 21#include"TFEL/Math/Bindings/atlas.hxx" 22#else 23#include"TFEL/Math/LUSolve.hxx" 24#endif /* LIB_TFEL_MATH_FACTORIZEDKRIGINGIXX */ 25#include"TFEL/Math/Kriging/KrigingErrors.hxx" 26 27namespace tfel 28{ 29 30 namespace math 31 { 32 33 template<unsigned short N,unsigned short M, 34 typename T,typename Model1,typename Model2> 35 FactorizedKriging<N,M,T,Model1,Model2>::FactorizedKriging(const Model1& m1_, 36 const Model2& m2_) 37 : m1(m1_),m2(m2_) 38 {} 39 40 template<unsigned short N,unsigned short M, 41 typename T,typename Model1,typename Model2> 42 T 43 FactorizedKriging<N,M,T,Model1,Model2>::operator()(const typename KrigingVariable<N,T>::type& xv1, 44 const typename KrigingVariable<M,T>::type& xv2) const 45 { 46 using namespace tfel::math::internals; 47 typedef typename vector<T>::difference_type diff; 48 auto p = a.begin()+static_cast<diff>(this->x1.size()); 49 T r(0); 50 for(typename vector<T>::size_type i=0;i!=x1.size();++i){ 51 r+=a[i]*(m1.covariance(xv1-this->x1[i])*m2.covariance(xv2-this->x2[i])); 52 } 53 ApplySpecificationDrifts<0,Model1::nb,N,T,Model1>::apply(r,p,xv1); 54 ApplySpecificationDrifts<0,Model2::nb,M,T,Model2>::apply(r,p,xv2); 55 return r; 56 } // end of FactorizedKriging<N,M,T,Model>::operator() 57 58 template<unsigned short N,unsigned short M, 59 typename T,typename Model1,typename Model2> 60 void 61 FactorizedKriging<N,M,T,Model1,Model2>::addValue(const typename KrigingVariable<N,T>::type& xv1, 62 const typename KrigingVariable<M,T>::type& xv2, 63 const T& fv) 64 { 65 this->x1.push_back(xv1); 66 this->x2.push_back(xv2); 67 this->f.push_back(fv); 68 } 69 70 template<unsigned short N,unsigned short M, 71 typename T,typename Model1,typename Model2> 72 void 73 FactorizedKriging<N,M,T,Model1,Model2>::buildInterpolation() 74 { 75 using namespace std; 76 using namespace tfel::math::internals; 77#ifdef HAVE_ATLAS 78 using namespace tfel::math::atlas; 79#endif /* LIB_TFEL_MATH_FACTORIZEDKRIGINGIXX */ 80 using tfel::math::vector; 81 if((x1.size()!=f.size())||(x2.size()!=f.size())){ 82 throw(KrigingErrorInvalidLength()); 83 } 84 if((x1.empty())||(x2.empty())){ 85 throw(KrigingErrorNoDataSpecified()); 86 } 87 if((x1.size()<=Model1::nb)||(x2.size()<=Model2::nb)){ 88 throw(KrigingErrorInsufficientData()); 89 } 90 matrix<T> m(x1.size()+Model1::nb+Model2::nb, 91 x1.size()+Model1::nb+Model2::nb,T(0)); // temporary matrix 92 this->a.resize(x1.size()+Model1::nb+Model2::nb,T(0)); // unknowns 93 copy(this->f.begin(),this->f.end(),this->a.begin()); // copy values 94 // filling the matrix 95 for(typename vector<T>::size_type i=0;i!=x1.size();++i){ 96 for(typename vector<T>::size_type j=0;j!=i;++j){ 97 m(i,j) = m(j,i) = m1.covariance(this->x1[i]-this->x1[j])*m2.covariance(this->x2[i]-this->x2[j]); 98 } 99 m(i,i) = T(0); 100 } 101 ApplySpecificationDrifts<0,Model1::nb,N,T,Model1>::apply(m,this->x1.size(),this->x1); 102 ApplySpecificationDrifts<0,Model2::nb,M,T,Model2>::apply(m,this->x1.size()+Model1::nb,this->x2); 103 // computing unknown values 104#ifdef HAVE_ATLAS 105 gesv(m,this->a); 106#else 107 LUSolve::exe(m,this->a); 108#endif /* LIB_TFEL_MATH_FACTORIZEDKRIGINGIXX */ 109 } 110 111 template<unsigned short N,unsigned short M, 112 typename T,typename Model1,typename Model2> 113 FactorizedKriging<N,M,T,Model1,Model2>::~FactorizedKriging() noexcept 114 {} 115 116 } // end of namespace math 117 118} // end of namespace tfel 119 120#endif /* LIB_TFEL_MATH_FACTORIZEDKRIGINGIXX */ 121 122