初始化未知维数的数组

Question

我很惊讶居然找不到这个问题。我试着将它一般化(使用一些不错的未经测试的代码)，使之成为每个人都能从中受益的东西。

假设我有一个多维的Point:

template <int dims> class Point { public: double data[dims]; };

现在我创建一个多维数组:

 template <int dims> void foobar(int count0, ...) {
    //Using variadic function.  Could also use variadic templates in C++ (arguably better)
    int counts[dims], total_count=count0; counts[0]=count0;
    va_list args; va_start(args,count0);
    for (int i=1;i<dims;++i) {
        int count = va_arg(args,int);
        counts[i] = count;
        total_count *= count;
    }
    va_end(args);
    Point<dims>* array = new Point<dims>[total_count];
    //...
}

可以看到，array是未知维数的多维数组，用一维数组表示。

---

我的问题:我如何才能干净地初始化这个数组到它的多维网格点?

这是我想要在1、2和3维度上的示例行为。很明显，我不想为我想要使用的每一个可能的维度写这个!目标是概括这个。

//Example: dim==1
for (int x=0; x<counts[0]; ++x) {
    Point<1>& point = array[x];
    point.data[0] = (x+0.5) / (double)counts[0];
}
//Example: dim==2
for (int y=0; y<counts[1]; ++y) {
    for (int x=0; x<counts[0]; ++x) {
        Point<2>& point = array[y*counts[0]+x];
        point.data[0] = (x+0.5) / (double)counts[0];
        point.data[1] = (y+0.5) / (double)counts[1];
    }
}
//Example: dim==3
for (int z=0; z<counts[2]; ++z) {
    for (int y=0; y<counts[1]; ++y) {
        for (int x=0; x<counts[0]; ++x) {
            Point<3>& point = array[(z*counts[1]+y)*counts[0]+x];
            point.data[0] = (x+0.5) / (double)counts[0];
            point.data[1] = (y+0.5) / (double)counts[1];
            point.data[2] = (z+0.5) / (double)counts[2];
        }
    }
}

再一次，我的问题是:以一种干净的方式将上述内容推广到任意数量的嵌套循环/维度。

注意:我已经想出了一些讨厌的方法，它们既不优雅又缓慢。特别是，如果可能的话，我希望避免递归，因为它将经常在高维小型数据集上被调用。注意:在C中有明显的相似之处，所以C或c++都可以。c++ 11优先。

Answer 1

回复评论和更新的问题编辑

如果你需要性能和"优雅"，我会:

• 放弃多维数组方法并将其扁平化(即一个数组维度)。没有new，没有指针，使用c++现代方法与std::vector或std::array。
• 为您的多dim数组提供一个抽象容器，使用方便的方法，如通用嵌套循环"生成器"
• 将可变参数替换为固定大小的数组(因为您在编译时知道dims)。

所以我找到了一个以下的解决方案，这是非常符合您的实现和需求，并尽量保持简单。

我已经设法以"现代c++ 11方式"重写了一个小MultiArray类。我认为count维度在编译时可能不知道，因此现在使用std::vector。当然，可以使用std::array获得更通用的编译时代码，请参阅下面我的原始答案。

#include <iostream>
#include <array>
#include <vector>
#include <numeric>
template<size_t DIMS>
class MultiArray {
public:
    // Point here is just an array
    using Point = std::array<double,DIMS>;
    // fill data_ with an init array
    // not that count is just a fix sized array here no variadic arguments needed
    MultiArray(const std::array<size_t, DIMS>& count) 
        : data_{init_array(count)} {}
private:   
    // the init functions are used for the constructor
    void init_point(Point& point, const std::array<size_t,DIMS>& coord, const std::array<size_t, DIMS>& count) {
        std::cout << " -> { ";
        for (size_t i = 0; i < DIMS; i ++) {
            point[i] = (coord[i] + 0.5) / count[i];
            std::cout << point[i] << ";";
        }
        std::cout << " }n";
    }
    std::vector<Point> init_array(const std::array<size_t, DIMS>& count) {
        std::vector<Point> data(std::accumulate(count.begin(), count.end(), 1, std::multiplies<int>())); // accumulate computes the prod of DIMS total_count
        std::array<size_t, DIMS> current{};
        size_t i=0;
        do {
             for (size_t i = 0; i < DIMS; i ++)
                 std::cout << current[i] << ";";
            init_point(data[i++],current,count);
        } while (increment(current, count));
        return data;
    }
    // the following function allows to imitate the nested loop by incrementing multidim coordinates
    bool increment( std::array<size_t, DIMS>& v, const std::array<size_t, DIMS>& upper) {
        for (auto i = v.size(); i-- != 0; ) {
            ++v[i];
            if (v[i] != upper[i]) {
                return true;
            }
            v[i] = 0;
        }
        return false;
    }   
private:
    std::vector<Point> data_; // A flatten multi dim vector of points
};

int main() {
   std::array<size_t, 3> count{{4,5,3}};
   MultiArray<3> test{count};
}

Live on Coliru

在结果中可以看到，data_可以为N维度初始化。如果你需要一个更高层次的抽象类，你可以检查下面我的原始答案，你可以执行一些方便的操作(即访问grid[{i,j,k}]来填充值)。

原始回答

我需要一个多维网格来满足我的需求，并且碰巧在代码审查时要求改进我的代码。这里有一个工作的例子，当然有些功能对你来说可能是不必要的…我的实现涉及模板和编译时计算。注意，维度大小必须在编译时已知。

简单地说，这个类看起来像这样:

template<typename T, size_t... DIMS> // variadic template here for the dimensions size
class MultiGrid {
   // Access from regular idx such as grid[64]
   T& operator[] (size_type idx)       { return values_[idx]; };
   // Access from multi dimensional coordinates such as grid[{6,3,4}]
   T& operator[] (const std::array<size_t,sizeof...(DIMS)>& coord)       { // can give code for runtime here };
private:
  std::array<T,sizeof...(DIMS)> data_;
}

然后你可以构造你的多维数组并以这种方式初始化它:

MultiGrid<float,DIM1,DIM2,DIM3> data; // 3D
// MultiGrid<float,DIM1,DIM2,DIM3,DIM4> data; // 4D
// etc...
// initialize it like this with nested arrays
for (size_t z=0; z < DIM3; z ++)
  for (size_t y=0; y < DIM2; y ++)
    for (size_t x=0; x < DIM1; x ++)
      data[{x,y,z}] = [...] // whatever
// or like this in C++11/14 way
for (auto &x : data) x = [...] // this is convenient to provide a container like approach since no nested arrays are needed here.

如果你需要为可变嵌套循环指定一个算法来填充值，你可以看看这里，用第一个答案这样做:

// here lower_bound is 0-filled vector
std::vector<int> current = lower_bound;
do {
   data[current] = [...] // fill in where current is a coordinate
} while (increment(current, lower_bound, upper_bound));

如果您需要我在实现中遗漏的东西，请随时询问。如果有人能指出改进的地方，我也会很高兴。

Answer 2

从X,Y,Z到扁平数组(F)，我们有以下公式

F=(Z*DimY+y)*DimX+X

或

F=Z*DimY*DimX+Y*DimX+X
X = F % DimX
Y = F % DimX*DimY/DimX
Z = F % DimX*DimY*DimZ/DimX*DimY

在一个7 × 3 × 5的数组中，Z=3, Y=1, x =2将会在3 *3*5 + 1 *5 + 2 = 45+5+2= 52

X = `52` % 5 = 2
Y = `52` % (5 * 3) / 5 = 7 / 5 = 1
Z = `52` % (7 * 5 * 3)/15 = 52/15 = 3

在一个7 × 3 × 5的数组中，Z=4, Y=2, x =3的值为:4 *3*5 + 2 *5 + 3 = 60+10+3= 73

X = `73` % 5 = 3
Y = `73` % (5 * 3) / 5 = 13 / 5 = 2
Z = `73` % (7 * 5 * 3)/15 = 73/15 = 4

如果我们将累积积保存在一个数组中，mult { 1, X, X*Y, X*Y*Z, ...}和所有点保存在一个数组中，val

指向平面数组:

F=sum(mult[i]*val[i]);

指向点的平面数组:

i[0]=F%mult[1]/mult[0];
i[1]=F%mult[2]/mult[1];
...

然后我们可以遍历F(平面数组)，从索引反向工程到平面数组的所有点:X,Y，…如上所述，并在泛型循环中执行所需的初始化:

给定mult为mult[0]=1; mult[d+1]=mult[d]*count[d];

for (int i = 0; i < total_count; ++i) {
    for (int d=0; d < dims; ++d) {
      int dim=(i%mult[d+1])/mult[d];
      point.data[d] = (dim+0.5) / (double)counts[d];
  }
}

Answer 3

这是我的解决方案，使用c++ 11可变模板和包扩展。

我的"foobar"作为constexpr函数编译，所以我想说这是相当有效的:p

至于优雅，你可以评判，但我认为它很好。

基本上，这个想法是用函数式编程的迭代方式取代for循环，我们只是寻求显式地构建我们想要首先迭代的集合。然后，该代码可以以更直接的方式推广到任意维度。

代码是完全独立的，保存<array>头。

它为我在c++ 11标准下编译，gcc 4.8.2和clang 3.6。

注意:如果您想在代码中使用相同的技术，其中维度是一个运行时参数，基本上您要做的是使用std::vector<std::vector<size_t>>之类的东西将下面的Cartesian_Product元函数重新实现为运行时函数。然后，您可以通过获取dim个数的笛卡尔积来构建迭代集，并使用单个for循环对其进行迭代以填充结果。

#include <array>
#include <iostream>
/////////////////////////////////////////////
// The point class from problem statement
/////////////////////////////////////////////
template <size_t dims> class Point { public: double data[dims]; };
// Some basic type list and meta function objects to work with
// List of indices
template<size_t... sl>
struct StList {
    static constexpr size_t length = sizeof...(sl);
};
// Metafunction to compute a volume
template <typename T>
struct Compute_Volume;
template<size_t s>
struct Compute_Volume<StList<s>> {
    static constexpr size_t value = s;
};
template <size_t s1, size_t s2, size_t... sl>
struct Compute_Volume<StList<s1, s2, sl...>> {
    static constexpr size_t value = s1 * Compute_Volume<StList<s2, sl...>>::value;
};
// Concatenate
template<typename TL, typename TR>
struct Concat;
template<size_t... SL, size_t... SR>
struct Concat<StList<SL...>, StList<SR...>> {
    typedef StList<SL..., SR...> type;
};
template<typename TL, typename TR>
using Concat_t = typename Concat<TL, TR>::type;
// Meta function to check if a typelist is component-wise less than another typelist
// For testing purposes
template <typename TL, typename TR>
struct all_less_than;
template <size_t l1, size_t... SL, size_t r1, size_t... SR>
struct all_less_than<StList<l1, SL...>, StList<r1, SR...>> {
    static constexpr bool value = (l1 < r1) && all_less_than<StList<SL...>, StList<SR...>>::value;
};
template<>
struct all_less_than<StList<>, StList<>> {
    static constexpr bool value = true;
};
/////////////////////////////////////////////////////////////////////////////
// constexpr template function for the point initializations you described
/////////////////////////////////////////////////////////////////////////////
template <typename index, typename dims>
struct point_maker;
template <size_t... idx, size_t... dim>
struct point_maker<StList<idx...>, StList<dim...>> {
    static_assert(sizeof...(idx) == sizeof...(dim), "misuse of 'point_maker' template, idx and dim must have same number of coordinates");
    static_assert(all_less_than<StList<idx...>, StList<dim...>>::value, "misuse of 'point_maker' template, idx is out of bounds");
    static constexpr Point<sizeof...(idx)> make_point() {
        return {{ ((idx + 0.5) / static_cast<double>(dim))... }};
    }
};
//////////////////////////////////////////////////////////////////////////////////////////
// Now we need to abstract the for loops. We need a little more infrastructure for this.
//////////////////////////////////////////////////////////////////////////////////////////
// A basic typelist object
template <typename... tl>
struct TypeList {
    static constexpr size_t length = sizeof...(tl);
};
// Specialization for the Concat metafunction
template<typename... TL, typename... TR>
struct Concat<TypeList<TL...>, TypeList<TR...>> {
    typedef TypeList<TL..., TR...> type;
};

// Metafunction to compute the cartesian product of two lists of lists, and evaluate the `Concat` metafunction for each pair.
template <typename S, typename T>
struct Cartesian_Product;
template <typename s1, typename s2, typename... sl, typename T>
struct Cartesian_Product<TypeList<s1, s2, sl...>, T> {
    typedef Concat_t<typename Cartesian_Product<TypeList<s1>, T>::type, typename Cartesian_Product<TypeList<s2, sl...>, T>::type> type;
};
template <typename s, typename t1, typename t2, typename... tl>
struct Cartesian_Product<TypeList<s>, TypeList<t1, t2, tl...>> {
    typedef Concat_t<typename Cartesian_Product<TypeList<s>, TypeList<t1>>::type, typename Cartesian_Product<TypeList<s>, TypeList<t2, tl...>>::type> type;
};
template <typename s, typename t>
struct Cartesian_Product<TypeList<s>, TypeList<t>> {
    typedef TypeList<Concat_t<s, t>> type;
};
template <typename S, typename T>
using Cartesian_Product_t = typename Cartesian_Product<S, T>::type;
// Some unit tests for the above :)
static_assert( std::is_same<Cartesian_Product_t<TypeList<StList<1>, StList<2>>, TypeList<StList<3>, StList<4>>>, TypeList<StList<1,3>, StList<1,4>, StList<2,3>, StList<2,4>>>::value , "cartesian product not working");
static_assert( std::is_same<Cartesian_Product_t<TypeList<StList<1,5>, StList<2>>, TypeList<StList<3>, StList<4>>>, TypeList<StList<1,5,3>, StList<1,5,4>, StList<2,3>, StList<2,4>>>::value , "cartesian product not working");
static_assert( std::is_same<Cartesian_Product_t<TypeList<StList<1,5>, StList<2>>, TypeList<StList<>>>, TypeList<StList<1,5>, StList<2>>>::value , "cartesian product not working");
// Metafunction to count from 0 to n
// Count from zero to n-1, and make a sequence of singleton sets containing the numbers
template<size_t s>
struct Count {
    typedef Concat_t<typename Count<s-1>::type, TypeList<StList<s-1>>> type;
};
template<>
struct Count<0> {
    typedef TypeList<> type;
};
template<size_t s>
using Count_t = typename Count<s>::type;

// Metafunction to abstract a series of for loops, generating the list of all the index tuples the collection of loops would generate
template <typename T>
struct Volume_Maker;
template <>
struct Volume_Maker<StList<>> {
    typedef TypeList<StList<>> type;
};
template <size_t s, size_t... sl>
struct Volume_Maker<StList<s, sl...>> {
    typedef Cartesian_Product_t<Count_t<s>, typename Volume_Maker<StList<sl...>>::type> type;
};
template <typename T>
using Volume_t = typename Volume_Maker<T>::type;
// Some more quick unit tests
template <typename T>
struct Volume_Test {
    static_assert( Volume_t<T>::length == Compute_Volume<T>::value, "volume computation mismatch");
};
Volume_Test<StList<1,2,3>> test1{};
Volume_Test<StList<1,1,1>> test2{};
Volume_Test<StList<1>> test3{};
Volume_Test<StList<7,6,8>> test4{};
/////////////////
// Grand finale
/////////////////
template <typename dim_list, typename idx_list = Volume_t<dim_list>>
struct foobar_helper;
template <size_t... dim, typename... idx>
struct foobar_helper<StList<dim...>, TypeList<idx...>> {
    typedef StList<dim...> dim_list;
    typedef std::array<Point<sizeof...(dim)>, sizeof...(idx)> array_type;
    static constexpr array_type make_array() {
        return {{ point_maker<idx, dim_list>::make_point()... }};
    }
};
template <size_t... dim>
constexpr std::array<Point<sizeof...(dim)>, Compute_Volume<StList<dim...>>::value> foobar() {
    return foobar_helper<StList<dim...>>::make_array();
}
/////////////////
// Example usage
/////////////////
template <size_t ndim>
std::ostream & operator << (std::ostream & s, const Point<ndim> & p) {
    s << "[ ";
    for (size_t i = 0; i < ndim; ++i) {
        if (i) { s << ", "; }
        s << p.data[i];
    }
    s << " ]";
    return s;
}
template<size_t ndim, size_t N>
void print_array(std::ostream & s, const std::array<Point<ndim>, N> & a) {
    for (size_t i = 0; i < N; ++i) {
        s << a[i] << std::endl;
    }
}
int main() {
    constexpr auto array1 = foobar<2,2,2>();
    print_array(std::cout, array1);
    constexpr auto array2 = foobar<3,3,3,3>();
    print_array(std::cout, array2);
}