Add Regression support (incomplete)

src/regression/poly_levelset_test.cpp

@@ -52,15 +52,15 @@ BOOST_AUTO_TEST_CASE ( PolyLevelset_Sphere )
     }    
     vd.map();
     
-    auto model = PolyLevelset<3>(vd, 1e-4);
-
+    // auto model = PolyLevelset<3>(vd, 1e-4);
+/*
     double max_err = -1.0;
     auto it2 = vd.getDomainIterator();
     while (it2.isNext())
     {
         auto key = it2.get();
         Point<3, double> pos = {vd.getPos(key)[0], vd.getPos(key)[1], vd.getPos(key)[2]};
-        vd.template getProp<mean_curvature>(key) = model.estimate_mean_curvature_at(pos);
+        // vd.template getProp<mean_curvature>(key) = model.estimate_mean_curvature_at(pos);
         // vd.template getProp<gauss_curvature>(key) = model->estimate_gauss_curvature_at(vd.getPos(key));
 
         double val = vd.getProp<mean_curvature>(key);
@@ -79,8 +79,9 @@ BOOST_AUTO_TEST_CASE ( PolyLevelset_Sphere )
     else
         check = false;
     std::cout<<"Max err (poly level) = "<<max_err<<"\n";
+    
     BOOST_TEST( check );
-
+    */
     // if(model)
     //     delete model;
 

src/regression/regression.hpp

@@ -11,31 +11,41 @@
 #include "Vector/map_vector.hpp"
 #include "Space/Shape/Point.hpp"
 #include "DMatrix/EMatrix.hpp"
-
+#include "DCPSE/SupportBuilder.hpp"
 #include "minter.h"
 
+/*
+
 template<typename vector_type, typename NN_type>
-class RegressionDomain
+class RegressionSupport
 {
 	openfpm::vector_std<size_t> keys;
 	
 public:
 
-	RegressionDomain(vector_type &vd, Point<vector_type::dims, typename vector_type::stype> pos, typename vector_type::stype size, NN_type &NN)
+	template<typename iterator_type>
+	RegressionSupport(vector_type &vd, iterator_type itPoint, unsigned int requiredSize, support_options opt, NN_type &cellList)
 	{
-		// Not efficient to compute cell list each time.
-        //auto NN = vd.getCellList(size);
-        auto Np = NN.template getNNIterator(NN.getCell(pos));
-     	
-        while(Np.isNext())
-        {
-        	auto key = Np.get();
-        	if (pos.distance(vd.getPos(key)) < size)
-            	keys.add(key);
-            ++Np;
-        }
-	}
+		// Get spatial position from point iterator
+        vect_dist_key_dx p = itPoint.get();
+        vect_dist_key_dx pOrig = itPoint.getOrig();
+        Point<vector_type::dims, typename vector_type::stype> pos = vd.getPos(p.getKey());
+
+		// Get cell containing current point and add it to the set of cell keys
+        grid_key_dx<vector_type::dims> curCellKey = cellList.getCellGrid(
+                pos); // Here get the key of the cell where the current point is
+        std::set<grid_key_dx<vector_type::dims>> supportCells;
+        supportCells.insert(curCellKey);
 
+        // Make sure to consider a set of cells providing enough points for the support
+        enlargeSetOfCellsUntilSize(supportCells, requiredSize + 1,
+                                   opt, cellList); // NOTE: this +1 is because we then remove the point itself
+
+        // Now return all the points from the support into a vector
+
+        keys = getPointsInSetOfCells(supportCells, p, pOrig, requiredSize, opt);
+	}
+	
 	auto getKeys()
 	{
 		return keys;
@@ -46,6 +56,7 @@ public:
 		return keys.size();
 	}
 };
+*/
 
 template<int spatial_dim, unsigned int prp_id, typename MatType = EMatrixXd, typename VecType = EVectorXd>
 class RegressionModel
@@ -158,7 +169,7 @@ public:
 	    		if (err > error)
 	    			error = err;
 	    	}
-	    	std::cout<<"Fit of degree "<<poly_degree<<" with error = "<<error<<std::endl;
+	    	// std::cout<<"Fit of degree "<<poly_degree<<" with error = "<<error<<std::endl;
 
 	    }while(error > tolerance);
 
@@ -235,6 +246,7 @@ public:
 
 		return res;
 	}
+
 };
 
 

src/regression/regression_test.cpp

@@ -14,8 +14,8 @@
 BOOST_AUTO_TEST_SUITE( Regression_test )
 
 
-/*
-BOOST_AUTO_TEST_CASE ( Regression_domain_initialization )
+
+BOOST_AUTO_TEST_CASE ( Regression_local )
 {
     Box<2,float> domain({0.0,0.0},{1.0,1.0});
     // Here we define the boundary conditions of our problem
@@ -44,13 +44,14 @@ BOOST_AUTO_TEST_CASE ( Regression_domain_initialization )
     }    
     vd.map();
 
-    double size = 0.075;
-    auto NN = vd.getCellList(size);
-    auto dom = new RegressionDomain<vectorType, decltype(NN)>(vd, {0.8,0.8}, size, NN);
+    auto it2 = vd.getDomainIterator();
+    auto NN = vd.getCellList(0.1);
+    /*
+    auto support = RegressionSupport<vectorType, decltype(NN)>(vd, it2, 10, N_PARTICLES, NN);
     
-    std::cout<<"Initialized domain with "<<dom->getNumParticles()<<" particles.\n";
-
-    auto model = new RegressionModel<2, 0, vectorType>(vd, dom, 6, 2.0);
+    std::cout<<"Initialized domain with "<<support.getNumParticles()<<" particles.\n";
+    
+    auto model = RegressionModel<2, 0, vectorType>(vd, dom, 6, 2.0);
 
     double max_err = -1.0;
     for(double x = 0.75; x < 0.85;x+=0.01)
@@ -73,15 +74,9 @@ BOOST_AUTO_TEST_CASE ( Regression_domain_initialization )
         check = false;
     std::cout<<"Max err = "<<max_err<<"\n";
     BOOST_TEST( check );
-
-    if(model)
-        delete model;
-    
-    if(dom)
-        delete dom;
-
+    */
 }
-*/
+
 
 
 BOOST_AUTO_TEST_CASE ( Regression_without_domain_initialization)
@@ -93,52 +88,94 @@ BOOST_AUTO_TEST_CASE ( Regression_without_domain_initialization)
     Ghost<2,float> g(0.01);
 
     using vectorType = vector_dist<2,float, aggregate<double> >;
-    vectorType vd(2048,domain,bc,g);
+    
+    vectorType vd_orig(1024,domain,bc,g);
+
+    vectorType vd_test(vd_orig.getDecomposition(), 200);
+    
+    Vcluster<> & v_cl = create_vcluster();
+    long int N_prc = v_cl.getProcessingUnits();
+
+    if (v_cl.getProcessUnitID() == 0)
+        std::cout<<"Running regression test with "<<N_prc<<" procs.\n";
+
     // the scalar is the element at position 0 in the aggregate
     const int scalar = 0;
 
-    auto it = vd.getDomainIterator();
-    while (it.isNext())
+    // Creating a grid with synthetic data
     {
-        auto key = it.get();
-        double posx = (double)rand() / RAND_MAX;
-        double posy = (double)rand() / RAND_MAX;
-
-        vd.getPos(key)[0] = posx;
-        vd.getPos(key)[1] = posy;
+        auto it = vd_orig.getDomainIterator();
+        while (it.isNext())
+        {
+            auto key = it.get();
+            double posx = (double)rand() / RAND_MAX;
+            double posy = (double)rand() / RAND_MAX;
+
+            vd_orig.getPos(key)[0] = posx;
+            vd_orig.getPos(key)[1] = posy;
+
+            // Use synthetic data x*y
+            vd_orig.template getProp<scalar>(key) = sin(posx*posy);
+            ++it;
+        }    
+        vd_orig.map();
+    }
 
-        // Use synthetic data x*y
-        vd.template getProp<scalar>(key) = sin(posx*posy);
-        ++it;
-    }    
-    vd.map();
+    // Creating test points
+    {
+        auto it = vd_test.getDomainIterator();
+        while (it.isNext())
+        {
+            auto key = it.get();
+            double posx = (double)rand() / RAND_MAX;
+            double posy = (double)rand() / RAND_MAX;
+
+            vd_test.getPos(key)[0] = posx;
+            vd_test.getPos(key)[1] = posy;
+
+            vd_test.template getProp<scalar>(key) = 0.0;
+            
+            ++it;
+        }    
+        vd_test.map();
+    }
     
-    auto model = RegressionModel<2, 0>(vd, 1e-6);
+    auto model = RegressionModel<2, 0>(vd_orig, 1e-6);
 
     double max_err = -1.0;
-    for(double x = 0.75; x < 0.85;x+=0.01)
+    // Checking the error
     {
-        for(double y = 0.75; y < 0.85; y+=0.01)
+        auto it = vd_test.getDomainIterator();
+        while (it.isNext())
         {
-            Point<2, double> pos{x,y};
+            auto key = it.get();
+            Point<2, double>  pos = {vd_test.getPos(key)[0], vd_test.getPos(key)[1]};
             double val = model.eval(pos);
-            double actual = sin(x*y);
+            double actual = sin(pos[0]*pos[1]);
             double err = std::abs(actual - val);
             if (err > max_err) max_err = err;
-        }
+           
+            vd_test.template getProp<scalar>(key) = val;
+        
+            ++it;
+        }    
+        vd_test.ghost_get<scalar>();
     }
 
+    v_cl.max(max_err);
+	v_cl.execute();
+	if (v_cl.getProcessUnitID() == 0)
+		std::cout << "Maximum error: " << max_err << "\n";
+
     double tolerance = 1e-5;
     bool check;
     if (std::abs(max_err) < tolerance)
         check = true;
     else
         check = false;
-    std::cout<<"Max err = "<<max_err<<"\n";
+
     BOOST_TEST( check );
 
-    // if(model)
-    //     delete model;
 
 }