Changed the way in which dimension reordering works.

src/combination/fubini.h

@@ -83,45 +83,26 @@ IntegratorFubini<IntegratorFirst,IntegratorRest,N> integrator_fubini(IntegratorF
 
 template<typename IntegratorFirst, typename IntegratorRest>
 auto integrator_fubini(const std::tuple<std::size_t>& dims, IntegratorFirst&& ifirst, IntegratorRest&& irest) {
- std::array<std::size_t,1> reorder{std::get<0>(dims)};
- return integrator_reorder(integrator_fubini<1>(std::forward<IntegratorFirst>(ifirst), std::forward<IntegratorRest>(irest)),reorder);
+ std::vector<std::tuple<size_t,size_t>> reorder{ {0, std::get<0>(dims)}};
+ return integrator_reorder(integrator_fubini<1>(std::forward<IntegratorFirst>(ifirst), std::forward<IntegratorRest>(irest)),{0,std::get<0>(dims)});
 } 
 
 template<typename IntegratorFirst, typename IntegratorRest>
 auto integrator_fubini(const std::tuple<std::size_t,std::size_t>& dims, IntegratorFirst&& ifirst, IntegratorRest&& irest) {
- std::array<std::size_t,2> reorder{std::get<0>(dims),std::get<1>(dims)};
+ std::vector<std::tuple<size_t,size_t>> reorder{ {0, std::get<0>(dims)}, {1, std::get<1>(dims)}};
  return integrator_reorder(integrator_fubini<2>(std::forward<IntegratorFirst>(ifirst), std::forward<IntegratorRest>(irest)),reorder);
 } 
 
 template<typename IntegratorFirst, typename IntegratorRest>
 auto integrator_fubini(const std::tuple<std::size_t,std::size_t,std::size_t>& dims, IntegratorFirst&& ifirst, IntegratorRest&& irest) {
- std::array<std::size_t,3> reorder{std::get<0>(dims),std::get<1>(dims),std::get<2>(dims)};
+ std::vector<std::tuple<size_t,size_t>> reorder{ {0, std::get<0>(dims)}, {1, std::get<1>(dims)}, {2, std::get<2>(dims)}};
  return integrator_reorder(integrator_fubini<3>(std::forward<IntegratorFirst>(ifirst), std::forward<IntegratorRest>(irest)),reorder);
 } 
 
 template<typename IntegratorFirst, typename IntegratorRest>
 auto integrator_fubini(const std::tuple<std::size_t,std::size_t,std::size_t,std::size_t>& dims, IntegratorFirst&& ifirst, IntegratorRest&& irest) {
- std::array<std::size_t,4> reorder{std::get<0>(dims),std::get<1>(dims),std::get<2>(dims),std::get<3>(dims)};
+ std::vector<std::tuple<size_t,size_t>> reorder{ {0, std::get<0>(dims)}, {1, std::get<1>(dims)}, {2, std::get<2>(dims)}, {3, std::get<3>(dims)}};
  return integrator_reorder(integrator_fubini<4>(std::forward<IntegratorFirst>(ifirst), std::forward<IntegratorRest>(irest)),reorder);
 } 
 
-template<typename IntegratorFirst, typename IntegratorRest>
-auto integrator_fubini(const std::tuple<std::size_t,std::size_t,std::size_t,std::size_t,std::size_t>& dims, IntegratorFirst&& ifirst, IntegratorRest&& irest) {
- std::array<std::size_t,5> reorder{std::get<0>(dims),std::get<1>(dims),std::get<2>(dims),std::get<3>(dims),std::get<4>(dims)};
- return integrator_reorder(integrator_fubini<5>(std::forward<IntegratorFirst>(ifirst), std::forward<IntegratorRest>(irest)),reorder);
-} 
-
-template<typename IntegratorFirst, typename IntegratorRest>
-auto integrator_fubini(const std::tuple<std::size_t,std::size_t,std::size_t,std::size_t,std::size_t,std::size_t>& dims, IntegratorFirst&& ifirst, IntegratorRest&& irest) {
- std::array<std::size_t,6> reorder{std::get<0>(dims),std::get<1>(dims),std::get<2>(dims),std::get<3>(dims),std::get<4>(dims),std::get<5>(dims)};
- return integrator_reorder(integrator_fubini<6>(std::forward<IntegratorFirst>(ifirst), std::forward<IntegratorRest>(irest)),reorder);
-} 
-
-template<typename IntegratorFirst, typename IntegratorRest>
-auto integrator_fubini(const std::tuple<std::size_t,std::size_t,std::size_t,std::size_t,std::size_t,std::size_t,std::size_t>& dims, IntegratorFirst&& ifirst, IntegratorRest&& irest) {
- std::array<std::size_t,7> reorder{std::get<0>(dims),std::get<1>(dims),std::get<2>(dims),std::get<3>(dims),std::get<4>(dims),std::get<5>(dims),std::get<6>(dims)};
- return integrator_reorder(integrator_fubini<7>(std::forward<IntegratorFirst>(ifirst), std::forward<IntegratorRest>(irest)),reorder);
-} 
-
-
 };

src/combination/reorder.h

@@ -11,37 +11,55 @@
 
 namespace viltrum {
 
-template<typename Integrator, std::size_t N>
+template<typename Integrator>
 class IntegratorReorder {
  Integrator integrator;
- std::array<std::size_t,N> first_dimensions;
+ std::vector<std::tuple<std::size_t,std::size_t>> to_from;
+
+ template<typename V>
+ std::vector<typename V::value_type> reorder(const V& v, const typename V::value_type& def = typename V::value_type()) const {
+
+ std::vector<bool> tos(v.size(),false), froms(v.size(),false);
+ for (auto [to, from] : to_from) {
+ if (from<v.size()) { 
+ if (to>=tos.size()) tos.resize(to+1,false);
+ froms[from]=true; tos[to]=true;
+ } 
+ }
+ std::size_t ntos = 0, nfroms = 0;
+ std::for_each(tos.begin(),tos.end(),[&ntos] (bool b) { if (b) ++ntos; });
+ std::for_each(froms.begin(),froms.end(),[&nfroms] (bool b) { if (b) ++nfroms; });
+
+ std::vector<typename V::value_type> sol(std::max(tos.size(),ntos + v.size()-nfroms),def);
+ for (auto [to, from] : to_from) {
+ if ((from<v.size()) && (to<sol.size())) {
+ sol[to]=v[from]; 
+ } 
+ }
+ std::size_t from, to = 0;
+ for (from = 0; from<v.size(); ++from) {
+ if (!froms[from]) { //We don't have this value yet
+ while ((to<sol.size()) && tos[to]) ++to; //Find where to put it
+ if (to<sol.size()) sol[to++] = v[from]; //Actually put it 
+ } 
+ } 
+ return sol; 
+ }
+
 public:
- IntegratorReorder(const Integrator& i, const std::array<std::size_t,N>& fd) :
- integrator(i), first_dimensions(fd) {}
+ IntegratorReorder(const Integrator& i, const std::vector<std::tuple<std::size_t,std::size_t>>& fd) :
+ integrator(i), to_from(fd) {}
 
  template<typename Bins, std::size_t DIMBINS, typename F, typename Float, std::size_t DIM, typename Logger>
  void integrate(Bins& bins, const std::array<std::size_t,DIMBINS>& bin_resolution,
  const F& f, const Range<Float,DIM>& range, Logger& logger) const {
 
- std::array<Float,DIM> reorder_min = range.min(); 
- std::array<Float,DIM> reorder_max = range.max();
- for (std::size_t i = 0; i<N; ++i) {
- if (first_dimensions[i]<DIM) {
- std::swap(reorder_min[i],reorder_min[first_dimensions[i]]);
- std::swap(reorder_max[i],reorder_max[first_dimensions[i]]);
- }
- }
- Range<Float,DIM> range_reordered(reorder_min,reorder_max);
-
-
+ auto range_reordered = viltrum::range<DIM>(reorder(range.min(),Float(0)),reorder(range.max(),Float(1)));
  const auto f_reordered = [&] (const std::array<Float,DIM>& x) {
- std::array<Float,DIM> x_reordered = x;
- for (std::size_t i = 0; i<N; ++i) {
- if (first_dimensions[i]<DIM) {
- std::swap(x_reordered[i],x_reordered[first_dimensions[i]]);
- }
- }
- return f(x_reordered);
+ auto re = reorder(x);
+ std::array<Float,DIM> x_re;
+ for (std::size_t i = 0; (i<DIM) && (i<re.size()); ++i) x_re[i] = re[i]; 
+ return (f(x_re));
  };
 
  integrator.integrate(bins,bin_resolution, f_reordered, range_reordered, logger);
@@ -51,38 +69,31 @@ class IntegratorReorder {
  void integrate(Bins& bins, const std::array<std::size_t,DIMBINS>& bin_resolution,
  const F& f, const RangeInfinite<Float>& range, Logger& logger) const {
 
- std::size_t size_reorder = 0; //For checking max size of the reorder vector
- for (std::size_t i : first_dimensions) if (size_reorder<=i) size_reorder = (i+1);
-
- std::vector<Float> reorder_min = range.min(); 
- if (reorder_min.size()<size_reorder) reorder_min.resize(size_reorder,0);
- std::vector<Float> reorder_max = range.max();
- if (reorder_max.size()<size_reorder) reorder_max.resize(size_reorder,1);
-
- for (std::size_t i = 0; i<N; ++i) {
- std::swap(reorder_min[i],reorder_min[first_dimensions[i]]);
- std::swap(reorder_max[i],reorder_max[first_dimensions[i]]);
- } 
- RangeInfinite<Float> range_reordered(reorder_min,reorder_max); 
-
+ RangeInfinite<Float> range_reordered(reorder(range.min(),Float(0)),reorder(range.max(),Float(1))); 
+ std::size_t size_reorder = 0;
+ for (auto [to, from] : to_from) {
+ if (size_reorder<=from) {
+ size_reorder = from + 1; 
+ } 
+ } 
  const auto f_reordered = [&] (const auto& x) {
- std::vector<Float> x_reordered(size_reorder);
+ std::vector<Float> xr(size_reorder);
  auto it = x.begin();
- for (Float& xr : x_reordered) { xr = (*it); ++it; }
-  
- for (std::size_t i = 0; i<N; ++i) {
-  std::swap(x_reordered[i],x_reordered[first_dimensions[i]]);
- }
- return f(concat(x_reordered,sequence_from_iterators(it,x.end())));
+ for (Float& v : xr) { v = (*it); ++it; }
+/* for (Float v : xr) std::cerr<<v<<" ";
+ std::cerr<<" -> ";
+ for (Float v : reorder(xr)) std::cerr<<v<<" ";
+ std::cerr<<std::endl; */
+ return f(concat(reorder(xr),sequence_from_iterators(it,x.end())));
  };
 
  integrator.integrate(bins,bin_resolution, f_reordered, range_reordered, logger);
  }
 };
 
-template<typename Integrator, std::size_t N>
-IntegratorReorder<Integrator,N> integrator_reorder(Integrator&& i, const std::array<std::size_t,N>& fd) {
- return IntegratorReorder<std::decay_t<Integrator>,N>(std::forward<Integrator>(i),fd);
+template<typename Integrator>
+IntegratorReorder<Integrator> integrator_reorder(Integrator&& i, const std::vector<std::tuple<std::size_t,std::size_t>>& fd) {
+ return IntegratorReorder<std::decay_t<Integrator>>(std::forward<Integrator>(i),fd);
 }
 
 };

src/range.h

@@ -161,6 +161,18 @@ Range<T,N> range_all(const T& va, const T& vb) {
  return range(a,b);
 }
 
+template<std::size_t N, typename T>
+Range<T,N> range(const std::vector<T>& a, const std::vector<T>& b) {
+ std::array<T,N> ar, br;
+ std::size_t i;
+ for (i = 0; (i<N) && (i<a.size()); ++i) ar[i] = a[i];
+ for (;i<N; ++i) ar[i]=0; 
+ for (i = 0; (i<N) && (i<b.size()); ++i) br[i] = b[i];
+ for (;i<N; ++i) br[i]=1;
+ return range(ar,br); 
+} 
+
+
 template<std::size_t N, typename T = float>
 Range<T,N> range_primary() {
  return range_all<N>(T(0),T(1));