1

我在写一个C++函数来计算边际PDF(概率密度函数)。这基本上意味着我可以获得沿着多个变量网格定义的多维数据(PDF)。我想将数据整合到未定义数量的维度上以保持函数的一般性。创建多个数组的任意视图

PDF的维度可以是任意的,边缘PDF的维度也可以是任意的。不可能定义输入数据的维数的顺序,所以我向该函数发送一个向量,该向量说明需要保留哪些变量。其他变量需要整合。

因此,例如: 否的变量:5(A,B,C,d,e)所示,PDF数据尺寸如图5所示,计算的边缘PDF(A,C,d)。这意味着变量/维度0,2,3需要保留,其他需要被集成(通过一个确定的积分)。 so:PDF [a] [b] [c] [d] [e] - > MARGPDF [a] [c] [d](包含其他值) for each [a] [c] [d] I需要对其他维度的数据执行操作[b] [e]。我可以通过制作视图来实现这一点,但是现在我不能动态地做到这一点。动态的我的意思是我希望维度的数量和维度可以自由选择。

基本上,我想要的是创建维度b和e中所有值的视图,并为a,c,d的每个(循环)值执行此操作。但是,我希望函数是一般的,以便输入可以是任何多数组,并且输出变量可以自由选择。所以它也可以是:PDF [a] [b] [c] [d] - > MARGPDF [c]或PDF [a] [b] [c] [d] [e] [f] - > MARGPDF [b ] [d]。

我有以下想法: 我由尺寸的PDF多数组进行排序,使得我可以创建尺寸的最后数目的视图,因此: PDF [A] [B] [C] [ d] [e]成为PDF [a] [c] [d] [b] [e]。然后,我循环遍历每个a,c,d并创建其余2个维度b和e的视图。我使用该视图执行计算并将值保存到MARGPDF [a] [c] [d]。

我需要知道的执行此操作的是: 如何切换boost :: multi_array的维度/索引的顺序? 如何在维度空闲时创建视图? 或者你有任何其他想法来完成同样的事情?

我的代码开始提供如下:

template<class DataType, int Dimension, int ReducedDimension> 
boost::multi_array<DataType, ReducedDimension> ComputeMarginalPDF(boost::multi_array<DataType, Dimension> PDF, 
                  std::vector< std::vector<DataType> > Variables , std::vector<int> VarsToKeep){ 
// check input dimensions 
if (VarsToKeep.size() != ReducedDimension){ 
    std::cout << "Dimensions do not match" << std::endl; 
} 

std::vector< std::vector<double> > NewVariables(0) ; 

// Construct reduced array with proper dimensions 
typedef boost::multi_array< DataType , ReducedDimension > ReducedArray ; 
boost::array< ReducedArray::index , ReducedDimension > dimensions; 

// get dimensions from array and insert into dimensions ; 
// set Marginal PDF dimensions 
for(int i = 0 ; i < VarsToKeep.size() ; i++){ 
    dimensions[i] = PDF.shape()[ VarsToKeep[i] ] ; 
    NewVariables.push_back(Variables[ VarsToKeep[i] ]); 
} 

ReducedArray Marginal(dimensions) ; 

// to be filled with code 

我希望我不要混淆。欢迎任何改善问题的建议。

回答

0

我修正了这个问题。我想我不能创建任意维度的boost :: multi_array,因为它需要维度作为模板参数,这在编译器时需要知道。这意味着我无法创建任意维度的视图。

因此,我做了以下操作: 我对PDF进行了排序,以便将要整合的尺寸是最后尺寸(最有可能不是最有效的方法)。 然后我逐个缩小了PDF的维数。每个循环只集成了一个维度,我保存在一个与初始数组大小相同的multi_array中(因为我无法使维度动态化)。 之后,我将这些值复制到缩小大小(已知)的multi_array中。

我独立使用下面的链接继续循环,尺寸,在尺寸:

// Dimension-independent loop over boost::multi_array?

1

几个月前我有类似的问题,但我只需要计算一维边缘。这是对我有用的解决方案的大纲,我想它也可以适应多维边缘:

我基本上是将pdf存储在一维数组/矢量中(使用任何你喜欢的):

double* pdf = new double[a*b*c*d*e]; 

然后我用您可以在二维阵列a[width][height]存储为一维阵列b[widht*height]和访问任何元素作为a[x][y]b[width*x + y]。您可以将此公式推广到任意维度,并且正确使用模/整数除法,您还可以计算反比。

使用模板,从一维索引到N维索引以及反之亦然的计算都非常简单。这使您可以将取决于维度的符号PDF[a][b][c][d][e]转换为易于扩展到任意维度的符号,如PDF(std::vector<size_t>{a,b,c,d,e}),因为您可以预先在循环中填充该向量。

如果您认为这种方法可能对您有所帮助,我可以尝试获取我的实施的一些关键功能并将其添加到此处。

编辑:

template <size_t DIM> 
inline void posToPosN(const size_t& pos, 
    const size_t* const size, 
    size_t* const posN){ 
    size_t r = pos; 

    for (size_t i = DIM; i > 0; --i){ 
     posN[i - 1] = r % size[i - 1]; 
     r /= size[i - 1]; 
    } 
} 

template <size_t DIM> 
inline void posNToPos(size_t& pos, 
    const size_t* const size, 
    const size_t* const posN){ 
    pos = 0; 
    size_t mult = 1; 

    for (size_t i = DIM; i > 0; --i){ 
     pos += mult * posN[i - 1]; 
     mult *= size[i - 1]; 
    } 
} 

template<typename type, size_t DIM> 
class Iterator{ 
private: 
    type* const _data; //pointer to start of Array 
    size_t _pos; //1-dimensional position 
    size_t _posN[DIM]; //n-dimensional position 
    size_t const * const _size; //pointer to the _size-Member of Array 
    size_t _total; 

private: 

public: 
    Iterator(type* const data, const size_t* const size, size_t total, size_t pos) 
     : _data(data), _pos(pos), _size(size), _total(total) 
    { 
     if (_pos > _total || _pos < 0) _pos = _total; 
     posToPosN<DIM>(_pos, _size, _posN); 
    } 

    bool operator!= (const Iterator& other) const 
    { 
     return _pos != other._pos; 
    } 

    type& operator*() const{ 
     if (_pos >= _total) 
      std::cout << "ERROR, dereferencing too high operator"; 
     return *(_data + _pos); 
    } 

    const Iterator& operator++() 
    { 
     ++_pos; 
     if (_pos > _total) _pos = _total; 

     posToPosN<DIM>(_pos, _size, _posN); 
     return *this; 
    } 

    Iterator& operator +=(const size_t& b) 
    { 
     _pos += b; 
     if (_pos > _total) _pos = _total; 

     posToPosN<DIM>(_pos, _size, _posN); 
     return *this; 
    } 

    const Iterator& operator--() 
    { 
     if (_pos == 0) 
      _pos = _total; 
     else 
      --_pos; 

     posToPosN<DIM>(_pos, _size, _posN); 
     return *this; 
    } 

    //returns position in n-th dimension 
    size_t operator[](size_t n){ 
     return _posN[n]; 
    } 

    //returns a new iterator, advanced by n steps in the dim Dimension 
    Iterator advance(size_t dim, int steps = 1){ 
     if (_posN[dim] + steps < 0 || _posN[dim] + steps >= _size[dim]){ 
      return Iterator(_data, _size, _total, _total); 
     } 

     size_t stride = 1; 
     for (size_t i = DIM - 1; i > dim; --i){ 
      stride *= _size[i]; 
     } 

     return Iterator(_data, _size, _total, _pos + steps*stride); 
    } 
}; 


template <typename type, size_t DIM> 
class Array{ 
    type* _data; 
    size_t _size[DIM]; 
    size_t _total; 

    void init(const size_t* const dimensions){ 
     _total = 1; 
     for (int i = 0; i < DIM; i++){ 
      _size[i] = dimensions[i]; 
      _total *= _size[i]; 
     } 

     _data = new type[_total]; 
    } 

public: 
    Array(const size_t* const dimensions){ 
     init(dimensions); 
    } 

    Array(const std::array<size_t, DIM>& dimensions){ 
     init(&dimensions[0]); 
    } 

    ~Array(){ 
     delete _data; 
    } 
    Iterator<type, DIM> begin(){ 
     return Iterator<type, DIM>(_data, _size, _total, 0); 
    } 
    Iterator<type, DIM> end(){ 
     return Iterator<type, DIM>(_data, _size, _total, _total); 
    } 
    const size_t* const size(){ 
     return _size; 
    } 
}; 


//for projections of the PDF 
void calc_marginals(size_t dir, double* p_xPos, double* p_yPos){ 
    assert(dir < N_THETA); 

    std::lock_guard<std::mutex> lock(calcInProgress); 

    //reset to 0 
    for (size_t i = 0; i < _size[dir]; ++i){ 
     p_yPos[i] = 0; 
    } 

    //calc projection 
    double sum = 0; 
    for (auto it = _p_theta.begin(); it != _p_theta.end(); ++it){ 

     p_yPos[it[dir]] += (*it); 
     sum += (*it); 
    } 

    if (abs(sum - 1) > 0.001){ cout << "Warning: marginal[" << dir << "] not normalized" << endl; } 
    //calc x-Axis 
    for (size_t i = 0; i < _size[dir]; ++i){ 
     p_xPos[i] = _p[dir].start + double(i)/double(_size[dir] - 1)*(_p[dir].stop - _p[dir].start); 
    } 
} 

该代码由几部分组成:

  • 两个函数posToPosN()posNToPos()其中做一维和DIM维坐标之间的变换提到。尺寸在这里作为模板参数DIM给出。 pos只是一维位置,posN一个指向DIM30维坐标的尺寸为​​的数组的指针,而size是一个尺寸为​​的数组,其宽度在不同的方向上(在您的情况下,{a,b, c,d,e})
  • Iterator是一个Iterator类,允许在DIM-Dimensional Array上使用基于范围或基于迭代器的for-loops。注意它返回DIM维坐标的第n个分量和advance()函数返回一个迭代器具有坐标的元素{posN[0], posN[1], ...,posN[dim] + steps , ... posN[DIM]}
  • Array应该是相当简单明了
  • calcMarginalsoperator[](size_t n)为i用于计算边缘人的函数。 dir这里是我想要计算边际的方向(记住:一维边际),并写入p_xPosp_yPos_p_thetaArray。注意基于迭代器的for循环,(*it)这里指的是存储在数组中的pdf的double值,就像通常的迭代器一样。另外,it[dim]返回方向dim中实际值的坐标。写入p_xPos的最后一个循环仅仅是因为我不想在这个数组中使用索引而是实际值。

我想如果你重新定义Iterator::operator[]采取维度索引的矢量/阵列并返回相应坐标的矢量/阵列并添加Array::operator[]随机访问,这需要矢量/阵列,以及你应该几乎完成了。

+0

谢谢您的回复!我认为这种方法可能会有所帮助。你能提供这些关键功能吗? – Rene

+0

@Rene将我的代码添加到它,希望它有帮助 – Anedar