Minor

CMakeLists.txt

@@ -48,7 +48,8 @@ endif()
 message(STATUS "CMake Build Type: ${CMAKE_BUILD_TYPE}")
 
 # Set RPATH options
-set(CMAKE_INSTALL_RPATH "${CMAKE_BINARY_DIR}/lib")
+set(CMAKE_BUILD_RPATH "@loader_path/../lib:/usr/lib::/usr/local/lib")
+set(CMAKE_INSTALL_RPATH "@loader_path/../lib:/usr/lib::/usr/local/lib")
 
 # ----------------- #
 # Find Dependencies #

cmake/BuildFFTW3.cmake

@@ -21,6 +21,7 @@ set(FFTW_CONFIG_FLAGS
   --enable-threads
   --enable-sse2
   --enable-openmp
+  LDFLAGS="-L/path/to/libomp -lomp"
 )
 
 # Download and build FFTW3
@@ -39,4 +40,21 @@ ExternalProject_Add(fftw3_build
 # Set FFTW3 variables
 set(FFTW3_LIBRARY_DIRS ${CMAKE_BINARY_DIR}/lib/ CACHE INTERNAL "")
 set(FFTW3_INCLUDE_DIRS ${CMAKE_BINARY_DIR}/include/ CACHE INTERNAL "")
-set(FFTW3_LIBRARIES -lfftw3 -lfftw3f CACHE INTERNAL "")
+
+# Define FFTW3::Float
+add_library(FFTW3::Float SHARED IMPORTED GLOBAL)
+set_target_properties(FFTW3::Float PROPERTIES
+  IMPORTED_LOCATION "${CMAKE_BINARY_DIR}/lib/libfftw3f.dylib"
+  INTERFACE_INCLUDE_DIRECTORIES "${CMAKE_BINARY_DIR}/include"
+  INTERFACE_LINK_LIBRARIES OpenMP::OpenMP_CXX
+)
+
+# Define FFTW3::Double
+add_library(FFTW3::Double SHARED IMPORTED GLOBAL)
+set_target_properties(FFTW3::Double PROPERTIES
+  IMPORTED_LOCATION "${CMAKE_BINARY_DIR}/lib/libfftw3.dylib"
+  INTERFACE_INCLUDE_DIRECTORIES "${CMAKE_BINARY_DIR}/include"
+  INTERFACE_LINK_LIBRARIES OpenMP::OpenMP_CXX
+)
+
+

cmake/BuildPybind.cmake

@@ -35,24 +35,27 @@ function(build_pybind_target target libs)
   # Extract the target name
   get_filename_component(tarname ${target} NAME_WLE)
 
-  # Add the target library
-  add_library(${tarname} SHARED ${target})
+  if(NOT TARGET ${tarname})
 
-  # Link libraries to target
-  target_link_libraries(${tarname} "${libs}" OpenMP::OpenMP_CXX pybind11::headers)
+     # Add the target library
+     add_library(${tarname} SHARED ${target})
 
-  # Set system dependend properties
-  if(APPLE)
-    set_target_properties(${tarname} PROPERTIES LINK_FLAGS "-undefined dynamic_lookup")
-  else(APPLE)
-    target_link_libraries(${tarname} ${Python_LIBRARIES})
-  endif(APPLE)
+     # Link libraries to target
+     target_link_libraries(${tarname} "${libs}" OpenMP::OpenMP_CXX pybind11::headers)
 
-  # Set target prefix and suffix properties
-  pybind11_extension(${tarname})
-  pybind11_strip(${tarname})
+     # Set system dependend properties
+     if(APPLE)
+       set_target_properties(${tarname} PROPERTIES LINK_FLAGS "-undefined dynamic_lookup")
+     else(APPLE)
+       target_link_libraries(${tarname} ${Python_LIBRARIES})
+     endif(APPLE)
 
-  # Set the install path for bindings
-  install(TARGETS ${tarname} DESTINATION ${PYBIND_INSTALL_PATH})
+     # Set target prefix and suffix properties
+     pybind11_extension(${tarname})
+     pybind11_strip(${tarname})
+
+     # Set the install path for bindings
+     install(TARGETS ${tarname} DESTINATION ${PYBIND_INSTALL_PATH})
+    endif()
 
 endfunction()

src/CMakeLists.txt

@@ -250,13 +250,12 @@ if(BUILD_PYBIND)
   # Build bindings
   build_pybind("${core_libs}" FALSE)
   # Note: Uncomment the following line when pymrs issues have been resolved
-  # if(ONLY_SPARSE)
-  #   build_pybind("${core_libs}" FALSE)
-  # else()
-  #   set(pybind_deps ${mrs_libs} ${HEALPIX_LIBRARIES} ${core_libs})
-  #   build_pybind("${pybind_deps}" TRUE)
-  # endif()
-
+  if(ONLY_SPARSE)
+     build_pybind("${core_libs}" FALSE)
+  else()
+    set(pybind_deps ${mrs_libs} ${HEALPIX_LIBRARIES} ${core_libs})
+    build_pybind("${pybind_deps}" TRUE)
+  endif()
 endif()
 
 # ----------#

src/python_binding/alm.hpp

@@ -76,10 +76,8 @@ C_ALM::C_ALM(bool verbose)
     // (the -1 is used to live one thread to the main process which ensures
     // better and more constant performances). - Fabrice Poupon 2013/03/09
     #ifdef _OPENMP
-        if (nb_procs <= 0)
-            this->nb_procs = omp_get_num_procs() - 1;
-        else
-            this->nb_procs = nb_procs;
+        int np = omp_get_num_procs();
+        this->nb_procs = (np > 1) ? np - 1 : 1;
         omp_set_num_threads(this->nb_procs);
     #endif
 }
@@ -136,14 +134,13 @@ py::array_t<xcomplex<REAL>> C_ALM::get_tabalm()
 {
     int Nl=alm.Lmax()+1;
     int nelem= Nl*Nl;
-    auto arr = py::array_t<xcomplex<REAL>>(nelem);
+    auto arr = py::array_t<xcomplex<REAL>>({Nl, Nl});
     auto buffer = arr.request();
     xcomplex<REAL> *pointer = (xcomplex<REAL> *) buffer.ptr;
 
     for (int l=0; l <= alm.Lmax(); l++)
     for (int m=0; m <= l; m++)
         pointer[l + m * Nl] = alm(l,m);
-    arr.resize({Nl, Nl});
     return arr;
 }
 
@@ -164,6 +161,8 @@ void C_ALM::wiener(py::array_t<float>& psn, py::array_t<float>& pss)
     int Nl = alm.Lmax()+1;
     auto buffer = pss.request();
     float *pointer = (float *) buffer.ptr;
+    Cls.alloc(Nl);
+    Cln.alloc(Nl);
     for (int i=0; i<Nl;i++) Cls(i) = pointer[i];
     buffer = psn.request();
     pointer = (float *) buffer.ptr;
@@ -179,16 +178,14 @@ py::array_t<float> C_ALM::alm2spec()
     PowSpec ps_data;
     alm.alm2powspec(ps_data);
     int Nl = alm.Lmax()+1;
-    auto arr = py::array_t<xcomplex<REAL>>(Nl);
+    py::array_t<REAL> arr(Nl);
     auto buffer = arr.request();
-    xcomplex<REAL> *pointer = (xcomplex<REAL> *) buffer.ptr;
+    REAL* pointer = (REAL*) buffer.ptr;
     for (int l=0; l < Nl; l++)
         pointer[l] = ps_data.tt(l);
     return arr;
 }
 
-// py::array_t<float> alm2powspec();
-
 
 /*
 

src/python_binding/mrs.hpp

@@ -95,9 +95,9 @@ MRS::MRS(bool verbose)
 
 // Destructor
 MRS::~MRS(){
-    mr_initialized=False;
-    nb_procs=0;
-    Verbose=False;
+     mr_initialized=False;
+     nb_procs=0;
+     Verbose=False;
 }
 
 // int MRS::get_lmax(int  & Lmax, int Nside, float ZeroPadding)
@@ -122,7 +122,7 @@ void  MRS::alloc(int Nside, int Nscale, int LmaxIn, int ALM_iter, bool Verb)
     float ZeroPadding=0.;
     int LM=LmaxIn;
     int Lmax = mrs_get_lmax (LM,  Nside, ZeroPadding);
-    cout << "ALLOC: " << Nside << " " << Nscale << " " << LmaxIn << " " << ALM_iter << endl;
+    // cout << "ALLOC: " << Nside << " " << Nscale << " " << LmaxIn << " " << ALM_iter << endl;
     WT.wt_alloc(Nside, Nscale, Lmax, nested);  //  DEF_MRS_ORDERING);
     WT.ALM_iter = ALM_iter;
     if (Verb) Verbose=Verb;
@@ -133,14 +133,11 @@ void  MRS::alloc(int Nside, int Nscale, int LmaxIn, int ALM_iter, bool Verb)
 py::array_t<float>  MRS::uwt(py::array_t<float>& arr, int Ns)
 {
     Verbose=True;
-    cout << "Input arr.ndim() = " << arr.ndim() << " " << arr.shape(0) << endl;
+    // cout << "Input arr.ndim() = " << arr.ndim() << " " << arr.shape(0) << endl;
     if (arr.ndim() != 1)
         throw std::runtime_error("Input should be 1-D NumPy array");
 
-
-    cout << "OK" << endl;
-    exit(-1);
-
+    // cout << "OK" << endl;
 
     auto buffer = arr.request();
     float *pointer = (float *) buffer.ptr;
@@ -156,7 +153,7 @@ py::array_t<float>  MRS::uwt(py::array_t<float>& arr, int Ns)
         throw std::runtime_error("Input map has not expected nside.");
     }
 
-    cout << "alloc ok. Nside = " << Nside << " " << Ns << ", Npix = " << Npix << endl;
+    // cout << "alloc ok. Nside = " << Nside << " " << Ns << ", Npix = " << Npix << endl;
 
     // Copy the numpy array into a healpix map
     Hdmap Map;
@@ -168,27 +165,28 @@ py::array_t<float>  MRS::uwt(py::array_t<float>& arr, int Ns)
 
     // Wavelet transform
     WT.transform(Map);
-    cout << "transform ok. " << endl;
+    // cout << "transform ok. " << endl;
 
     // Recopy to numpy array
-    auto arr1 = py::array_t<float>(Npix*WT.nscale());
+    auto arr1 = py::array_t<float>({WT.nscale(), Npix});
     auto buf1 = arr1.request();
-    pointer = (float *) buf1.ptr;
+    pointer = static_cast<float*>(buf1.ptr);
+
     for(int s=0; s < WT.nscale(); s++)
     for (int i=0; i<Npix; i++)
     {
         pointer[i + s * Npix] = WT.WTTrans(i,s);
     }
-    cout << "copy ok. Npix =" << Npix << endl;
+    // cout << "copy ok. Npix =" << Npix << endl;
 
-    arr1.resize({WT.nscale(), Npix});
-    cout << "resize dim = " << arr1.ndim() << " " << arr1.shape(1)<< " " <<   arr1.shape(0) << endl;
     return arr1;
 
 }
 
 py::array_t<float> MRS::get_tabnorm()
 {
+    int n = WT.TabNorm.nx();
+    if (n <= 0) throw std::runtime_error("TabNorm has invalid size");
     auto arr1 = py::array_t<float>(WT.TabNorm.nx());
     auto buf1 = arr1.request();
     float *pointer = (float *) buf1.ptr;

src/sparse/libsparse2d/FCur.cc

@@ -1236,7 +1236,7 @@ void FCUR::get_wedges(Icomplex_f * &TabWT)
 
 void FCUR::cur_trans(Ifloat &Data)
 {
-	Verbose = True;
+	// Verbose = True;
 	if (Verbose == True) cout << "Transform WT ... " << endl;
 	if (ModifSize == False) transform(Data);
 	else