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