创建多阵列的任意视图
Creating arbitrary views of multi array
我正在编写一个 c++ 函数来计算边际 PDF（概率密度函数）。这基本上意味着我获得了沿多个变量的网格定义的多维数据 (PDF)。我想在未定义的维度上集成数据以保持功能通用。
PDF的尺寸可以任意，边缘PDF的尺寸也可以任意。无法定义输入数据的维度顺序，因此我向函数发送了一个向量，其中说明需要保留哪些变量。其他变量需要整合。
例如：
变量数：5(a,b,c,d,e)，PDF 数据维度 5，计算 (a,c,d) 的边际 PDF。这意味着 variables/dimensions 0,2,3 需要保留，其他需要积分出来（通过定积分）。
所以：PDF[a][b][c][d][e] -> MARGPDF[a][c][d]（包含其他值）
对于每个 [a][c][d]，我需要对其他维度 [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 的 dimensions/indices 的顺序？
尺寸自由时如何创建视图？
或者你有什么其他的想法来完成同样的事情吗？
我的代码开头如下：
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
希望我没有混淆。欢迎提出任何改进问题的建议。
几个月前我遇到了类似的问题，但我只需要计算一维边缘。这是对我有用的解决方案的概述，我想它也可以适用于多维边缘：
我基本上将 pdf 存储在一维 array/vector 中（使用任何你喜欢的）：
double* pdf = new double[a*b*c*d*e];
然后我用到你可以将二维数组 a[width][height] 存储为一维数组 b[widht*height] 并访问任何元素 a[x][y] 作为 b[width*x + y]。您可以将此公式推广到任意维度，并且通过正确使用 modulo/integer 除法，您还可以计算倒数。
从一维索引到 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 指向大小为 DIM 的数组的指针，指的是 DIM 维坐标，size 是大小为 [=22] 的数组=] 包含不同方向的宽度（在你的情况下类似于 {a,b,c,d,e}）
Iterator 是一个迭代器 class，允许在 DIM 维数组上进行基于范围或基于迭代器的 for 循环。请注意 operator[](size_t n) return 是 DIM 维坐标的第 n 个分量和 advance() 函数 return 是坐标为 [=28 的元素的迭代器=]Array 应该很直接calcMarginals 是我用来计算边缘的函数。 dir这里是我要计算边际的方向（记住：一维边际）写到p_xPos和p_yPos，_p_theta是一个Array.注意基于迭代器的 for 循环，(*it) 在这里指的是存储在数组中的 pdf 的双精度值，就像通常的迭代器一样。另外，it[dim]returns是dim方向实际值的坐标。写入 p_xPos 的最后一个循环只是因为我不想要这个数组中的索引而是实际值。
我想如果你重新定义 Iterator::operator[] 以获取维度索引的 vector/array 和 return 适当坐标的 vector/array 并添加 Array::operator[] 对于需要 vector/array 的随机访问，你应该已经完成了。
我解决了这个问题。我想我不能创建任意维度的 boost::multi_array，因为它需要维度作为模板参数，这需要在编译器时知道。这意味着我无法创建任意维度的视图。
因此，我做了以下事情：
我对 PDF 进行了排序，以便将要整合的维度是最后一个维度（很可能不是最有效的方法）。
然后我一个一个地缩小了PDF的尺寸。每个循环我只集成了 1 个维度，我将其保存在一个 multi_array 中，它与初始数组具有相同的大小（因为我无法使维度动态化）。
之后，我将值复制到缩小尺寸（已知）的 multi_array。
我使用以下 link 在维度上独立标注循环：
// Dimension-independent loop over boost::multi_array?
我正在编写一个 c++ 函数来计算边际 PDF（概率密度函数）。这基本上意味着我获得了沿多个变量的网格定义的多维数据 (PDF)。我想在未定义的维度上集成数据以保持功能通用。
PDF的尺寸可以任意，边缘PDF的尺寸也可以任意。无法定义输入数据的维度顺序，因此我向函数发送了一个向量，其中说明需要保留哪些变量。其他变量需要整合。
例如： 变量数：5(a,b,c,d,e)，PDF 数据维度 5，计算 (a,c,d) 的边际 PDF。这意味着 variables/dimensions 0,2,3 需要保留，其他需要积分出来（通过定积分）。 所以：PDF[a][b][c][d][e] -> MARGPDF[a][c][d]（包含其他值） 对于每个 [a][c][d]，我需要对其他维度 [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 的 dimensions/indices 的顺序？ 尺寸自由时如何创建视图？ 或者你有什么其他的想法来完成同样的事情吗？
我的代码开头如下：
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
希望我没有混淆。欢迎提出任何改进问题的建议。
几个月前我遇到了类似的问题，但我只需要计算一维边缘。这是对我有用的解决方案的概述，我想它也可以适用于多维边缘：
我基本上将 pdf 存储在一维 array/vector 中（使用任何你喜欢的）：
double* pdf = new double[a*b*c*d*e];
然后我用到你可以将二维数组 a[width][height] 存储为一维数组 b[widht*height] 并访问任何元素 a[x][y] 作为 b[width*x + y]。您可以将此公式推广到任意维度，并且通过正确使用 modulo/integer 除法，您还可以计算倒数。
从一维索引到 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指向大小为DIM的数组的指针，指的是 DIM 维坐标，size是大小为 [=22] 的数组=] 包含不同方向的宽度（在你的情况下类似于 {a,b,c,d,e}） Iterator是一个迭代器 class，允许在 DIM 维数组上进行基于范围或基于迭代器的 for 循环。请注意operator[](size_t n)return 是 DIM 维坐标的第 n 个分量和advance()函数 return 是坐标为 [=28 的元素的迭代器=]Array应该很直接calcMarginals是我用来计算边缘的函数。dir这里是我要计算边际的方向（记住：一维边际）写到p_xPos和p_yPos，_p_theta是一个Array.注意基于迭代器的 for 循环，(*it)在这里指的是存储在数组中的 pdf 的双精度值，就像通常的迭代器一样。另外，it[dim]returns是dim方向实际值的坐标。写入 p_xPos 的最后一个循环只是因为我不想要这个数组中的索引而是实际值。
我想如果你重新定义 Iterator::operator[] 以获取维度索引的 vector/array 和 return 适当坐标的 vector/array 并添加 Array::operator[] 对于需要 vector/array 的随机访问，你应该已经完成了。
我解决了这个问题。我想我不能创建任意维度的 boost::multi_array，因为它需要维度作为模板参数，这需要在编译器时知道。这意味着我无法创建任意维度的视图。
因此，我做了以下事情： 我对 PDF 进行了排序，以便将要整合的维度是最后一个维度（很可能不是最有效的方法）。 然后我一个一个地缩小了PDF的尺寸。每个循环我只集成了 1 个维度，我将其保存在一个 multi_array 中，它与初始数组具有相同的大小（因为我无法使维度动态化）。 之后，我将值复制到缩小尺寸（已知）的 multi_array。
我使用以下 link 在维度上独立标注循环：
// Dimension-independent loop over boost::multi_array?