compute.cpp 4.62 KB
Newer Older
1 2
#include "compute.hpp"

Lars Kuehne's avatar
Lars Kuehne committed
3 4 5
#include <cassert>
#include <cmath>

6
#include <algorithm>
Lars Kuehne's avatar
Lars Kuehne committed
7 8 9 10 11 12
#include <chrono>
#include <fstream>
#include <iostream>
#include <string>

#include "picture.hpp"
13

Lars Kuehne's avatar
Lars Kuehne committed
14 15 16 17 18
constexpr float PI = 3.14159265358979323846264338327950288419716939f;
constexpr float TWO_PI = 2 * PI;

float compute(float alpha, float beta, const aligned_vector<float>& seed_x,
              const aligned_vector<float>& seed_y, int num_iterations,
19
              float threshold) {
20 21 22 23 24 25 26 27 28 29 30
  int num_seeds = seed_x.size();
  assert(seed_y.size() == num_seeds);

  aligned_vector<float> x = seed_x;
  aligned_vector<float> y = seed_y;

  float* xp = x.data();
  float* yp = y.data();

  float d = 0.0;

31
  for (int i = 0; i < num_iterations && d <= threshold; i++) {
32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
#pragma omp simd aligned(xp, yp : 64)
    for (int s = 0; s < num_seeds; s++) {
      yp[s] = yp[s] + beta * std::sin(TWO_PI * xp[s]);
    }
#pragma omp simd aligned(xp, yp : 64)
    for (int s = 0; s < num_seeds; s++) {
      xp[s] = xp[s] + alpha * std::sin(TWO_PI * yp[s]);
    }

    for (int s = 0; s < num_seeds; s++) {
      d = std::max(d, std::abs(yp[s]));
    }
  }

  return d;
}
48

Lars Kuehne's avatar
Lars Kuehne committed
49 50 51 52
void compute_all(int num_iterations, float threshold, float alphamin,
                 float alphamax, int alpha_num_intervals, float betamin,
                 float betamax, int beta_num_intervals, int num_seedpoints,
                 bool output_csv, std::vector<float> seedpoints) {
53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72
  // these are computed
  int alpha_num_params = alpha_num_intervals + 1;
  float alpha_interval_size = (alphamax - alphamin) / (alpha_num_intervals);
  int beta_num_params = beta_num_intervals + 1;
  float beta_interval_size = (betamax - betamin) / (beta_num_intervals);

  // fill prametervectors alpha and beta
  aligned_vector<float> alphas(alpha_num_params);
  float* alphasp = alphas.data();
#pragma omp simd aligned(alphasp : 64)
  for (int i = 0; i < alpha_num_params; i++) {
    alphasp[i] = alphamin + i * alpha_interval_size;
  }
  aligned_vector<float> betas(beta_num_params);
  float* betasp = betas.data();
#pragma omp simd aligned(betasp : 64)
  for (int i = 0; i < beta_num_params; i++) {
    betasp[i] = betamin + i * beta_interval_size;
  }

Lars Kuehne's avatar
Lars Kuehne committed
73
  // Initialization and output of seedpoints
74 75 76 77
  std::cout << "Following seedpoints are used for computation:" << std::endl;
  aligned_vector<float> x_start(num_seedpoints);
  aligned_vector<float> y_start(num_seedpoints);
  for (int i = 1; i < num_seedpoints - seedpoints.size() + 1; ++i) {
Lars Kuehne's avatar
Lars Kuehne committed
78 79 80
    x_start[i - 1] =
        0.5f * static_cast<float>(i) / (num_seedpoints - seedpoints.size() + 1);
    //    y_start[i - 1] = 0;
81 82 83 84
    std::cout << x_start[i - 1] << ", ";
  }
  for (int i = 0; i < seedpoints.size(); ++i) {
    x_start[num_seedpoints - seedpoints.size() + i] = seedpoints[i];
Lars Kuehne's avatar
Lars Kuehne committed
85
    //    y_start[num_seedpoints - seedpoints.size() + i] = 0;
86 87 88 89 90 91 92 93 94
    std::cout << x_start[num_seedpoints - seedpoints.size() + i] << ", ";
  }
  std::cout << '\n';

  // Initialization pixel values and color vectors
  aligned_vector<float> result(alpha_num_params * beta_num_params);

  auto time_start = std::chrono::system_clock::now();

Lars Kuehne's avatar
Lars Kuehne committed
95
  // Computation
96 97 98 99 100 101 102 103 104 105 106 107 108
#pragma omp parallel for schedule(dynamic)
  for (int b = beta_num_params - 1; b >= 0; b--) {
    for (int a = 0; a < alpha_num_params; a++) {
      result[(beta_num_params - b - 1) * alpha_num_params + a] = compute(
          alphas[a], betas[b], x_start, y_start, num_iterations, threshold);
    }
  }
  auto time_end = std::chrono::system_clock::now();
  float elapsed_seconds =
      std::chrono::duration<float>(time_end - time_start).count();
  std::cout << "TIME for computation: " << elapsed_seconds << std::endl;

  time_start = std::chrono::system_clock::now();
Lars Kuehne's avatar
Lars Kuehne committed
109 110
  write_png("picture.png", result.data(), threshold, alpha_num_params,
            beta_num_params);
111
  time_end = std::chrono::system_clock::now();
Lars Kuehne's avatar
Lars Kuehne committed
112
  elapsed_seconds = std::chrono::duration<float>(time_end - time_start).count();
113 114
  std::cout << "TIME for picture: " << elapsed_seconds << std::endl;

Lars Kuehne's avatar
Lars Kuehne committed
115
  // Generate output
116 117
  if (output_csv) {
    time_start = std::chrono::system_clock::now();
Lars Kuehne's avatar
Lars Kuehne committed
118
    // Output result into .csv
119 120 121 122 123 124 125 126 127 128
    std::string file_result = "result.csv";
    std::ofstream ostrm_csv(file_result);
    ostrm_csv << "alpha beta value\n";
    for (int b = beta_num_params - 1; b >= 0; b--) {
      for (int a = 0; a < alpha_num_params; a++) {
        ostrm_csv << alphas[a] << ' ' << betas[b] << ' '
                  << result[(beta_num_params - b - 1) * alpha_num_params + a]
                  << '\n';
      }
    }
Lars Kuehne's avatar
Lars Kuehne committed
129 130 131 132
    time_end = std::chrono::system_clock::now();
    elapsed_seconds =
        std::chrono::duration<float>(time_end - time_start).count();
    std::cout << "TIME for csv: " << elapsed_seconds << std::endl;
133 134
  }
}