Merge branch 'transpose' into 'ducc_0_1'

add array transposition

See merge request !16

README.md

@@ -55,7 +55,7 @@ this is not a problem, as multiple DUCC versions can be installed
 simultaneously.
 The latest patch levels of a given DUCC version will always be available at the
 HEAD of the git branch with the respective name. In other words, if you need
-the latest incarnation of DUCC 0.1, this will be in branch "ducc_0_1" of the
+the latest incarnation of DUCC 0.1, this will be on branch "ducc_0_1" of the
 git repository, and it will be installed as the package "ducc_0_1".
 Later versions (like ducc_0_2 or ducc_1_0) will be maintained on new branches
 and will be installed as "ducc_0_2" and "ducc_1_0", so that there will be no
@@ -124,7 +124,7 @@ Library for high-accuracy 4pi convolution on the sphere, which generates a
 total convolution data cube from a set of sky and beam `a_lm` and computes
 interpolated values for a given list of detector pointings.
 This code has evolved from the original `totalconvolver` algorithm described
-in <https://arxiv.org/abs/astro-ph/0008227> vie the `conviqt` code
+in <https://arxiv.org/abs/astro-ph/0008227> via the `conviqt` code
 (<https://arxiv.org/abs/1002.1050>).
 
 

python/misc.cc

@@ -14,8 +14,6 @@
  *  You should have received a copy of the GNU General Public License
  *  along with DUCC; if not, write to the Free Software
  *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
- *
- *  For more information about HEALPix, see http://healpix.sourceforge.net
  */
 
 /*
@@ -32,6 +30,7 @@
 #include <vector>
 
 #include "mr_util/infra/mav.h"
+#include "mr_util/infra/transpose.h"
 #include "mr_util/math/fft.h"
 #include "mr_util/math/constants.h"
 #include "mr_util/math/gl_integrator.h"
@@ -155,6 +154,57 @@ py::array py_upsample_to_cc(const py::array &in, size_t nrings_out, bool has_np,
   return move(out);
   }
 
+template<typename T> py::array tphelp(const py::array &in)
+  {
+  auto in2 = to_fmav<T>(in, false);
+  auto out = make_Pyarr<T>(in2.shape());
+  auto out2 = to_fmav<T>(out, true);
+  transpose(in2, out2, [](const T &in, T &out){out=in;});
+  return out;
+  }
+
+py::array py_ascontiguousarray(const py::array &in)
+  {
+  if (isPyarr<float>(in))
+    return tphelp<float>(in);
+  if (isPyarr<double>(in))
+    return tphelp<double>(in);
+  if (isPyarr<complex<float>>(in))
+    return tphelp<complex<float>>(in);
+  if (isPyarr<complex<double>>(in))
+    return tphelp<complex<double>>(in);
+  if (isPyarr<int>(in))
+    return tphelp<int>(in);
+  if (isPyarr<long>(in))
+    return tphelp<long>(in);
+  MR_fail("unsupported datatype");
+  }
+
+template<typename T> py::array tphelp2(const py::array &in, py::array &out)
+  {
+  auto in2 = to_fmav<T>(in, false);
+  auto out2 = to_fmav<T>(out, true);
+  transpose(in2, out2, [](const T &in, T &out){out=in;});
+  return out;
+  }
+
+py::array py_transpose(const py::array &in, py::array &out)
+  {
+  if (isPyarr<float>(in))
+    return tphelp2<float>(in, out);
+  if (isPyarr<double>(in))
+    return tphelp2<double>(in, out);
+  if (isPyarr<complex<float>>(in))
+    return tphelp2<complex<float>>(in, out);
+  if (isPyarr<complex<double>>(in))
+    return tphelp2<complex<double>>(in, out);
+  if (isPyarr<int>(in))
+    return tphelp2<int>(in, out);
+  if (isPyarr<long>(in))
+    return tphelp2<long>(in, out);
+  MR_fail("unsupported datatype");
+  }
+
 
 const char *misc_DS = R"""(
 Various unsorted utilities
@@ -174,6 +224,9 @@ void add_misc(py::module &msup)
 
   m.def("upsample_to_cc",&py_upsample_to_cc, "in"_a, "nrings_out"_a,
     "has_np"_a, "has_sp"_a, "out"_a=py::none());
+
+  m.def("ascontiguousarray",&py_ascontiguousarray, "in"_a);
+  m.def("transpose",&py_transpose, "in"_a, "out"_a);
   }
 
 }

python/sht.cc

 /*
- *  This file is part of pyHealpix.
+ *  This file is part of DUCC.
  *
- *  pyHealpix is free software; you can redistribute it and/or modify
+ *  DUCC is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
  *  the Free Software Foundation; either version 2 of the License, or
  *  (at your option) any later version.
  *
- *  pyHealpix is distributed in the hope that it will be useful,
+ *  DUCC is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU General Public License for more details.
  *
  *  You should have received a copy of the GNU General Public License
- *  along with pyHealpix; if not, write to the Free Software
+ *  along with DUCC; if not, write to the Free Software
  *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
- *
- *  For more information about HEALPix, see http://healpix.sourceforge.net
  */
 
 /*
- *  pyHealpix is being developed at the Max-Planck-Institut fuer Astrophysik
+ *  DUCC is being developed at the Max-Planck-Institut fuer Astrophysik
  */
 
 /*

src/mr_util/infra/transpose.h

+/*
+ *  This file is part of DUCC.
+ *
+ *  DUCC is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  DUCC is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with DUCC; if not, write to the Free Software
+ *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+ */
+
+/*
+ *  DUCC is being developed at the Max-Planck-Institut fuer Astrophysik
+ */
+
+/*
+ *  Copyright (C) 2020 Max-Planck-Society
+ *  Author: Martin Reinecke
+ */
+
+#ifndef DUCC_TRANSPOSE_H
+#define DUCC_TRANSPOSE_H
+
+#include <algorithm>
+#include "mr_util/infra/error_handling.h"
+#include "mr_util/infra/mav.h"
+
+namespace mr {
+
+namespace detail_transpose {
+
+using namespace std;
+using shape_t=fmav_info::shape_t;
+using stride_t=fmav_info::stride_t;
+
+template<typename T1, typename T2> inline void rearrange(T1 &v, const T2 &idx)
+  {
+  T1 tmp(v);
+  for (size_t i=0; i<idx.size(); ++i)
+    v[i] = tmp[idx[i]];
+  }
+
+auto prep(const fmav_info &in, const fmav_info &out)
+  {
+  MR_assert(in.shape()==out.shape(), "shape mismatch");
+  shape_t shp;
+  stride_t si, so;
+  for (size_t i=0; i<in.ndim(); ++i)
+    if (in.shape(i)!=1) // remove axes of length 1
+      {
+      shp.push_back(in.shape(i));
+      si.push_back(in.stride(i));
+      so.push_back(out.stride(i));
+      }
+  // sort dimensions in order of descending output stride
+  vector<size_t> idx(shp.size());
+  iota(idx.begin(), idx.end(), 0);
+  sort(idx.begin(), idx.end(),
+    [&so](size_t i1, size_t i2) {return so[i1] > so[i2];});
+  rearrange(shp, idx);
+  rearrange(si, idx);
+  rearrange(so, idx);
+  // try merging dimensions
+//  [...]
+  if (shp.size()>1)
+    {
+    // move axis with smallest remaining input stride to second-to-last place
+    auto iminstr = min_element(si.begin(), si.end()-1) - si.begin();
+    size_t s2l = si.size()-2;
+    swap (shp[iminstr], shp[s2l]);
+    swap (si[iminstr], si[s2l]);
+    swap (so[iminstr], so[s2l]);
+    }
+
+  return make_tuple(shp, si, so);
+  }
+
+template<typename T, typename Func> void sthelper1(const T *in, T *out,
+  size_t s0, ptrdiff_t sti0, ptrdiff_t sto0, Func func)
+  {
+  for (size_t i=0; i<s0; ++i, in+=sti0, out+=sto0)
+    func(*in, *out);
+  }
+
+bool critical(ptrdiff_t s)
+  {
+  s = (s>=0) ? s : -s;
+  return (s>4096) && ((s&(s-1))==0);
+  }
+
+template<typename T, typename Func> void sthelper2(const T * MRUTIL_RESTRICT in,
+  T * MRUTIL_RESTRICT out, size_t s0, size_t s1, ptrdiff_t sti0, ptrdiff_t sti1,
+  ptrdiff_t sto0, ptrdiff_t sto1, Func func)
+  {
+  if ((sti0<=sti1) && (sto0<=sto1)) // no need to block
+    {
+    for (size_t i1=0; i1<s1; ++i1, in+=sti1, out+=sto1)
+      {
+      auto pi0=in;
+      auto po0=out;
+      for (size_t i0=0; i0<s0; ++i0, pi0+=sti0, po0+=sto0)
+        func(*pi0, *po0);
+      }
+    return;
+    }
+  if ((sti0>=sti1) && (sto0>=sto1)) // no need to block
+    {
+    for (size_t i0=0; i0<s0; ++i0, in+=sti0, out+=sto0)
+      {
+      auto pi1=in;
+      auto po1=out;
+      for (size_t i1=0; i1<s1; ++i1, pi1+=sti1, po1+=sto1)
+        func(*pi1, *po1);
+      }
+    return;
+    }
+  // OK, we have to do a real transpose
+  // select blockig sizes depending on critical strides
+  bool crit0 = critical(sizeof(T)*sti0) || critical(sizeof(T)*sto0);
+  bool crit1 = critical(sizeof(T)*sti1) || critical(sizeof(T)*sto1);
+  size_t bs0 = crit0 ? 8 : 8;
+  size_t bs1 = crit1 ? 8 : 8;
+  // make sure that the smallest absolute stride goes in the innermost loop
+  if (min(abs(sti0),abs(sto0))<min(abs(sti1),abs(sto1)))
+    {
+    swap(s0,s1);
+    swap(sti0,sti1);
+    swap(sto0,sto1);
+    swap(bs0, bs1);
+    }
+  for (size_t ii0=0; ii0<s0; ii0+=bs0)
+    {
+    size_t ii0e = min(s0, ii0+bs0);
+    for (size_t ii1=0; ii1<s1; ii1+=bs1)
+      {
+      size_t ii1e = min(s1, ii1+bs1);
+      for (size_t i0=ii0; i0<ii0e; ++i0)
+        for (size_t i1=ii1; i1<ii1e; ++i1)
+          func(in[i0*sti0+i1*sti1], out[i0*sto0+i1*sto1]);
+      }
+    }
+  }
+
+template<typename T, typename Func> void iter(const fmav<T> &in,
+  fmav<T> &out, size_t dim, ptrdiff_t idxin, ptrdiff_t idxout, Func func)
+  {
+  size_t ndim = in.ndim();
+  if (dim+2==ndim)
+    sthelper2(in.data()+idxin, out.vdata()+idxout, in.shape(ndim-2), in.shape(ndim-1),
+      in.stride(ndim-2), in.stride(ndim-1), out.stride(ndim-2),
+      out.stride(ndim-1), func);
+  else
+    for (size_t i=0; i<in.shape(dim); ++i)
+      iter(in, out, dim+1, idxin+i*in.stride(dim), idxout+i*out.stride(dim), func);
+  }
+
+template<typename T, typename Func> void transpose(const fmav<T> &in,
+  fmav<T> &out, Func func)
+  {
+  auto [shp, si, so] = prep(in, out);
+  fmav<T> in2(in, shp, si), out2(out, shp, so);
+  if (in2.ndim()==1)  // 1D, just iterate
+    {
+    sthelper1(in2.data(), out2.vdata(), in2.shape(0), in2.stride(0), out2.stride(0), func);
+    return;
+    }
+  iter(in2, out2, 0, 0, 0, func);
+  }
+
+}
+
+using detail_transpose::transpose;
+
+}
+
+#endif