From e85b713cf31697a3309e12f30ba5759fee1cd3cc Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Fri, 20 Mar 2026 16:08:25 +0100 Subject: Trying out n linear interpolation --- lib/core/c++/math/n_linear.hpp | 44 +++++++++++++++++++++++++++++------------- 1 file changed, 31 insertions(+), 13 deletions(-) (limited to 'lib/core/c++/math') diff --git a/lib/core/c++/math/n_linear.hpp b/lib/core/c++/math/n_linear.hpp index 62906cd..4f47965 100644 --- a/lib/core/c++/math/n_linear.hpp +++ b/lib/core/c++/math/n_linear.hpp @@ -6,24 +6,42 @@ namespace kel { namespace lbm { template -void n_linear_interpolate(const saw::data& field, const saw::data>& pos){ +auto n_linear_interpolate(const saw::data& field, const saw::data>& pos){ + // Pos auto pos_bound = pos; + + // Dimensions auto meta = field.dims(); - saw::data ind; + + // Lower Index + saw::data> ind; + for(saw::data i{0u}; i < saw::data{D}; ++i){ - pos_bound.at(i).set( - std::max( - std::min( - meta.at(i).template cast_to().get(), - pos.at(i).get() + 1.5 - ), - 1, - ) - -1 - ); - ind.at(i) = pos_bound.at(i).template cast_to(); + // Native Positive i + auto npos_i = pos.at({i}).get(); + + { + // Ok I want to stay in bounds + npos_i = std::min(npos_i,meta.at(i).get()-1.0); + npos_i = std::max(npos_i,1.0); + } + + // Native Index i + auto nind_i = static_cast(std::floor(npos_i))-1ul; + + // Set index to i + ind.at(i).set(nind_i); } + saw::data> pos_frac; + for(saw::data i{0u}; i < saw::data{D}; ++i){ + pos_frac.at({i}) = pos_bound.at({i}) - ind.at(i).template cast_to(); + } + + // Base value + auto res = field.at(ind); + + constexpr uint64_t d_corners = 1ull << D; } } } -- cgit v1.2.3 From 16c1198f8fb401c3a98d927053fb2d29c2ce5f91 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Fri, 20 Mar 2026 18:17:55 +0100 Subject: Too much time spent on downcasts, but well. Here I am. --- lib/core/c++/math/n_linear.hpp | 76 ++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 73 insertions(+), 3 deletions(-) (limited to 'lib/core/c++/math') diff --git a/lib/core/c++/math/n_linear.hpp b/lib/core/c++/math/n_linear.hpp index 4f47965..bb7653c 100644 --- a/lib/core/c++/math/n_linear.hpp +++ b/lib/core/c++/math/n_linear.hpp @@ -1,12 +1,51 @@ #pragma once #include "../common.hpp" +#include "../iterator.hpp" namespace kel { namespace lbm { +namespace impl { +template +struct n_linear_interpolate_helper final { + template + auto apply(const saw::data& field, const saw::data>& pos){ + return pos; + } +}; +} + +template +saw::data,sch::Vector>> position_to_index_and_fraction(const saw::data>& pos){ + saw::data,sch::Vector>> sep; + + auto& ind = sep.template get<0u>(); + auto& frac = sep.template get<1u>(); + + auto pos_cpy = pos; + // Guarantee that the pos is at least 0 + for(uint64_t i = 0u; i < D; ++i){ + pos_cpy.at({{i}}).set(std::max(pos.at({{i}}).get(), static_cast::type>(0))); + } + + // Now we can cast to uint64_t + for(uint64_t i = 0u; i < D; ++i){ + ind.at({{i}}) = pos_cpy.at({{i}}).template cast_to(); + } + + frac = pos - ind.template cast_to(); + + return sep; +} + +template +auto floor_index_from_position(const saw::data>& pos){ + return position_to_index_and_fraction(pos).template get<0u>(); +} template -auto n_linear_interpolate(const saw::data& field, const saw::data>& pos){ +auto n_linear_interpolate( + const saw::data& field, const saw::data>& pos){ // Pos auto pos_bound = pos; @@ -39,9 +78,40 @@ auto n_linear_interpolate(const saw::data& field, const saw::data res; - constexpr uint64_t d_corners = 1ull << D; + // constexpr uint64_t d_corners = 1ul << D; + + saw::data> ones_ind; + for(saw::data i{0u}; i < saw::data{D}; ++i){ + ones_ind.at({i}).set(1u); + } + + iterator::apply([&](auto ind){ + // Iterates over (0,0,0) to (1,1,1) + saw::data weight{1.0}; + + for(saw::data d{0u}; d < saw::data{D}; ++d){ + + saw::data t = pos_frac.at({d}); + + if(ind.at(d).get() == 0u){ + weight = weight * (saw::data{1} - t); + }else{ + weight = weight * t; + } + } + + }, {}, ones_ind); +} + +template +saw::data> bilinear_interpolate(const saw::data& field, const saw::data>& pos){ + saw::data> res; + + { + } + return {}; } } } -- cgit v1.2.3 From 2a5261cd719547686f283b2aeae004b589b69e32 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Fri, 20 Mar 2026 18:23:39 +0100 Subject: Fraction and downcast looks ok. Need upper bound now --- lib/core/c++/math/n_linear.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'lib/core/c++/math') diff --git a/lib/core/c++/math/n_linear.hpp b/lib/core/c++/math/n_linear.hpp index bb7653c..ac02da8 100644 --- a/lib/core/c++/math/n_linear.hpp +++ b/lib/core/c++/math/n_linear.hpp @@ -33,7 +33,7 @@ saw::data,sch::Vector>> position_to_i ind.at({{i}}) = pos_cpy.at({{i}}).template cast_to(); } - frac = pos - ind.template cast_to(); + frac = pos_cpy - ind.template cast_to(); return sep; } -- cgit v1.2.3 From 14320632fbf39237b5f2254ec1d312ce3c12d879 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Fri, 20 Mar 2026 18:33:31 +0100 Subject: Wanted to finish n linearity today, but I guess I should be happy with this --- lib/core/c++/math/n_linear.hpp | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) (limited to 'lib/core/c++/math') diff --git a/lib/core/c++/math/n_linear.hpp b/lib/core/c++/math/n_linear.hpp index ac02da8..8fb0600 100644 --- a/lib/core/c++/math/n_linear.hpp +++ b/lib/core/c++/math/n_linear.hpp @@ -43,6 +43,25 @@ auto floor_index_from_position(const saw::data>& pos){ return position_to_index_and_fraction(pos).template get<0u>(); } +template +saw::data,sch::Vector>> position_to_index_and_fraction_bounded( + const saw::data>& pos, + const saw::data>& bound) +{ + auto infr = position_to_index_and_fraction(pos); + auto& ind = infr.template get<0u>(); + auto& fra = infr.template get<1u>(); + for(uint64_t i = 0u; i < D; ++i){ + // If index is higher than bound. Set to bound and reset fraction + if((ind.at({{i}}).get()+1u) >= bound.at({{i}}).get()){ + ind.at({{i}}).set(bound.at({{i}}).get()-1u); + fra.at({{i}}) = {}; + } + } + return infr; +} + + template auto n_linear_interpolate( const saw::data& field, const saw::data>& pos){ -- cgit v1.2.3