Interpolates its argument to its hypercube rank must equal arg->rank(), and xvectors type must match. More...

#include <interpolateHypercube.h>

Inheritance diagram for civita::InterpolateHC:

Collaboration diagram for civita::InterpolateHC:

Classes
struct	WeightedIndex
	structure for referring to an argument index and its weight More...

struct	WeightedIndexVector

Public Member Functions
void	setArgument (const TensorPtr &a, const ITensor::Args &args={"", 0}) override

double	operator[] (std::size_t) const override
	return or compute data at a location More...

Timestamp	timestamp () const override
	timestamp indicating how old the dependendent data might be. Used in CachedTensorOp to determine when to invalidate the cache More...

Public Member Functions inherited from civita::ITensor
	CLASSDESC_ACCESS (ITensor)

	ITensor ()

	ITensor (const Hypercube &hc)

	ITensor (Hypercube &&hc)

	ITensor (const std::vector< unsigned > &dims)

	ITensor (const ITensor &)=default

	ITensor (ITensor &&)=default

ITensor &	operator= (const ITensor &)=default

ITensor &	operator= (ITensor &&)=default

virtual	~ITensor ()

virtual const Hypercube &	hypercube () const
	information describing the axes, types and labels of this tensor More...

virtual const Hypercube &	hypercube (const Hypercube &hc)

virtual const Hypercube &	hypercube (Hypercube &&hc)

std::size_t	rank () const

std::vector< unsigned >	shape () const

void	imposeDimensions (const Dimensions &dimensions)
	impose dimensions according to dimension map dimensions More...

virtual const Index &	index () const
	the index vector - assumed to be ordered and unique More...

virtual std::size_t	size () const
	return number of elements in tensor - maybe less than hypercube.numElements if sparse More...

double	atHCIndex (std::size_t hcIdx) const
	returns the data value at hypercube index hcIdx, or NaN if More...

template<class T >
std::size_t	hcIndex (const std::initializer_list< T > &indices) const

template<class T >
double	operator() (const std::initializer_list< T > &indices) const

virtual void	setArguments (const TensorPtr &, const TensorPtr &, const ITensor::Args &args={})

virtual void	setArguments (const std::vector< TensorPtr > &a, const ITensor::Args &args={"", 0})

virtual void	setArguments (const std::vector< TensorPtr > &a1, const std::vector< TensorPtr > &a2, const ITensor::Args &args={"", 0})

Private Member Functions
void	sortAndAdd (const XVector &)

vector< std::size_t >	splitAndRotate (std::size_t) const

WeightedIndexVector	bodyCentredNeighbourhood (std::size_t idx) const
	computes the neighbourhood around a target argument index when the target index is not on the hypercube More...

Private Attributes
TensorPtr	arg

Hypercube	interimHC
	hypercube that's been rotated to match the arguments hypercube More...

std::vector< std::size_t >	strides
	strides along each dimension of this->hypercube() More...

std::vector< std::size_t >	rotation
	permutation of axes of interimHC and this->hypercube() More...

std::vector< std::pair< XVector, std::vector< std::size_t > > >	sortedArgHC

vector< WeightedIndexVector >	weightedIndices
	map from this tensor's index into the argument tensor More...

Additional Inherited Members
Public Types inherited from civita::ITensor
using	Timestamp = std::chrono::time_point< std::chrono::high_resolution_clock >

Protected Member Functions inherited from civita::ITensor
void	notImpl () const

Protected Attributes inherited from civita::ITensor
Hypercube	m_hypercube

Index	m_index

Detailed Description

Interpolates its argument to its hypercube rank must equal arg->rank(), and xvectors type must match.

Definition at line 29 of file interpolateHypercube.h.

Member Function Documentation

◆ bodyCentredNeighbourhood()

InterpolateHC::WeightedIndexVector civita::InterpolateHC::bodyCentredNeighbourhood ( std::size_t idx ) const

private

computes the neighbourhood around a target argument index when the target index is not on the hypercube

Definition at line 140 of file interpolateHypercube.cc.

References arg, civita::ITensor::index(), interimHC, civita::ITensor::rank(), anonymous_namespace{interpolateHypercube.cc}::sorted(), sortedArgHC, and splitAndRotate().

Referenced by setArgument().

   {
     // note this agorithm is limited in rank (typically 32 dims on 32bit machine, or 64 dims on 64bit)
     if (rank()>sizeof(size_t)*8)
       throw runtime_error("Ranks > "+to_string(sizeof(size_t)*8)+" not supported");
     size_t numNeighbours=size_t(1)<<rank();
     vector<size_t> iIdx=splitAndRotate(destIdx);
     const auto& argHC=arg->hypercube();
     // loop over the nearest neighbours in argument hypercube space of
     // this point in interimHypercube space
     double sumWeight=0;
     WeightedIndexVector r;
     auto argIndexVector=arg->index();
     // multivariate interpolation - eg see Abramowitz & Stegun 25.2.66
     for (size_t nbr=0; nbr<numNeighbours; ++nbr)
       {
         size_t argIdx=0;
         double weight=1;
         size_t idx=0;
         for (size_t dim=0, stride=1; dim<rank(); stride*=argHC.xvectors[dim].size(), ++dim)
           {
             const auto& x=sortedArgHC[dim].first;
             assert(!x.empty());
             assert(sorted(x.begin(),x.end()));
             auto v=interimHC.xvectors[dim][iIdx[dim]];
             auto lesserIt=std::upper_bound(x.begin(), x.end(), v, AnyLess());
             if (lesserIt!=x.begin()) --lesserIt; // find greatest value <= v, 
             boost::any lesser=*lesserIt, greater;
             if (diff(lesser, v)>=0)
               greater=lesser;  // one sided interpolation for beginning or exact match 
             else if (lesserIt+1==x.end())
               greater=x.back(); // one sided interpolation for end
             else
               greater=*(lesserIt+1);
 
             bool greaterSide = nbr&(size_t(1)<<dim); // on higher side of hypercube
             double d=diff(greater,lesser);
             if (d==0 && greaterSide) goto nextNeighbour; // already taken lesserVal at weight 1.
 
             idx += sortedArgHC[dim].second[(lesserIt-x.begin()+greaterSide)]*stride;
             assert(idx<argHC.numElements());
             if (d>0)
               weight *= greaterSide? diff(v,lesser): d-diff(v,lesser);
           }
         if (index().empty())
           {
             r.emplace_back(WeightedIndex{idx,weight});
             sumWeight+=weight;
           }
         else
           {
             auto indexV=index();
             if (binary_search(indexV.begin(), indexV.end(), idx))
               {
                 r.emplace_back(WeightedIndex{idx,weight});
                 sumWeight+=weight;
               }
           }
       nextNeighbour:;
       }
 
     // normalise weights
     for (auto& i: r) i.weight/=sumWeight;
     return r;
   }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ operator[]()

double civita::InterpolateHC::operator[] ( std::size_t ) const

overridevirtual

return or compute data at a location

Implements civita::ITensor.

Definition at line 124 of file interpolateHypercube.cc.

References arg, and weightedIndices.

   {
     if (idx<weightedIndices.size())
       {
         if (weightedIndices[idx].empty()) return nan("");
         double r=0;
         for (const auto& i: weightedIndices[idx])
           {
             assert(i.index<arg->hypercube().numElements());
             r+=i.weight * arg->atHCIndex(i.index);
           }
         return r;
       }
     return nan("");
   }

◆ setArgument()

void civita::InterpolateHC::setArgument	(	const TensorPtr &	a,
		const ITensor::Args &	args = `{"",0}`
	)

overridevirtual

Reimplemented from civita::ITensor.

Definition at line 49 of file interpolateHypercube.cc.

References arg, bodyCentredNeighbourhood(), civita::ITensor::hypercube(), civita::ITensor::index(), interimHC, civita::ITensor::rank(), rotation, civita::ITensor::size(), sortAndAdd(), strides, and weightedIndices.

   {
     arg=a;
     if (rank()!=arg->rank())
       throw runtime_error("Rank of interpolated tensor doesn't match its argument");
     // reorder hypercube for type and name
     interimHC.xvectors.clear();
     size_t stride=1;
     const auto& targetHC=hypercube().xvectors;
     rotation.clear();
     rotation.resize(rank(), rank());
     // ensure rotation vector will contain unique indices
     std::map<size_t,size_t> tmpRotation;
     // construct interim hypercube and its rotation permutation
     for (size_t i=0; i<rank(); ++i)
       {
         const auto& src=arg->hypercube().xvectors[i];
         sortAndAdd(src);
         if (src.name.empty())
           {
             const auto& dst=targetHC[i];
             if (!dst.name.empty() || src.dimension.type!=dst.dimension.type ||
                 src.dimension.units!=dst.dimension.units) //TODO handle conversion between different units
               throw runtime_error("mismatch between unnamed dimensions");
             interimHC.xvectors.push_back(dst);
             tmpRotation.emplace(make_pair(i,i));
           }
         else // find matching names dimension
           {
             auto dst=find_if(targetHC.begin(), targetHC.end(),
                              [&](const XVector& i){return i.name==src.name;});
             if (dst==targetHC.end())
             {
               // possible alternative when targetHC has no xvectors or undefined ones. for feature 147
               interimHC.xvectors.push_back(src);
               tmpRotation.emplace(make_pair(i,i));
             }
             else
               {
                 interimHC.xvectors.push_back(*dst);
                 tmpRotation.emplace(make_pair(dst-targetHC.begin(),i));
               }
           }
         strides.push_back(stride);
         stride*=targetHC[i].size();
       }
     if (tmpRotation.size()!=rank()) throw runtime_error("rotation of indices is not a permutation"); 
     for (auto& i: tmpRotation) rotation[i.first]=i.second;
 #ifndef NDEBUG
     for (auto& i: rotation) assert(i<rank()); // check that no indices have been doubly assigned.
     // Now we're sure rotation is a permutation
 #endif
     if (index().empty())
       for (size_t i=0; i<size(); ++i)
         weightedIndices.push_back(bodyCentredNeighbourhood(i));
     else
       for (auto i: index())
         weightedIndices.push_back(bodyCentredNeighbourhood(i));
   }

Here is the call graph for this function:

◆ sortAndAdd()

void civita::InterpolateHC::sortAndAdd ( const XVector & xv )

private

Definition at line 109 of file interpolateHypercube.cc.

References anonymous_namespace{interpolateHypercube.cc}::sorted(), and sortedArgHC.

Referenced by setArgument().

   {
     sortedArgHC.emplace_back();
     auto& s=sortedArgHC.back();
     for (size_t i=0; i<xv.size(); ++i)
       s.second.push_back(i);
     AnyLess less;
     sort(s.second.begin(), s.second.end(),
          [&](size_t i, size_t j){return less(xv[i],xv[j]);});
     for (auto i: s.second)
       s.first.push_back(xv[i]);
     assert(sorted(s.first.begin(), s.first.end()));
   }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ splitAndRotate()

vector< size_t > civita::InterpolateHC::splitAndRotate ( std::size_t ) const

private

Definition at line 39 of file interpolateHypercube.cc.

References civita::ITensor::hcIndex(), civita::ITensor::hypercube(), civita::ITensor::rank(), rotation, and strides.

Referenced by bodyCentredNeighbourhood().

   {
     vector<size_t> r(rank());
     const auto& h=hypercube();
     for (size_t dim=0; dim<rank(); ++dim)
       r[rotation[dim]] = (hcIndex/strides[dim]) % h.xvectors[dim].size();
     return r;
   }