gsSpaceTimeFitter.h

Go to the documentation of this file.

#pragma once

namespace gismo
{

template<size_t domainDim, class T>
class gsSpaceTimeFitter
{
    typedef typename gsTensorBSpline<domainDim+1,real_t>::BoundaryGeometryType slice_t;


public:
  gsSpaceTimeFitter  ( const std::vector<gsMatrix<T>> & solutionCoefs,
                      const gsVector<T> & times,
                      const gsVector<T> & ptimes,
                      const gsMultiBasis<T> & spatialBasis,
                      const index_t deg = 2)
  :
  m_data(solutionCoefs),
  m_times(times),
  m_ptimes(ptimes),
  m_bases(spatialBasis),
  m_deg(deg)
  {
    GISMO_ENSURE(m_data.size() != 0,"No data provided!");
    GISMO_ENSURE(m_data.size() == (size_t)m_times.size(),"Solution coefs and times should have the same size");
    GISMO_ENSURE((size_t)m_ptimes.size() == (size_t)m_times.size(),"(Parametric)times should have the same size");
    m_targetDim = m_data.at(0).cols();
  }

  void setDegree(index_t deg) {m_deg = deg;}

  void addDataPoint(gsMatrix<T> & solution, T time, T ptime, index_t continuity )
  {
    index_t N = m_ptimes.rows();
    m_ptimes.conservativeResize(N+1);
    m_ptimes.row(N) = ptime;
    gsDebugVar(m_ptimes);
    std::sort(begin(m_ptimes),end(m_ptimes));
    gsDebugVar(m_ptimes);
  }

protected:

  gsTensorBSplineBasis<domainDim+1,T> _basis(const gsKnotVector<T> &kv, index_t nsteps)
  {
    return _basis_impl<domainDim+1>(kv,nsteps);
  }

  template<index_t _domainDim>
  typename std::enable_if<_domainDim==2, gsTensorBSplineBasis<_domainDim,T>>::type
  _basis_impl(const gsKnotVector<T> &kv, index_t nsteps)
  {
    gsTensorBSplineBasis<_domainDim,T> tbasis(
                                            static_cast<gsBSplineBasis<T> *>(&m_bases.basis(0).component(1))->knots(),
                                            kv
                                            );
    return tbasis;
  }

  template<index_t _domainDim>
  typename std::enable_if<_domainDim==3, gsTensorBSplineBasis<_domainDim,T>>::type
  _basis_impl(const gsKnotVector<T> &kv, index_t nsteps)
  {
    gsTensorBSplineBasis<_domainDim,T> tbasis(
                                            static_cast<gsBSplineBasis<T> *>(&m_bases.basis(0).component(0))->knots(),
                                            static_cast<gsBSplineBasis<T> *>(&m_bases.basis(0).component(1))->knots(),
                                            kv
                                            );
    return tbasis;
  }

  template<index_t _domainDim>
  typename std::enable_if<_domainDim==4, gsTensorBSplineBasis<_domainDim,T>>::type
  _basis_impl(const gsKnotVector<T> &kv, index_t nsteps)
  {
    gsTensorBSplineBasis<_domainDim,T> tbasis(
                                            static_cast<gsBSplineBasis<T> *>(&m_bases.basis(0).component(0))->knots(),
                                            static_cast<gsBSplineBasis<T> *>(&m_bases.basis(0).component(1))->knots(),
                                            static_cast<gsBSplineBasis<T> *>(&m_bases.basis(0).component(2))->knots(),
                                            kv
                                            );
    return tbasis;
  }

public:

  std::pair<T,gsGeometry<T> *> slice(T xi)
  {
    slice_t res;
    m_fit.slice(domainDim,xi,res);
    gsGeometry<T> * geom = res.clone().release();
    T load = geom->coefs()(0,domainDim+1);
    geom->embed(m_targetDim);
    return std::make_pair(load,geom);
  }

  void compute()
  {
    GISMO_ASSERT(m_data.size()==(size_t)m_times.rows(),"Number of time and solution steps should match! "<<m_data.size()<<"!="<<m_times.rows());
    // GISMO_ASSERT(m_data.at(0).cols() == dim,"Is the dimension correct?"<<m_data.at(0).cols()<<"!="<<dim);
    index_t nsteps = m_times.rows();
    index_t bsize = m_data.at(0).rows();

    // Prepare fitting basis
    gsKnotVector<> kv(m_ptimes.minCoeff(),m_ptimes.maxCoeff(),nsteps-(m_deg+1),m_deg+1);
    gsBSplineBasis<T> lbasis(kv);


    m_basis = _basis(kv,nsteps);

    gsMatrix<T> rhs(m_times.size(),m_targetDim*bsize+1);
    gsVector<T> ones; ones.setOnes(bsize);

    for (index_t lam = 0; lam!=nsteps; ++lam)
    {
      rhs.block(lam,0,1,m_targetDim * bsize) = m_data.at(lam).reshape(1,m_targetDim * bsize);
      rhs(lam,m_targetDim*bsize) = m_times.at(lam);
    }

    // get the Greville Abcissae (anchors)
    gsMatrix<T> anchors = lbasis.anchors();

    // Get the collocation matrix at the anchors
    gsSparseMatrix<T> C = lbasis.collocationMatrix(anchors);

    gsSparseSolver<>::LU solver;
    solver.compute(C);

    gsMatrix<T> sol;
    m_coefs.resize((nsteps)*bsize,m_targetDim+1);

    sol = solver.solve(rhs);

    for (index_t lam = 0; lam!=nsteps; ++lam)
    {
      gsMatrix<> tmp = sol.block(lam,0,1,m_targetDim * bsize);
      m_coefs.block(lam * bsize,0,bsize,m_targetDim) = tmp.reshape(bsize,m_targetDim);
      m_coefs.block(lam * bsize,m_targetDim,bsize,1) = sol(lam,m_targetDim*bsize) * ones;
    }

    // gsTensorBSpline<3,T> tspline = tbasis.makeGeometry(give(coefs)).release();
    // gsTensorBSpline<dim,T> tspline(m_basis,give(m_coefs));

    m_fit = gsTensorBSpline<domainDim+1,T>(m_basis,give(m_coefs));
  }

protected:
  index_t m_targetDim;
  std::vector<gsMatrix<T>> m_data;
  gsVector<T> m_times;
  gsVector<T> m_ptimes;
  gsMultiBasis<T> m_bases;
  index_t m_deg;

  mutable gsTensorBSplineBasis<domainDim+1,T> m_basis;
  gsMatrix<T> m_coefs;

  mutable gsTensorBSpline<domainDim+1,T> m_fit;

};


} // namespace gismo

// #ifndef GISMO_BUILD_LIB
// #include GISMO_HPP_HEADER(gsSpaceTimeHierarchy.hpp)
// #endif