added shit

Makefile.test

@@ -3,13 +3,17 @@ CUDA_INCLUDE = -I./cuda_util_include #most of it comes with nvcc, this is a conv
 BOOST_INCLUDE = -I/usr/local/include/
 DELAUNAY_INCLUDE = -I./gDel3D/GPU/
 CUB_INCLUDE = -I/data/tools/cub-1.8.0/
-
+EIGEN_INCLUDE = -I/data/tools/lorm/3rdparty/include/eigen3/
 CFLAGS=-I. -dc $(CUDA_INCLUDE) -g -G 
 
 BOOST_LINKS = -L/usr/local/lib -lboost_unit_test_framework
 CUDA_LINKS = -L/usr/local/cuda-10.0/lib64/ -lcusolver -lcudart -lcublas -lcuda -lcudadevrt
 
-LFLAGS=-L. $(CUDA_LINKS) $(BOOST_LINKS) 
+LFLAGS=-L. $(CUDA_LINKS) $(BOOST_LINKS)
+
+LINK_LORM = -L/data/tools/lorm/build/lib -llorm -lLORM_NFFTwrapper_obj
+LINK_NFFT3 = -L/usr/local/lib -lnfft3
+
 
 eval: GpuDelaunay.o KerDivision.o Splaying.o Star.o PredWrapper.o RandGen.o InputCreator.o predicates.o ThrustWrapper.o KerPredicates.o CRotationGrid.o DelaunayChecker.o
 	$(NVCC) --gpu-architecture=sm_61 --include-path=./,./gDel3D/GPU/ CRotationGrid.o InputCreator.o RandGen.o DelaunayChecker.o KerPredicates.o Splaying.o Star.o predicates.o ThrustWrapper.o PredWrapper.o KerDivision.o GpuDelaunay.o TransformationGrid.cu -o eval
@@ -24,16 +28,20 @@ ccl2: CMrcReader.o
 	$(NVCC) --gpu-architecture=sm_61 --device-link Sampler.o PdbReader.o Particles.o Labeler.o Kernels.o Density.o ccl_test.o --output-file link.o
 	g++ Sampler.o PdbReader.o Particles.o Labeler.o Kernels.o Density.o CMrcReader.o ccl_test.o link.o -o ccl_test -L. $(CUDA_LINKS) $(BOOST_LINKS)
 
-weights:
+weights: 
 	g++ -c -Wall -I. -I/usr/local/cuda-10.0/include/ SO3_weights.cpp
 	g++ SO3_weights.o cg.o -lm -o weights
 	
-grid:
-	g++ -O2 -o grid search_grid_alpha.cpp
 
 sort:
 	$(NVCC) -g -G -dc --gpu-architecture=sm_61 --include-path=./ $(CUB_INCLUDE) example_block_radix_sort.cu
 	$(NVCC)  example_block_radix_sort.o -lcudadevrt -o sort
+
+grid: QuadratureErrorSO3.o
+	g++ $(EIGEN_INCLUDE) $(LINK_LORM) $(LINK_NFFT3) QuadratureErrorSO3.o -O2 -o grid search_grid_alpha.cpp
+
+QuadratureErrorSO3.o: /data/tools/lorm/src/ObjectiveFunctions/QuadratureErrorSO3.cpp
+	g++ $(EIGEN_INCLUDE) $(LINK_LORM) -O2 -c -o $@ $<
 	
 ccl_test.o: ccl_test.cu
 	$(NVCC) -c -o $@ $< $(CFLAGS) $(BOOST_INCLUDE) 

search_grid_alpha.cpp

@@ -62,6 +62,12 @@
 #include <string>
 #include <limits>
 
+#include "/data/tools/lorm/src/Manifolds/RotationGroup/RotationGroupPointArray.hpp"
+#include "/data/tools/lorm/src/ObjectiveFunctions/QuadratureErrorSO3.hpp"
+#include "/data/tools/lorm/src/Manifolds/Sphere/Sphere.hpp"
+#include "/data/tools/lorm/src/OptimizationMethods/ArmijoMethod.hpp"
+#include "/data/tools/lorm/src/OptimizationMethods/ConjugateGradientMethod.hpp"
+
 // Windows doesn't define M_PI in the standard header?
 #if !defined(M_PI)
 #define M_PI 3.1415926535897932384626433832795028841971694
@@ -261,111 +267,170 @@ int find_max(float d){
 	return std::max(2*((int)f),2*((int)(f-0.5))+1);
 }
 
-int main(int argc, char* argv[], char*[]) {
-  bool euler = false;
-//  if (argc > 1 && string(argv[1]) == "-e")
-//    euler = true;
-//  assert(cin.good());
-  string line;
-//  while (cin.peek() == '#') {
-//    getline(cin, line);
-//    cout << line << endl;
-//  }
-  float d0 = 0.15167;
-  float d = d0;
-  int max_i = find_max(d0);
-  int umax = max_i%2==0?max_i-1:max_i;
-  int emax = max_i%2==0?max_i:max_i-1;
-
-  std::list<std::tuple<int,int,int>> u_indeces;
-  std::list<std::tuple<int,int,int>> e_indeces;
-
-  printf("Umax: %i\n",umax);
-  printf("Emax: %i\n",emax);
-
-  for (int i=1;i<umax+1;i+=2){
-	 for (int j=1;j<=i+1;j+=2){
-		 for (int k=1;k<=j+1;k+=2){
-			 if ((i+j+k)*d/2<1){
-				 u_indeces.push_back(std::make_tuple(i,j,k));
-				 printf("%i %i %i\n",i,j,k);}
+int get_multiplicity(int i,int j,int k){
+	//really lazy, but not performance critical, so pardon me
+	if (i%2 != 0){
+		//uneven index tuple
+		if (i != j && j != k && i != k) return 48;
+		if ((i != j && j == k) || (i != k && j == i) || (j != k && k == i)) return 24;
+		return 8;
+	}
+	//even index tuple
+	int number_zeros = 0;
+	if (i == 0) number_zeros +=1;
+	if (j == 0) number_zeros +=1;
+	if (k == 0) number_zeros +=1;
+	if (number_zeros == 3) return 1;
+	if (number_zeros == 2) return 6;
+	if (number_zeros == 1) {
+		if (i == j || j ==k || i==k )return 12;
+		return 24;
+	}
+	if (number_zeros == 0){
+		if (i != j && j != k && i != k) return 48;
+		if ((i != j && j == k) || (i != k && j == i) || (j != k && k == i)) return 24;
+		return 8;
+	}
+	//if 0 --> problem
+	return 0;
+}
+
+void create_orientations(float d,int degree_refinement){
+	  int max_i = find_max(d);
+	  int umax = max_i%2==0?max_i-1:max_i;
+	  int emax = max_i%2==0?max_i:max_i-1;
+	  bool euler = true;
+
+	  typedef std::tuple<int,int,int,int> four_int;
+	  std::list<four_int> u_indeces;
+	  std::list<four_int> e_indeces;
+	  std::list<four_int> indeces;
+//
+//	  printf("Umax: %i\n",umax);
+//	  printf("Emax: %i\n",emax);
+
+	  for (int i=1;i<umax+1;i+=2){
+		 for (int j=1;j<=i+1;j+=2){
+			 for (int k=1;k<=j+1;k+=2){
+				 if ((i+j+k)*d/2<1)u_indeces.push_back(std::make_tuple(i,j,k,get_multiplicity(i,j,k)));
+			 }
 		 }
-	 }
-  }
+	  }
 
-  for (int i=0;i<emax+1;i+=2){
-	 for (int j=0;j<=i+1;j+=2){
-		 for (int k=0;k<=j+1;k+=2){
-			 if ((i+j+k)*d/2<1){
-				 e_indeces.push_back(std::make_tuple(i,j,k));
-				 printf("%i %i %i\n",i,j,k);
+	  for (int i=0;i<emax+1;i+=2){
+		 for (int j=0;j<=i+1;j+=2){
+			 for (int k=0;k<=j+1;k+=2){
+				 if ((i+j+k)*d/2<1)e_indeces.push_back(std::make_tuple(i,j,k,get_multiplicity(i,j,k)));
 			 }
 		 }
-	 }
-  }
+	  }
 
+	  for (four_int _tuple:u_indeces) indeces.push_back(_tuple);
+	  if (degree_refinement ==0) for (four_int _tuple:e_indeces) indeces.push_back(_tuple);
 
-  std::exit(0);
-  assert(cin.good());
-  getline(cin, line);
-  assert(line == "format grid");
-  cout << "format " << (euler ? "euler" : "quaternion") << endl;
-  double delta, sigma, maxrad, coverage;
-  size_t ncell, ntot, nent;
-  cin >> delta >> sigma >> ntot >> ncell >> nent >> maxrad >> coverage;
-  // Use extra digit of precision with weights and radii.  This also
-  // triggers a memory minimizing expansion.
-  const bool fine = delta < 0.05;
-  cout << ntot << " " << fixed
-       << setprecision(fine ? 3 : 2) << maxrad << " "
-       << setprecision(5) << coverage << endl;
-  PackSet s;
-  size_t ncell1 = 0;
-  for (size_t n = 0; n < nent; ++n) {
-    int k, l, m;
-    size_t mult;
-    double r, w;
-    assert(cin.good());
-    cin >> k >> l >> m >> w >> r >> mult;
-    Permute p(Triple(k, l, m));
-    assert(mult == p.Number());
-    for (size_t i = 0; i < mult; ++i) {
-      Triple t = p.Member(i);
-      s.Add(Quaternion(1.0,
-		       pind(0.5 * t.a, delta, sigma),
-		       pind(0.5 * t.b, delta, sigma),
-		       pind(0.5 * t.c, delta, sigma)),
-	    w);
-    }
-    ncell1 += mult;
-    if (fine) {
-      // Skip n = 0; that's already included.
-      for (size_t n = 1; n < 24; ++n) {
-	Quaternion q(CubeSyms[n][0], CubeSyms[n][1],
-		     CubeSyms[n][2], CubeSyms[n][3]);
-	for (size_t i = 0; i < mult; ++i)
-	  s.Add(q.Times(s.Orientation(i)), s.Weight(i));
-      }
-      s.Print(cout, euler, fine ? 7 : 6);
-      s.Clear();
-    }
-  }
-  assert(cin.good());
-  assert(ncell1 == ncell);
-  if (!fine) {
-    size_t nc = s.Number();
-    assert(nc == ncell);
-    for (size_t n = 1; n < 24; ++n) {
-      Quaternion q(CubeSyms[n][0], CubeSyms[n][1],
-		   CubeSyms[n][2], CubeSyms[n][3]);
-      for (size_t i = 0; i < nc; ++i)
-	s.Add(q.Times(s.Orientation(i)), s.Weight(i));
+	  indeces.sort([](const four_int & tuple1, const four_int & tuple2)
+	  {
+		  if(std::get<0>(tuple1) != std::get<0>(tuple2)) return std::get<0>(tuple1) < std::get<0>(tuple2);
+		  if(std::get<1>(tuple1) != std::get<1>(tuple2)) return std::get<1>(tuple1) < std::get<1>(tuple2);
+		  return std::get<2>(tuple1) < std::get<2>(tuple2);
+	  });
+	  size_t ncell, ntot, nent;
+	  nent = indeces.size();
+	  ncell = 0;
+	  for(four_int _tuple:indeces) ncell += std::get<3>(_tuple);
+	//  for(four_int _tuple:indeces) printf("%i %i %i %i\n",std::get<0>(_tuple),std::get<1>(_tuple),std::get<2>(_tuple),std::get<3>(_tuple));
+	  ntot = 24*ncell;
+
+	  printf("#%i %u %u %u\n",degree_refinement,ntot,ncell,nent);
+	//  std::exit(0);
+	//  double delta, sigma, maxrad, coverage;
+
+	//  cin >> delta >> sigma >> ntot >> ncell >> nent >> maxrad >> coverage;
+	//  // Use extra digit of precision with weights and radii.  This also
+	//  // triggers a memory minimizing expansion.
+	  const bool fine = d < 0.05;
+	//  cout << ntot << " " << fixed
+	//       << setprecision(fine ? 3 : 2) << maxrad << " "
+	//       << setprecision(5) << coverage << endl;
+	  PackSet s;
+	  size_t ncell1 = 0;
+	  for(four_int _tuple:indeces){
+	//  for (size_t n = 0; n < nent; ++n) {
+		  int k = std::get<0>(_tuple);
+		  int l = std::get<1>(_tuple);
+		  int m = std::get<2>(_tuple);
+		  size_t mult = (size_t) std::get<3>(_tuple);;
+	//    double r, w;
+	//    cin >> k >> l >> m >> w >> r >> mult;
+		  Permute p(Triple(k, l, m));
+	    for (size_t i = 0; i < mult; ++i) {
+	      Triple t = p.Member(i);
+	      s.Add(Quaternion(1.0,
+			       pind(0.5 * t.a, d, 0.0),
+			       pind(0.5 * t.b, d, 0.0),
+			       pind(0.5 * t.c, d, 0.0)),
+		    1.0);
+	    }
+	//    ncell1 += mult;
+	    if (fine) {
+	//      // Skip n = 0; that's already included.
+	      for (size_t n = 1; n < 24; ++n) {
+		Quaternion q(CubeSyms[n][0], CubeSyms[n][1],
+			     CubeSyms[n][2], CubeSyms[n][3]);
+		for (size_t i = 0; i < mult; ++i)
+		  s.Add(q.Times(s.Orientation(i)), s.Weight(i));
+	      }
+	      s.Print(cout, euler, fine ? 7 : 6);
+	      s.Clear();
+	    }
+	  }
+	//  assert(cin.good());
+	//  assert(ncell1 == ncell);
+	  if (!fine) {
+	    size_t nc = s.Number();
+	    assert(nc == ncell);
+	    for (size_t n = 1; n < 24; ++n) {
+	      Quaternion q(CubeSyms[n][0], CubeSyms[n][1],
+			   CubeSyms[n][2], CubeSyms[n][3]);
+	      for (size_t i = 0; i < nc; ++i)
+		s.Add(q.Times(s.Orientation(i)), s.Weight(i));
+	    }
+	    assert(s.Number() == ntot);
+	    s.Print(cout, euler, fine ? 7 : 6);
+	    s.Clear();
+	  }
+
+}
+
+inline int ipow(int base, int exp)
+{
+    int result = 1;
+    for (;;)
+    {
+        if (exp & 1)
+            result *= base;
+        exp >>= 1;
+        if (!exp)
+            break;
+        base *= base;
     }
-    assert(s.Number() == ntot);
-    s.Print(cout, euler, fine ? 7 : 6);
-    s.Clear();
-  }
-  return 0;
+
+    return result;
+}
+
+
+int main(int argc, char* argv[], char*[]) {
+	float d0 = 0.41422;
+	int degree_refinement = 2;
+	for (int i=0;i<=degree_refinement;i++) create_orientations(d0/(ipow(2,i)),i);
+
+
+	  // creates a new instance of the quadrature error function on the rotation
+	  // group SO(3) for polynomial degree N and M points with weights
+	int M=4, N=1;
+  	LORM::ObjectiveFunctions::QuadratureErrorSO3 error(M, N);
+
+	return 0;
 }
 
 void Quaternion::Print(ostream& s) const {
@@ -412,6 +477,9 @@ void Quaternion::PrintEuler(ostream& s) const {
     a = atan2(-m01, m11);
     c = 0;
   }
+  //make it so that
+  a += M_PI;
+  c += M_PI;
   s << fixed << setprecision(9) << setw(12) << a << " "
     << setw(12) << b << " " << setw(12) << c;
 
@@ -422,6 +490,6 @@ void Quaternion::PrintEuler(ostream& s) const {
 	  Times(Quaternion(cos(c/2), 0, 0, sin(c/2)))); // c about z
   // and check that q is parallel to *this.
   double t = abs(q.w * w + q.x * x + q.y * y + q.z * z);
-  assert(t > 1 - 16 * numeric_limits<double>::epsilon());
+//  assert(t > 1 - 16 * numeric_limits<double>::epsilon());
 #endif
 }