Commit ed8ec511 authored by ESTERIE Pierre's avatar ESTERIE Pierre

Some work on tensorial interpolation

parent c3250587
......@@ -126,12 +126,12 @@ if(MORSE_DISTRIB_DIR OR EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/modules/morse_cmake"
endif(NOT OpenMP_CXX_FOUND)
#
# XSIMD
# -----
if(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/modules/xsimd")
set(xsimd_DIR "${CMAKE_CURRENT_BINARY_DIR}/modules/xsimd")
# Find XSIMD properly
add_subdirectory(modules/xsimd)
# XTL, XSIMD & XTENSOR
# --------------------
#if(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/modules/xsimd")
# set(xsimd_DIR "${CMAKE_CURRENT_BINARY_DIR}/modules/xsimd")
# # Find XSIMD properly
# add_subdirectory(modules/xsimd)
find_package(xsimd CONFIG REQUIRED)
if(xsimd_FOUND)
message(STATUS "XSIMD found")
......@@ -139,7 +139,30 @@ if(MORSE_DISTRIB_DIR OR EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/modules/morse_cmake"
else(xsimd_FOUND)
message(FATAL_ERROR "Can't find XSIMD !")
endif(xsimd_FOUND)
endif(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/modules/xsimd")
#endif(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/modules/xsimd")
find_package(xtl CONFIG REQUIRED)
if(xtl_FOUND)
message(STATUS "XTL found")
target_link_libraries(${CMAKE_PROJECT_NAME} INTERFACE xtl)
else(xtl_FOUND)
message(FATAL_ERROR "Can't find XTL !")
endif(xtl_FOUND)
find_package(xtensor CONFIG REQUIRED)
if(xtensor_FOUND)
message(STATUS "XTENSOR found")
target_link_libraries(${CMAKE_PROJECT_NAME} INTERFACE xtensor)
else(xtensor_FOUND)
message(FATAL_ERROR "Can't find XTENSOR !")
endif(xtensor_FOUND)
#set(xtensor-blas_DIR $ENV{xtensor_blas_DIR})
#find_package(xtensor-blas CONFIG REQUIRED)
#if(xtensor-blas_FOUND)
# message(STATUS "XTENSOR-BLAS found")
# target_link_libraries(${CMAKE_PROJECT_NAME} INTERFACE xtensor-blas)
#else(xtensor-blas_FOUND)
# message(FATAL_ERROR "Can't find XTENSOR-BLAS !")
#endif(xtensor-blas_FOUND)
#
# Module inria
......@@ -157,7 +180,6 @@ if(MORSE_DISTRIB_DIR OR EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/modules/morse_cmake"
$<INSTALL_INTERFACE:include>
)
#
# Specific Debug flags
# --------------------
......@@ -517,7 +539,7 @@ if(MORSE_DISTRIB_DIR OR EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/modules/morse_cmake"
#
# Export Library
# --------------
install(TARGETS ${CMAKE_PROJECT_NAME} xsimd
install(TARGETS ${CMAKE_PROJECT_NAME} #xsimd
EXPORT ${CMAKE_PROJECT_NAME}-targets
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib
......
......@@ -4,6 +4,8 @@
# List of source files
set(source_tests_files
sandbox.cpp
test-blas.cpp
test-xtensor.cpp
)
# Add execs - 1 cpp = 1 exec
......@@ -27,6 +29,7 @@ foreach(exec ${source_tests_files})
set_target_properties(${execname} PROPERTIES ENABLE_EXPORTS TRUE
RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BUILD_TYPE}
)
target_compile_definitions(${execname} PUBLIC -DXTENSOR_USE_XSIMD)
target_link_libraries( ${execname} ${CMAKE_PROJECT_NAME})
# TODO: Remove these include derectories, ust temporary to reuse old components
target_include_directories( ${execname} PRIVATE ${CMAKE_SOURCE_DIR}/../include/)
......
//---------------------
// Experimental example
//---------------------
#include <iosfwd>
#include <string>
//#include <iosfwd>
//#include <string>
#include <vector>
#include <iostream>
#include <tuple>
#include <random>
#include <algorithm>
#include <inria/tcli/help_descriptor.hpp>
#include <inria/tcli/tcli.hpp>
......@@ -22,14 +23,19 @@
#include <range/v3/functional/concepts.hpp>
#include <range/v3/range_fwd.hpp>
#include <xsimd/xsimd.hpp>
#include <xtensor/xtensor.hpp>
#include <scalfmm/container/particle_container.hpp>
#include <scalfmm/container/point.hpp>
#include <scalfmm/container/particle.hpp>
#include <scalfmm/interpolation/uniform.hpp>
#include <scalfmm/tree/tree.hpp>
#include <scalfmm/tree/box.hpp>
#include <scalfmm/functional/utils.hpp>
#include <scalfmm/tools/colorized.hpp>
namespace args
{
struct nb_particle : inria::tcli::required_tag
......@@ -54,6 +60,13 @@ namespace args
type def{1};
};
struct order : inria::tcli::required_tag
{
inria::tcli::str_vec flags = {"--order", "-o"};
const char* description = "Precision order.";
using type = std::size_t;
};
struct file : inria::tcli::required_tag
{
inria::tcli::str_vec flags = {"--input-file", "-fin"};
......@@ -73,7 +86,8 @@ namespace args
nb_particle{},
tree_height{}, thread_count{},
file{}, inria::tcli::help{},
chunk_size{}
chunk_size{},
order{}
);
} // namespace args
......@@ -86,10 +100,12 @@ struct params {
std::size_t m_nb_threads{}; ///< Maximum thread count (when used).
std::string m_filename{}; ///< Particles file.
std::size_t m_group_size{}; ///< Group tree group size
std::size_t m_order{}; ///< Group tree group size
params() = default;
params(const decltype(args::cli)& cli)
{
m_order = cli.get<args::order>();
m_nb_particle = cli.get<args::nb_particle>();
m_tree_height = cli.get<args::tree_height>();
m_nb_threads = cli.get<args::thread_count>();
......@@ -101,70 +117,97 @@ struct params {
using namespace scalfmm;
//CPP_template(class Iter, class Sent, class Fun)
// (requires /* ranges::Invocable<Fun&> && ranges::OutputIterator<Iter,ranges::invoke_result_t<Fun&>> && ranges::Sentinel<Sent, Iter>*/
// ranges::Writable<Iter,ranges::invoke_result_t<Fun&>> )
//void my_algorithm(Iter first, Sent last, Fun f)
//{
// // ...
//}
CPP_template(class Iter, class Sent, class Fun)
(requires //ranges::Invocable<Fun&> && ranges::OutputIterator<Iter,ranges::invoke_result_t<Fun&>> && ranges::Sentinel<Sent, Iter>
ranges::Writable<Iter,ranges::invoke_result_t<Fun&>> )
void my_algorithm(Iter first, Sent last, Fun f)
{
// ...
}
// Quick type displayer
template<typename T>
class TD;
int main(int argc, char** argv)
{
args::cli.parse(argc,argv);
params parameters(args::cli);
const std::size_t tree_height{parameters.m_tree_height};
interpolation::uniform_interpolator<double,4> cheb_interpolator(parameters.m_order ,parameters.m_tree_height);
using source_t = container::particle<double,3,double,std::size_t>;
using position_t = container::point<double,3>;
using container_source_t = container::particle_container<source_t>;
const std::size_t nb_of_part{parameters.m_nb_particle};
//variad<int,float> ot{{8,9},{1.3,2.6}};
container_source_t cs{nb_of_part};
//variad<int, float> copy{ot};
//std::cout << colors::blue << std::get<0>(copy)[0] << " " << colors::green << std::get<1>(copy)[0] << colors::reset << '\n';
std::random_device rd;
std::mt19937 gen(rd());
std::uniform_real_distribution<double> dis(-1.0, 1.0);
auto random_r = [&dis, &gen](){ return dis(gen); };
//variad<int, float, double> v_{{1,2}, {3.3,4.4}, {5.5,6.6}};
//variad<int, float, double> v_(10);
//v_.push_back(std::make_tuple(2,1.1,2.2));
//std::cout << std::get<0>(v_).size() << '\n';
//std::cout << colors::blue << std::get<0>(v_)[0] << " " << colors::green << std::get<1>(v_)[0] << " " << colors::yellow << std::get<2>(v_)[0] << colors::reset << '\n';
auto make_particle = [&tree_height, &random_r]()
{
position_t pos = {random_r(), random_r(), random_r()};
source_t part( pos
, random_r()
, index::get_morton_index(pos, component::box<position_t>{2.0, {0.0,0.0,0.0}}, tree_height)
);
//const std::size_t tree_height{parameters.m_tree_height};
//using source_t = container::particle<double,3,double,std::size_t>;
//using position_t = container::point<double,3>;
//using container_source_t = container::particle_container<source_t>;
//const std::size_t nb_of_part{parameters.m_nb_particle};
//container_source_t cs{nb_of_part};
//std::random_device rd;
//std::mt19937 gen(rd());
//std::uniform_real_distribution<double> dis(-1.0, 1.0);
//auto random_r = [&dis, &gen](){ return dis(gen); };
//auto make_particle = [&tree_height, &random_r]()
//{
// position_t pos = {random_r(), random_r(), random_r()};
// source_t part( pos
// , random_r()
// , index::get_morton_index(pos, component::box<position_t>{2.0, {0.0,0.0,0.0}}, tree_height)
// );
// return part.as_tuple();
//};
//auto print_particle = [](source_t p)
//{
// std::cout << colors::yellow << "=================" << colors::reset << '\n';
// std::cout << "x = " << std::get<0>(p) << '\n';
// std::cout << "y = " << std::get<1>(p) << '\n';
// std::cout << "z = " << std::get<2>(p) << '\n';
// std::cout << "p = " << std::get<3>(p) << '\n';
// std::cout << "x.p = " << p.position()[0] << '\n';
// std::cout << "norm2 = " << utils::norm2(p.position()) << '\n';
// std::cout << "morton = " << std::get<4>(p) << '\n';
//};
//auto pred = [](auto a, auto b)
//{
// return std::get<4>(a) < std::get<4>(b);
//};
//std::vector<std::tuple<int,float>> v = {{8,3.0},{9,2.1},{10,9.6}};
//my_algorithm(cs.begin(), cs.end(), make_particle);
//std::generate(cs.begin(), cs.end(), make_particle);
//std::sort(cs.begin(), cs.end(), pred);
//std::for_each(cs.begin(), cs.end(), print_particle);
return part.as_tuple();
};
auto print_particle = [](source_t p)
{
std::cout << colors::yellow << "=================" << colors::reset << '\n';
std::cout << "x = " << std::get<0>(p) << '\n';
std::cout << "y = " << std::get<1>(p) << '\n';
std::cout << "z = " << std::get<2>(p) << '\n';
std::cout << "p = " << std::get<3>(p) << '\n';
std::cout << "x.p = " << p.position()[0] << '\n';
std::cout << "norm2 = " << utils::norm2(p.position()) << '\n';
std::cout << "morton = " << std::get<4>(p) << '\n';
};
auto pred = [](auto a, auto b)
{
return std::get<4>(a) < std::get<4>(b);
};
ranges::generate(cs.begin(), cs.end(), make_particle);
ranges::sort(cs.begin(), cs.end(), pred);
ranges::for_each(cs.begin(), cs.end(), print_particle);
//component::tree t{};
return 0;
}
......
// --------------------------------
// See LICENCE file at project root
// File : simple_variadic.cpp
// --------------------------------
#ifndef SCALFMM_CONTAINER_SIMPLE_VARIADIC_CONTAINER_HPP
#define SCALFMM_CONTAINER_SIMPLE_VARIADIC_CONTAINER_HPP
#include <vector>
#include <tuple>
#include <type_traits>
#include <utility>
#include <xsimd/config/xsimd_config.hpp>
namespace scalfmm::container
{
template<typename Types, typename Indices>
class variad_impl;
//template<typename Allocator = memory::aligned_allocator<XSIMD_DEFAULT_ALIGNMENT, Ts>>
template<typename... Ts, std::size_t... Indices>
class variad_impl<std::tuple<Ts...>, std::integer_sequence<std::size_t, Indices...>> : public std::tuple<std::vector<Ts, XSIMD_DEFAULT_ALLOCATOR(Ts)>...>
{
private:
//Discard fold expression results
struct noop_t
{
template<typename... Types>
noop_t(const Types&... ts) {}
};
public:
using base_type = std::tuple< std::vector<Ts, XSIMD_DEFAULT_ALLOCATOR(Ts)>... >;
using value_type = std::tuple< Ts... >;
using allocator_type = std::tuple<XSIMD_DEFAULT_ALLOCATOR(Ts)...>;
using size_type = std::size_t;
using difference_type = std::ptrdiff_t;
using reference = std::tuple<Ts&...>&;
using const_reference = const std::tuple<const Ts&...>&;
using pointer = std::tuple<Ts*...>;
using const_pointer = const std::tuple<const Ts*...>;
private:
allocator_type m_allocator{};
public:
variad_impl() = default;
variad_impl(const variad_impl&) = default;
variad_impl(variad_impl&&) noexcept = default;
variad_impl& operator=(const variad_impl&) = default;
variad_impl& operator=(variad_impl&&) noexcept = default;
~variad_impl() = default;
// TODO
//explicit variad( const allocator_type& alloc){}
explicit variad_impl( size_type count
, const value_type& value
, const allocator_type& alloc = allocator_type() )
: m_allocator(alloc)
{
reserve(count);
for(std::size_t i=0; i<count; ++i)
{
push_back(value);
}
}
explicit variad_impl(size_type count)
{
reserve(count);
for(std::size_t i=0; i<count; ++i)
{
push_back(value_type{});
}
}
inline void reserve(size_type size)
{
noop_t{(std::get<Indices>(*this).reserve(size),0) ...};
}
inline void push_back(const value_type& value)
{
noop_t{(std::get<Indices>(*this).push_back(std::get<Indices>(value)),0) ...};
}
};
template<typename... Ts>
struct variad : public variad_impl<std::tuple<Ts...>, std::make_index_sequence<sizeof...(Ts)>>
{
using base_type = variad_impl<std::tuple<Ts...>, std::make_index_sequence<sizeof...(Ts)>>;
using base_type::base_type;
};
}
#endif // SCALFMM_CONTAINER_SIMPLE_VARIADIC_CONTAINER_HPP
......@@ -1297,6 +1297,21 @@ namespace scalfmm::container
static_assert(sizeof...(Types) >= 1, "The vector must be instanciated with one or more types.");
};
template<typename... Types>
struct tuple_wrapper
{
std::tuple<Types...> data;
std::tuple<Types...>& get() const noexcept {return data;}
};
template<typename T>
T& stay(T&& t) {return t;}
template<typename... Types>
void swap(std::tuple<Types&...> a, std::tuple<Types&...> b)
{
//std::get<>
}
template<typename... Types, std::size_t... Indices>
class variadic_container_iterator<std::tuple<Types...>, inria::index_sequence<Indices...> >
......@@ -1320,9 +1335,13 @@ class variadic_container_iterator<std::tuple<Types...>, inria::index_sequence<In
#else
using typename base_t::tuple;
#endif
variadic_container_iterator(base_t tup_in) : base_t(tup_in){}
variadic_container_iterator(Types*... val) : base_t(val...){}
variadic_container_iterator(base_t tup_in, base_seq t) : base_t(tup_in){}
variadic_container_iterator(variadic_container_iterator const& other) = default;
variadic_container_iterator() : base_t(){}
/*
*/
// variadic_container_iterator(base_t tup_in):(tup_in){}
......@@ -1457,7 +1476,7 @@ class variadic_container_iterator<std::tuple<Types...>, inria::index_sequence<In
*
* \return A tuple of references to the values of the pointer element
*/
auto operator*() {
reference operator*() {
return reference(*std::get<Indices>(*this)...);
}
......
// See LICENCE file at project root
//
#ifndef SCALFMM_FUNCTIONAL_UTILS_HPP
#define SCALFMM_FUNCTIONAL_UTILC_HPP
#define SCALFMM_FUNCTIONAL_UTILS_HPP
#include <cmath>
......
// --------------------------------
// See LICENCE file at project root
// File : interpolation/builders.hpp
// --------------------------------
#ifndef SCALFMM_INTERPOLATION_BUILDERS_HPP
#define SCALFMM_INTERPOLATION_BUILDERS_HPP
#include <scalfmm/meta/const_functions.hpp>
#include <xtensor/xtensor.hpp>
#include <xtensor-blas/xlinalg.hpp>
namespace scalfmm::interpolation
{
template<std::size_t, typename T>
constexpr auto id(T t){ return std::forward<T>(t); };
template<typename Gen, std::size_t... I>
constexpr auto get_generator(Gen&& gen, std::index_sequence<I...>)
{
return std::move(xt::stack(std::make_tuple(id<I>(gen)...)));
};
template<std::size_t dim>
constexpr auto get(const std::size_t order)
{
return get_generator(xt::linspace(std::size_t{0},order-1, order), std::make_index_sequence<dim>{});
}
// TODO for nodes
//auto lin = xt::stack(xt::xtuple(xt::linspace(std::size_t{0},order-1,order)));
//xt::xarray<std::size_t> saved;
//for(std::size_t i=0; i<dimension; ++i)
//{
// auto r = xt::stack(xt::xtuple(xt::ones<std::size_t>({meta::pow(order,(dimension-1-i))})));
// auto tmp = xt::linalg::kron(lin,r);
// auto l = xt::stack(xt::xtuple(xt::ones<std::size_t>({meta::pow(order,i)})));
// auto reshaped = xt::linalg::kron(l,tmp);
// reshaped.reshape({1,-1});
// auto res = xt::stack(xt::xtuple(reshaped,saved));
// saved = res;
//}
// std::cout << "ids: " << saved << '\n';
template<std::size_t dim> struct nodes{ static_assert(dim<4, "Dimension for interpolation node not supported."); };
template<>
struct nodes<1>
{
inline static auto get(std::size_t order)
{
return xt::linspace(std::size_t{0}, order-1, order);
}
};
template<>
struct nodes<2>
{
inline static auto get(std::size_t order)
{
const auto s{meta::pow(order,2)};
const xt::xarray<std::size_t>::shape_type shape{2,s};
xt::xarray<std::size_t> ns(shape);
for(std::size_t i=0; i<s; ++i)
{
ns(0,i) = i%order;
ns(1,i) = (i/order)%order;
}
return ns;
}
};
template<>
struct nodes<3>
{
inline static auto get(std::size_t order)
{
const auto s{meta::pow(order,3)};
const xt::xarray<std::size_t>::shape_type shape{3,s};
xt::xarray<std::size_t> ns(shape);
for(std::size_t i=0; i<s; ++i)
{
ns(0,i) = i%order;
ns(1,i) = (i/order)%order;
ns(2,i) = i/(order*order);
}
return ns;
}
};
template<>
struct nodes<4>
{
inline static auto get(std::size_t order)
{
const auto s{meta::pow(order,4)};
const xt::xarray<std::size_t>::shape_type shape{4,s};
xt::xarray<std::size_t> ns(shape);
for(std::size_t i=0; i<s; ++i)
{
ns(0,i) = i%order;
ns(1,i) = (i/(order))%order;
ns(2,i) = (i/(order*order))%order;
ns(3,i) = i/(order*order*order);
}
return ns;
}
};
// Builders for relative center
}
#endif // SCALFMM_INTERPOLATION_BUILDERS_HPP
// --------------------------------
// See LICENCE file at project root
// File : interpolation/chebyshev_interpolator.hpp
// --------------------------------
#ifndef SCALFMM_INTERPOLATION_CHEBYSHEV_INTERPOLATOR_HPP
#define SCALFMM_INTERPOLATION_CHEBYSHEV_INTERPOLATOR_HPP
#include <array>
#include <cmath>
#include <type_traits>
#include <scalfmm/interpolation/chebyshev/chebyshev_root_factory.hpp>
#include <scalfmm/meta/const_functions.hpp>
#include <xtensor/xtensor.hpp>
#include <xtensor/xio.hpp>
namespace scalfmm::interpolation
{
template<std::size_t, typename T>
constexpr auto id(T t){ return std::forward<T>(t); };
template<typename Gen, std::size_t... I>
constexpr auto get_generator(Gen&& gen, std::index_sequence<I...>)
{
return std::move(xt::stack(std::make_tuple(id<I>(gen)...)));
};
template<std::size_t dim>
constexpr auto node_ids(const std::size_t order)
{
return get_generator(xt::linspace(std::size_t{0},order-1, order), std::make_index_sequence<dim>{});
}
template<typename ValueType, std::size_t Dimension>
struct chebyshev_interpolator
{
public:
using value_type = ValueType;
static constexpr std::size_t dimension = Dimension;
chebyshev_interpolator() = default;
chebyshev_interpolator(chebyshev_interpolator const &) = default;
chebyshev_interpolator(chebyshev_interpolator &&) = default;
chebyshev_interpolator& operator=(chebyshev_interpolator const&) = default;
chebyshev_interpolator& operator=(chebyshev_interpolator &&) = default;
explicit chebyshev_interpolator(const std::size_t order, std::size_t level)
: m_order(order), m_level(level), m_nnodes(meta::pow(order,dimension))
{
//constexpr unsigned int size_of_coef = std::pow(m_order, dimension);
//for(auto a : m_shape) { a = size_of_coef; }
//m_permutations.reshape(m_shape);
m_ts.resize({m_order,m_order});
//m_flat_ts.resize({(m_order)*m_order});
// Initializing roots
chebyshev_root_factory<value_type> cheb_fact{};
auto cheb_roots = cheb_fact.get(m_order);
// Calculating Ts
for( std::size_t o = 1; o < m_order; ++o )
for( std::size_t j = 0; j < m_order; ++j )
{ m_ts(o,j) = chebyshev::ts(o,cheb_roots(j)); }
auto lin = node_ids<dimension>(order);
std::cout << lin.size() << '\n';
std::cout << lin.dimension() << '\n';
std::cout << lin << '\n';
}
private:
const std::size_t m_order{};
const std::size_t m_level{};
const std::size_t m_nnodes{};
const std::array<std::size_t, dimension> m_shape{};
xt::xtensor<unsigned int, dimension> m_permutations{m_shape};
xt::xtensor<value_type, 2> m_ts{};
xt::xtensor<std::size_t, dimension> m_node_ids{};
//xt::xarray<ValueType,Dimension>
};
}
#endif // SCALFMM_INTERPOLATION_CHEBYSHEV_INTERPOLATOR_HPP
// --------------------------------
// See LICENCE file at project root
// File : interpolation/chebychev/chebyshev_root_factory.hpp
// --------------------------------
#ifndef SCALFMM_INTERPOLATION_CHEBYSHEV_CHEBYSHEV_ROOT_FACTORY_HPP
#define SCALFMM_INTERPOLATION_CHEBYSHEV_CHEBYSHEV_ROOT_FACTORY_HPP
#include <vector>
#include <exception>
#include <cmath>
#include <xtensor/xtensor.hpp>
#include <xtensor/xview.hpp>
namespace scalfmm::interpolation
{
namespace chebyshev
{
template<typename T>
constexpr T ts(const unsigned int n, T x)
{
assert(std::abs(x)-1.<10.*std::numeric_limits<T>::epsilon());
if (std::abs(x)>1.)
{
x = (x > T( 1.) ? T( 1.) : x);
x = (x < T(-1.) ? T(-1.) : x);
}
return T(cos(n * acos(x)));
}
template<typename T>
constexpr T us(const unsigned int n, T x)
{
assert(std::abs(x)-1.<10.*std::numeric_limits<T>::epsilon());
if (std::abs(x)>1.)
{
x = (x > T( 1.) ? T( 1.) : x);
x = (x < T(-1.) ? T(-1.) : x);
}
return T(n * (sin(n * acos(x))) / sqrt(1 - x*x));
}
} // chebyshev
template<typename ValueType, typename Enable = void>
struct chebyshev_root_factory
{
static_assert( std::is_floating_point<ValueType>::value