Create Pointset on a Sphere or a Klein bottle

genPointsShpereKleinbottle.cpp

This code is based on the code of Clement Maria from Inria. You basically have two functions generating either a sphere (of the dimension you like) or a Klein bottle. Then, it translates it randomly in a space of the dimension you like.

#include <cmath>
#include <vector>
#include <iostream>
#include <fstream>

/*
* norm squared
*/
double normPoint(std::vector<double>& point, int dim) {
  double norm = 0;
  for (int i = 0; i < dim; i++)
    norm += point[i] * point[i];
  return norm;
}
double ranf() { return (double)rand() / (double)RAND_MAX; }
double randGauss(double sigma) {
  double x1, x2, w, y1, y2;
  do {
    x1 = 2.0 * ranf() - 1.0;
    x2 = 2.0 * ranf() - 1.0;
    w = x1 * x1 + x2 * x2;
  } while (w >= 1.0);

  w = sqrt((-2.0 * log(w)) / w);
  y1 = x1 * w;
  y2 = x2 * w;

  return y1 * sqrt(sigma);
}
/** Generate nb_P points uniformly random on an (intr_dim)-sphere
  * embedded in (emb_dim)-dimensional Euclidean space.
  *
  * noise is the amplitude of the gaussian noise added to each coordinate
  */
void generate_rand_sphere(int intr_dim,
                          int emb_dim,
                          int nb_P,
                          double noise,
                          std::string filename) {
  std::vector<std::vector<double> > pointset;
  pointset.reserve(nb_P);
  // Generate uniform random points on a intr_dim sphere contained exactly in
  // the intr_dim+1 first coordinates
  std::cout << "Generate Sphere: \n";

  for (int i = 0; i < nb_P; ++i) {
    std::vector<double> v(intr_dim + 1, 0.);
    double norm = 0.;
    for (int d = 0; d < intr_dim + 1; ++d)  // all other coordinates are 0
    {
      v[d] = randGauss(1.);
      norm += v[d] * v[d];
    }
    norm = sqrt(norm);
    for (int d = 0; d < intr_dim + 1; ++d) {
      v[d] = v[d] / norm + randGauss((double)1) * noise;
    }
    pointset.push_back(v);
  }
  // Scale and translate anywhere in space
  double trans_ampl = 100;
  double scale_factor = randGauss(1.) * 0.5;

  std::vector<double> translation_vector(emb_dim);
  for (int d = 0; d < emb_dim; ++d) {
    translation_vector[d] = randGauss(1.) * trans_ampl;
  }

  std::ofstream outfile(filename.c_str(), std::ios::trunc);
  if (!outfile) {
    std::cerr << "Can't open " << filename << std::endl;
  }

  std::cout << "Write in file: \n";

  for (int i = 0; i < nb_P; ++i) {
    for (int d = 0; d < intr_dim + 1; ++d) {
      outfile << scale_factor* pointset[i][d] + translation_vector[d] << " ";
    }
    for (int d = intr_dim + 1; d < emb_dim; ++d) {
      outfile << translation_vector[d] << " ";
    }

    outfile << " \n";
  }
  outfile.close();

  std::cout << "Generate " << nb_P << " points on a sphere. \n";
}

/** Generate "about" nb_P points on a Klein bottle
  * embedded in (emb_dim)-dimensional Euclidean space.
  *
  * noise is the amplitude of the gaussian noise added to each
  * coordinate
  */
void generate_Klein_bottle(int emb_dim,
                           int nb_P,
                           double noise,
                           std::string filename) {
  double sqrt_nb_P = sqrt((double)nb_P);
  double Pi = 3.141592;
  double a = .5;
  int count = 0;

  std::vector<std::vector<double> > pointset;
  pointset.reserve(nb_P);

  std::cout << "Generate Klein bottle: \n";

  for (double u = 0; u < 2 * Pi - 0.00001; u += 2 * Pi / sqrt_nb_P) {
    for (double v = 0; v < 2 * Pi - 0.00001; v += 2 * Pi / sqrt_nb_P) {
      std::vector<double> point(5, 0.);

      point[0] =
          (a + cos(u / 2.) * sin(v) - sin(u / 2.) * sin(2. * v)) * cos(u) +
          randGauss((double)1) * noise;
      point[1] =
          (a + cos(u / 2.) * sin(v) - sin(u / 2.) * sin(2. * v)) * sin(u) +
          randGauss((double)1) * noise;
      point[2] = sin(u / 2.) * sin(v) + cos(u / 2.) * sin(2. * v) +
                 randGauss((double)1) * noise;
      point[3] = sin(u) + randGauss((double)1) * noise;
      point[4] = cos(v) + randGauss((double)1) * noise;

      ++count;

      pointset.push_back(point);
    }
  }

  // Scale and translate anywhere in space
  double trans_ampl = 100;
  double scale_factor = randGauss(1.) * 0.5;

  std::vector<double> translation_vector(emb_dim);
  for (int d = 0; d < emb_dim; ++d) {
    translation_vector[d] = randGauss(1.) * trans_ampl;
  }

  std::ofstream outfile(filename.c_str(), std::ios::trunc);
  if (!outfile) {
    std::cerr << "Can't open " << filename << std::endl;
  }

  std::cout << "Write in file: \n";

  for (int i = 0; i < count; ++i) {
    for (int d = 0; d < 5; ++d) {
      outfile << scale_factor* pointset[i][d] + translation_vector[d] << " ";
    }
    for (int d = 5; d < emb_dim; ++d) {
      outfile << translation_vector[d] << " ";
    }

    outfile << " \n";
  }
  outfile.close();

  std::cout << "Generate " << count << " points on a Klein bottle. \n";
}

int main() {

  generate_Klein_bottle(
      10000, 100000, 0.05, "Klein_nb_P100000_dim10000.txt");

  generate_rand_sphere(
      6, 10000, 100000, 0.1, "Sphere_nb_P100000_dim10000.txt");

  generate_Klein_bottle(
      100, 10000000, 0.05, "Klein_nb_P10000000_dim100.txt");

  generate_rand_sphere(
      6, 100, 10000000, 0.1, "Sphere_nb_P10000000_dim100.txt");

  generate_Klein_bottle(
      10000, 1000000, 0.05, "Klein_nb_P1000000_dim10000.txt");

  generate_rand_sphere(
      6, 10000, 1000000, 0.1, "Sphere_nb_P1000000_dim10000.txt");

  return 0;
}

Have questions about this code? Comments? Did you find a bug? Let me know!😀
Page created by G. (George) Samaras (DIT)