diff --git a/variants/variant_AVX512_1_flopcount.cpp b/variants/variant_AVX512_1_flopcount.cpp deleted file mode 100644 index 718ee9f..0000000 --- a/variants/variant_AVX512_1_flopcount.cpp +++ /dev/null @@ -1,568 +0,0 @@ -// # This file attempts to use AVX-512 intrinsics to speed up computation. -// # Predecessor: variant_inlining_7.cpp -#ifdef __AVX512F__ - -#include "../headers/board_reduced.hpp" -#include "../headers/constants.hpp" -#include "../headers/pedestrian_reduced.hpp" - -#include "immintrin.h" -#include "sleef.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#define SEED 42 -#define BOARD_WIDTH 10.0 -#define BOARD_LENGTH 50.0 -#define BOARD_LENGTH_MINUS_ONE BOARD_LENGTH - 1.0f -#define BOARD_LENGTH_MINUS_ONE_POINT_FIVE BOARD_LENGTH - 1.5f -#define WALL_EPS 1e-3 -#define BOARD_WIDTH_MINUS_WALL_EPS 9.999 - -#ifndef RANDOM_SEED -#define RANDOM_SEED 0 -#endif - -#ifndef BENCHMARK -#define BENCHMARK 0 -#endif - -using namespace std; - -Pedestrian *init_pedestrians(int n_p) { -#if RANDOM_SEED - random_device rd; - mt19937 gen(rd()); -#else - mt19937 gen(SEED); -#endif - - normal_distribution velocity_gen(v0_mean, v0_std_dev); - uniform_real_distribution y_gen(0.3, BOARD_WIDTH - 0.3); - uniform_real_distribution x_gen_left(2.0, BOARD_LENGTH / 2.0); - uniform_real_distribution x_gen_right(BOARD_LENGTH / 2.0, BOARD_LENGTH - 2.0); - - Pedestrian *P = new Pedestrian[n_p]; - - int n_p_half = n_p / 2; - - for (int i = 0; i < n_p_half; i++) { - double x = x_gen_left(gen); - double y = y_gen(gen); - - double start_pos_x = x; - double start_pos_y = y; - - double destination_x = BOARD_LENGTH + 30.0; - double destination_y = y; - - double init_e_x = destination_x - start_pos_x; - double init_e_y = destination_y - start_pos_y; - double init_e_norm = sqrt(init_e_x * init_e_x + init_e_y * init_e_y); - double init_e_unit_x = init_e_x / init_e_norm; - double init_e_unit_y = init_e_y / init_e_norm; - - double desired_v0; - do { - desired_v0 = velocity_gen(gen); - } while (desired_v0 <= 0.1); - - double v_max = 1.3 * desired_v0; - - P[i] = Pedestrian(0, 0, 0, start_pos_x, start_pos_y, init_e_unit_x, init_e_unit_y, v_max, desired_v0, destination_x, destination_y); - } - - for (int i = n_p_half; i < n_p; i++) { - double x = x_gen_right(gen); - double y = y_gen(gen); - - double start_pos_x = x; - double start_pos_y = y; - double destination_x = -30; - double destination_y = y; - - double init_e_x = destination_x - start_pos_x; - double init_e_y = destination_y - start_pos_y; - double init_e_norm = sqrt(init_e_x * init_e_x + init_e_y * init_e_y); - double init_e_unit_x = init_e_x / init_e_norm; - double init_e_unit_y = init_e_y / init_e_norm; - - double desired_v0; - do { - desired_v0 = velocity_gen(gen); - } while (desired_v0 <= 0.1); - - double v_max = 1.3 * desired_v0; - - P[i] = Pedestrian(1, 0, 0, start_pos_x, start_pos_y, init_e_unit_x, init_e_unit_y, v_max, desired_v0, destination_x, destination_y); - } - - return P; -} - -Board init_board_without_obstacles() { return Board(BOARD_LENGTH, BOARD_WIDTH, 0, nullptr); } - -int run_simulation(Pedestrian *pedestrians, Board &board, int n_p, int n_iter, const char *dump_path) { - -#if BENCHMARK - (void)dump_path; -#endif - -#if !BENCHMARK - const int dump_interval = 1; - - vector> samples; - - if (dump_path) { - samples.reserve(static_cast(n_iter / dump_interval + 1) * n_p); - - for (int i = 0; i < n_p; i++) { - Pedestrian &p = pedestrians[i]; - samples.emplace_back(0, i, p.curr_pos_x, p.curr_pos_y, p.e_x, p.e_y, p.color); - } - } -#endif - - // # ADDs: n_iter * n_p * 4 + - // # n_iter * n_p * ceil(((double) (n_p - 7)) / 8) * (1 / (ceil(((double) (n_p - 7)) / 8))) * 7 * 14 + - // # n_iter * n_p * ceil(((double) (n_p - 7)) / 8) * (1 - 1 / (ceil(((double) (n_p - 7)) / 8))) * 224 + - // # n_iter * n_p * (n_p - floor((double) n_p / 8) * 8) * 14 + - // # n_iter * n_p * board.n_obstacles * 12 + - // # n_iter * n_p * 7 - - // # MULTs: n_iter * n_p * 4 + - // # n_iter * n_p * ceil(((double) (n_p - 7)) / 8) * (1 / (ceil(((double) (n_p - 7)) / 8))) * 7 * 24 + - // # n_iter * n_p * ceil(((double) (n_p - 7)) / 8) * (1 - 1 / (ceil(((double) (n_p - 7)) / 8))) * 352 + - // # n_iter * n_p * (n_p - floor((double) n_p / 8) * 8) * 24 + - // # n_iter * n_p * board.n_obstacles * 19 + - // # n_iter * n_p * 17 - - // # DIVs: n_iter * n_p * ceil(((double) (n_p - 7)) / 8) * (1 / (ceil(((double) (n_p - 7)) / 8))) * 7 + - // # n_iter * n_p * ceil(((double) (n_p - 7)) / 8) * (1 - 1 / (ceil(((double) (n_p - 7)) / 8))) * 8 + - // # n_iter * n_p * (n_p - floor((double) n_p / 8) * 8) + - // # n_iter * n_p * board.n_obstacles + - // # n_iter * n_p * 2 - - // # SQRTs: n_iter * n_p * ceil(((double) (n_p - 7)) / 8) * (1 / (ceil(((double) (n_p - 7)) / 8))) * 7 * 3 + - // # n_iter * n_p * ceil(((double) (n_p - 7)) / 8) * (1 - 1 / (ceil(((double) (n_p - 7)) / 8))) * 24 + - // # n_iter * n_p * (n_p - floor((double) n_p / 8) * 8) * 3 + - // # n_iter * n_p * board.n_obstacles + - // # n_iter * n_p * 3 - - // # EXPs: n_iter * n_p * ceil(((double) (n_p - 7)) / 8) * (1 / (ceil(((double) (n_p - 7)) / 8))) * 7 + - // # n_iter * n_p * (n_p - floor((double) n_p / 8) * 8) + - // # n_iter * n_p * board.n_obstacles - - - for (int t = 0; t < n_iter; t++) { - for (int i = 0; i < n_p; i++) { - Pedestrian &p = pedestrians[i]; - - // # 2 MULTs - double desired_velocity_x = p.e_x * p.desired_v0; - double desired_velocity_y = p.e_y * p.desired_v0; - - // # 2 ADDs, 2 MULTs - double f_alpha_dt_x = (desired_velocity_x - p.curr_vel_x) * dt_times_tau_alpha_inv; - double f_alpha_dt_y = (desired_velocity_y - p.curr_vel_y) * dt_times_tau_alpha_inv; - - __m512d desired_velocity_x_vec = _mm512_set1_pd(desired_velocity_x); - __m512d desired_velocity_y_vec = _mm512_set1_pd(desired_velocity_y); - __m512d f_alpha_dt_x_vec = _mm512_set1_pd(f_alpha_dt_x); - __m512d f_alpha_dt_y_vec = _mm512_set1_pd(f_alpha_dt_y); - - __m512d alpha_curr_pos_x_vec = _mm512_set1_pd(p.curr_pos_x); - __m512d alpha_curr_pos_y_vec = _mm512_set1_pd(p.curr_pos_y); - __m512d alpha_e_x_vec = _mm512_set1_pd(p.e_x); - __m512d alpha_e_y_vec = _mm512_set1_pd(p.e_y); - - int j; - - // # ############################################################################################################################################################################################ - - for (j = 0; j < n_p - 7; j += 8) { - - if (j <= i && i <= j + 7) { - - for (int jj = j; jj < j + 8; jj++) { - - if (jj == i) { - continue; - } - - // # 2 ADDs - double r_ab_x = p.curr_pos_x - pedestrians[jj].curr_pos_x; - double r_ab_y = p.curr_pos_y - pedestrians[jj].curr_pos_y; - - // # 1 SQRT, 2 MULTs, 1 ADD - double r_ab_norm = sqrt(r_ab_x * r_ab_x + r_ab_y * r_ab_y); - - // # 2 ADDs - double r_ab2_x = r_ab_x - pedestrians[jj].s_b_e_x; - double r_ab2_y = r_ab_y - pedestrians[jj].s_b_e_y; - - // # 1 SQRT, 2 MULTs, 1 ADD - double r_ab2_norm = sqrt(r_ab2_x * r_ab2_x + r_ab2_y * r_ab2_y); - - // # 2 ADDs, 1 MULT - double r_ab_norm_sum = r_ab_norm + r_ab2_norm; - double two_b_sq = r_ab_norm_sum * r_ab_norm_sum - pedestrians[jj].s_b_sq; - - // # 1 SQRT - double two_b = sqrt(two_b_sq); - - // # 5 MULTs, 1 DIV, 1 EXP - double scalar_mult_common = (r_ab_norm_sum / (two_b * r_ab_norm * r_ab2_norm)) * (NEG_V0_AB_div_TWO_SIGMA * exp(two_b * NEG_TWO_SIGMA_INV)); - - // # 2 MULTs - double scalar_mult_r_ab = scalar_mult_common * r_ab2_norm; - double scalar_mult_r_ab2 = scalar_mult_common * r_ab_norm; - - // # 4 MULTs, 2 ADDs - double f_alpha_beta_neg_x = r_ab_x * scalar_mult_r_ab + r_ab2_x * scalar_mult_r_ab2; - double f_alpha_beta_neg_y = r_ab_y * scalar_mult_r_ab + r_ab2_y * scalar_mult_r_ab2; - - double w_val; - - // # 3 MULTs, 1 ADD - double a = (p.e_x * f_alpha_beta_neg_x + p.e_y * f_alpha_beta_neg_y) * cos_phi_inv; - - // # 2 MULTs, 1 ADD - double b_sq = f_alpha_beta_neg_x * f_alpha_beta_neg_x + f_alpha_beta_neg_y * f_alpha_beta_neg_y; - - w_val = dt_neg; - // # 1 MULT - if (a >= 0 && (a * a > b_sq)) { - w_val = c_neg_times_dt; - } - - // # 2 MULTs, 2 ADDs - f_alpha_dt_x = f_alpha_dt_x + f_alpha_beta_neg_x * w_val; - f_alpha_dt_y = f_alpha_dt_y + f_alpha_beta_neg_y * w_val; - } - - } else { - __m512d beta_curr_pos_x_vec = _mm512_setr_pd(pedestrians[j].curr_pos_x, pedestrians[j + 1].curr_pos_x, pedestrians[j + 2].curr_pos_x, pedestrians[j + 3].curr_pos_x, pedestrians[j + 4].curr_pos_x, pedestrians[j + 5].curr_pos_x, pedestrians[j + 6].curr_pos_x, pedestrians[j + 7].curr_pos_x); - __m512d beta_curr_pos_y_vec = _mm512_setr_pd(pedestrians[j].curr_pos_y, pedestrians[j + 1].curr_pos_y, pedestrians[j + 2].curr_pos_y, pedestrians[j + 3].curr_pos_y, pedestrians[j + 4].curr_pos_y, pedestrians[j + 5].curr_pos_y, pedestrians[j + 6].curr_pos_y, pedestrians[j + 7].curr_pos_y); - - // # 16 ADDs - __m512d r_ab_x_vec = _mm512_sub_pd(alpha_curr_pos_x_vec, beta_curr_pos_x_vec); - __m512d r_ab_y_vec = _mm512_sub_pd(alpha_curr_pos_y_vec, beta_curr_pos_y_vec); - - // # 8 MULTs - __m512d r_ab_y_mul_r_ab_y_vec = _mm512_mul_pd(r_ab_y_vec, r_ab_y_vec); - - // # 8 MULTs, 8 ADDs - __m512d r_ab_norm_sq_vec = _mm512_fmadd_pd(r_ab_x_vec, r_ab_x_vec, r_ab_y_mul_r_ab_y_vec); - - // # 8 SQRTs - __m512d r_ab_norm_vec = _mm512_sqrt_pd(r_ab_norm_sq_vec); - - __m512d beta_s_b_e_x_vec = _mm512_setr_pd(pedestrians[j].s_b_e_x, pedestrians[j + 1].s_b_e_x, pedestrians[j + 2].s_b_e_x, pedestrians[j + 3].s_b_e_x, pedestrians[j + 4].s_b_e_x, pedestrians[j + 5].s_b_e_x, pedestrians[j + 6].s_b_e_x, pedestrians[j + 7].s_b_e_x); - - // # 8 ADDs - __m512d r_ab2_x_vec = _mm512_sub_pd(r_ab_x_vec, beta_s_b_e_x_vec); - - __m512d beta_s_b_e_y_vec = _mm512_setr_pd(pedestrians[j].s_b_e_y, pedestrians[j + 1].s_b_e_y, pedestrians[j + 2].s_b_e_y, pedestrians[j + 3].s_b_e_y, pedestrians[j + 4].s_b_e_y, pedestrians[j + 5].s_b_e_y, pedestrians[j + 6].s_b_e_y, pedestrians[j + 7].s_b_e_y); - - // # 8 ADDs - __m512d r_ab2_y_vec = _mm512_sub_pd(r_ab_y_vec, beta_s_b_e_y_vec); - - // # 8 MULTs - __m512d r_ab2_y_mul_r_ab2_y_vec = _mm512_mul_pd(r_ab2_y_vec, r_ab2_y_vec); - - // # 8 MULTs, 8 ADDs - __m512d r_ab2_norm_sq_vec = _mm512_fmadd_pd(r_ab2_x_vec, r_ab2_x_vec, r_ab2_y_mul_r_ab2_y_vec); - - // # 8 SQRTs - __m512d r_ab2_norm_vec = _mm512_sqrt_pd(r_ab2_norm_sq_vec); - - // # 8 ADDs - __m512d r_ab_norm_sum_vec = _mm512_add_pd(r_ab_norm_vec, r_ab2_norm_vec); - - __m512d beta_s_b_sq_vec = _mm512_setr_pd(pedestrians[j].s_b_sq, pedestrians[j + 1].s_b_sq, pedestrians[j + 2].s_b_sq, pedestrians[j + 3].s_b_sq, pedestrians[j + 4].s_b_sq, pedestrians[j + 5].s_b_sq, pedestrians[j + 6].s_b_sq, pedestrians[j + 7].s_b_sq); - - // # 8 MULTs, 8 ADDs - __m512d two_b_sq_vec = _mm512_fmsub_pd(r_ab_norm_sum_vec, r_ab_norm_sum_vec, beta_s_b_sq_vec); - - // # 8 SQRTs - __m512d two_b_vec = _mm512_sqrt_pd(two_b_sq_vec); - - // # 8 MULTs - __m512d two_b_mul_r_ab_norm_vec = _mm512_mul_pd(two_b_vec, r_ab_norm_vec); - - // # 8 MULTs - __m512d denom_vec = _mm512_mul_pd(two_b_mul_r_ab_norm_vec, r_ab2_norm_vec); - - // # 8 DIVs - __m512d scalar_mult_common_left_vec = _mm512_div_pd(r_ab_norm_sum_vec, denom_vec); - - // # 8 MULTs - __m512d two_b_mul_NEG_TWO_SIGMA_INV_vec = _mm512_mul_pd(two_b_vec, NEG_TWO_SIGMA_INV_VEC); - - // # 160 MULTs, 112 ADDs - // objdump -d --demangle /usr/lib/x86_64-linux-gnu/libsleef.so | grep -A 67 "Sleef_expd8_u10@@Base" - __m512d e_vec = Sleef_expd8_u10(two_b_mul_NEG_TWO_SIGMA_INV_vec); - - // # 16 MULTs - __m512d scalar_mult_common_right_vec = _mm512_mul_pd(NEG_V0_AB_DIV_TWO_SIGMA_VEC, e_vec); - __m512d scalar_mult_common_vec = _mm512_mul_pd(scalar_mult_common_left_vec, scalar_mult_common_right_vec); - - // # 16 MULTs - __m512d scalar_mult_r_ab_vec = _mm512_mul_pd(scalar_mult_common_vec, r_ab2_norm_vec); - __m512d scalar_mult_r_ab2_vec = _mm512_mul_pd(scalar_mult_common_vec, r_ab_norm_vec); - - // # 16 MULTs - __m512d r_ab2_x_mul_scalar_mult_r_ab2_vec = _mm512_mul_pd(r_ab2_x_vec, scalar_mult_r_ab2_vec); - __m512d r_ab2_y_mul_scalar_mult_r_ab2_vec = _mm512_mul_pd(r_ab2_y_vec, scalar_mult_r_ab2_vec); - - // # 16 MULTs, 16 ADDs - __m512d f_alpha_beta_neg_x_vec = _mm512_fmadd_pd(r_ab_x_vec, scalar_mult_r_ab_vec, r_ab2_x_mul_scalar_mult_r_ab2_vec); - __m512d f_alpha_beta_neg_y_vec = _mm512_fmadd_pd(r_ab_y_vec, scalar_mult_r_ab_vec, r_ab2_y_mul_scalar_mult_r_ab2_vec); - - // # 8 MULTs - __m512d rgt_vec = _mm512_mul_pd(alpha_e_y_vec, f_alpha_beta_neg_y_vec); - - // # 8 MULTs, 8 ADDs - __m512d lft_add_rgt_vec = _mm512_fmadd_pd(alpha_e_x_vec, f_alpha_beta_neg_x_vec, rgt_vec); - - // # 8 MULTs - __m512d a_vec = _mm512_mul_pd(lft_add_rgt_vec, cos_phi_inv_vec); - - // # 8 MULTs - __m512d f_alpha_beta_neg_y_sq_vec = _mm512_mul_pd(f_alpha_beta_neg_y_vec, f_alpha_beta_neg_y_vec); - - // # 8 MULTs, 8 ADDs - __m512d b_sq_vec = _mm512_fmadd_pd(f_alpha_beta_neg_x_vec, f_alpha_beta_neg_x_vec, f_alpha_beta_neg_y_sq_vec); - - __m512d w_val_vec = dt_neg_vec; - __mmask8 a_ge_zero_mask = _mm512_cmp_pd_mask(a_vec, ZERO_VEC, _CMP_GE_OQ); - - // # 8 MULTs - __m512d a_sq_vec = _mm512_mul_pd(a_vec, a_vec); - - __mmask8 a_sq_gt_b_sq_mask = _mm512_cmp_pd_mask(a_sq_vec, b_sq_vec, _CMP_GT_OQ); - __mmask8 cond_mask = _kand_mask8(a_ge_zero_mask, a_sq_gt_b_sq_mask); - w_val_vec = _mm512_mask_mov_pd(w_val_vec, cond_mask, c_neg_times_dt_vec); - - // # 16 MULTs - __m512d x_term_vec = _mm512_mul_pd(f_alpha_beta_neg_x_vec, w_val_vec); - __m512d y_term_vec = _mm512_mul_pd(f_alpha_beta_neg_y_vec, w_val_vec); - - // # 16 ADDs - f_alpha_dt_x += _mm512_reduce_add_pd(x_term_vec); - f_alpha_dt_y += _mm512_reduce_add_pd(y_term_vec); - } - } - - // # ############################################################################################################################################################################################ - - for (; j < n_p; j++) { - - if (j == i) { - continue; - } - - // # 2 ADDs - double r_ab_x = p.curr_pos_x - pedestrians[j].curr_pos_x; - double r_ab_y = p.curr_pos_y - pedestrians[j].curr_pos_y; - - // # 1 SQRT, 2 MULTs, 1 ADD - double r_ab_norm = sqrt(r_ab_x * r_ab_x + r_ab_y * r_ab_y); - - // # 2 ADDs - double r_ab2_x = r_ab_x - pedestrians[j].s_b_e_x; - double r_ab2_y = r_ab_y - pedestrians[j].s_b_e_y; - - // # 1 SQRT, 2 MULTs, 1 ADD - double r_ab2_norm = sqrt(r_ab2_x * r_ab2_x + r_ab2_y * r_ab2_y); - - // # 1 ADD - double r_ab_norm_sum = r_ab_norm + r_ab2_norm; - - // # 1 MULT, 1 ADD - double two_b_sq = r_ab_norm_sum * r_ab_norm_sum - pedestrians[j].s_b_sq; - - // # 1 SQRT - double two_b = sqrt(two_b_sq); - - // # 5 MULTs, 1 DIV, 1 EXP - double scalar_mult_common = (r_ab_norm_sum / (two_b * r_ab_norm * r_ab2_norm)) * (NEG_V0_AB_div_TWO_SIGMA * exp(two_b * NEG_TWO_SIGMA_INV)); - - // # 2 MULTs - double scalar_mult_r_ab = scalar_mult_common * r_ab2_norm; - double scalar_mult_r_ab2 = scalar_mult_common * r_ab_norm; - - // # 4 MULTs, 2 ADDs - double f_alpha_beta_neg_x = r_ab_x * scalar_mult_r_ab + r_ab2_x * scalar_mult_r_ab2; - double f_alpha_beta_neg_y = r_ab_y * scalar_mult_r_ab + r_ab2_y * scalar_mult_r_ab2; - - double w_val; - - // # 3 MULTs, 1 ADD - double a = (p.e_x * f_alpha_beta_neg_x + p.e_y * f_alpha_beta_neg_y) * cos_phi_inv; - - // # 2 MULTs, 1 ADD - double b_sq = f_alpha_beta_neg_x * f_alpha_beta_neg_x + f_alpha_beta_neg_y * f_alpha_beta_neg_y; - - w_val = dt_neg; - // # 1 MULT - if (a >= 0 && (a * a > b_sq)) { - w_val = c_neg_times_dt; - } - - // # 2 MULTs, 2 ADDs - f_alpha_dt_x = f_alpha_dt_x + f_alpha_beta_neg_x * w_val; - f_alpha_dt_y = f_alpha_dt_y + f_alpha_beta_neg_y * w_val; - } - - // # ############################################################################################################################################################################################ - - for (int k = 0; k < board.n_obstacles; k++) { - - // # 2 ADDs - double v1_x = p.curr_pos_x - board.obstacles[k].start_x; - double v1_y = p.curr_pos_y - board.obstacles[k].start_y; - - // # 3 MULTs, 1 ADD - double projection = (v1_x * board.obstacles[k].obstacle_v_x + v1_y * board.obstacles[k].obstacle_v_y) * board.obstacles[k].length_sq_inv; - - if (projection < 0.0) - projection = 0.0; - else if (projection > 1.0) - projection = 1.0; - - // # 2 MULTs, 4 ADDs - double r_aB_x = p.curr_pos_x - (board.obstacles[k].start_x + board.obstacles[k].obstacle_v_x * projection); - double r_aB_y = p.curr_pos_y - (board.obstacles[k].start_y + board.obstacles[k].obstacle_v_y * projection); - - // # 1 SQRT, 2 MULTs, 1 ADD - double r_aB_norm = sqrt(r_aB_x * r_aB_x + r_aB_y * r_aB_y); - - // # 2 MULTs, 1 EXP, 1 DIV - double scalar_mult = (NEG_U0_AB_div_R * exp(r_aB_norm * NEG_R_INV)) / r_aB_norm; - - // # 2 MULTs - double f_alpha_B_neg_x = r_aB_x * scalar_mult; - double f_alpha_B_neg_y = r_aB_y * scalar_mult; - - double w_val; - - // # 3 MULTs, 1 ADD - double a = (p.e_x * f_alpha_B_neg_x + p.e_y * f_alpha_B_neg_y) * cos_phi_inv; - - // # 2 MULTs, 1 ADD - double b_sq = f_alpha_B_neg_x * f_alpha_B_neg_x + f_alpha_B_neg_y * f_alpha_B_neg_y; - - w_val = dt_neg; - - // # 1 MULT - if (a >= 0 && (a * a > b_sq)) { - w_val = c_neg_times_dt; - } - - // # 2 MULTs - double f_alpha_B_dt_x = f_alpha_B_neg_x * w_val; - double f_alpha_B_dt_y = f_alpha_B_neg_y * w_val; - - // # 2 ADDs - f_alpha_dt_x += f_alpha_B_dt_x; - f_alpha_dt_y += f_alpha_B_dt_y; - } - - // # 2 ADDs - p.w_x += f_alpha_dt_x; - p.w_y += f_alpha_dt_y; - } - - // # ############################################################################################################################################################################################ - - for (int i = 0; i < n_p; i++) { - Pedestrian &p = pedestrians[i]; - - // # 1 SQRT, 2 MULTs, 1 ADD - double w_norm = sqrt(p.w_x * p.w_x + p.w_y * p.w_y); - - // # 1 DIV - double correction_factor = min(1.0, p.v_max / w_norm); - - // # 2 MULTs - p.curr_vel_x = p.w_x * correction_factor; - p.curr_vel_y = p.w_y * correction_factor; - - // # 2 MULTs, 2 ADDs - p.curr_pos_x = p.curr_pos_x + p.curr_vel_x * dt; - p.curr_pos_y = p.curr_pos_y + p.curr_vel_y * dt; - - if (p.color == 0 && p.curr_pos_x > BOARD_LENGTH_MINUS_ONE_POINT_FIVE) { - p.curr_pos_x = 1.0; - p.destination_y = p.curr_pos_y; - p.w_x = p.desired_v0; - p.w_y = 0.0; - p.curr_vel_x = p.w_x; - p.curr_vel_y = p.w_y; - } - if (p.color == 1 && p.curr_pos_x < 1.5) { - p.curr_pos_x = BOARD_LENGTH_MINUS_ONE; - p.destination_y = p.curr_pos_y; - p.w_x = -p.desired_v0; - p.w_y = 0.0; - p.curr_vel_x = p.w_x; - p.curr_vel_y = p.w_y; - } - - p.curr_pos_y = clamp(p.curr_pos_y, WALL_EPS, BOARD_WIDTH_MINUS_WALL_EPS); - - // # 2 ADDs - double e_x = p.destination_x - p.curr_pos_x; - double e_y = p.destination_y - p.curr_pos_y; - - // # 1 DIV, 1 SQRT, 2 MULTs, 1 ADD - double e_norm_inv = 1.0 / sqrt(e_x * e_x + e_y * e_y); - - // # 2 MULTs - p.e_x = e_x * e_norm_inv; - p.e_y = e_y * e_norm_inv; - - // # 2 MULTs, 1 ADD - double velocity_norm_sq = p.curr_vel_x * p.curr_vel_x + p.curr_vel_y * p.curr_vel_y; - - // # 1 SQRT - double velocity_norm = sqrt(velocity_norm_sq); - - // # 5 MULTs - p.s_b_sq = velocity_norm_sq * delta_t_sq; - p.s_b_e_x = p.e_x * (velocity_norm * delta_t); - p.s_b_e_y = p.e_y * (velocity_norm * delta_t); - -#if !BENCHMARK - if (dump_path && (t % dump_interval == 0)) { - samples.emplace_back(t + 1, i, p.curr_pos_x, p.curr_pos_y, p.e_x, p.e_y, p.color); - } -#endif - } - } - -#if !BENCHMARK - if (dump_path) { - ofstream out(dump_path); - out << "time,id,x,y,p.e.x,p.e.y,color\n"; - - for (auto &[t, id, x, y, e_x, e_y, c] : samples) { - out << t << ',' << id << ',' << x << ',' << y << ',' << e_x << ',' << e_y << ',' << c << '\n'; - } - } -#endif - - return 0; -} - -#endif diff --git a/variants/variant_baseline_flopcount.cpp b/variants/variant_baseline_flopcount.cpp deleted file mode 100644 index 6f4f90e..0000000 --- a/variants/variant_baseline_flopcount.cpp +++ /dev/null @@ -1,440 +0,0 @@ -// This file is merged from the forces.cpp and social_forces_baseline.cpp -// The goal is that we from now on only have one file as our baseline -// This should make it more easy in the future to implement optimizations -// Our naming convention from now on is: variant_.cpp (Please make sure to keep that consistent, otherwise our Makefile wont work anymore) - -#include "../headers/board.hpp" -#include "../headers/constants.hpp" -#include "../headers/pedestrian.hpp" -#include "../headers/vector2d.hpp" - -#include -#include -#include -#include -#include // keep — forces code uses unqualified sqrt/exp -#include -#include -#include - -#define SEED 42 -#define BOARD_WIDTH 10.0 -#define BOARD_LENGTH 50.0 -#define WALL_EPS 1e-3 -#define BOARD_WIDTH_MINUS_WALL_EPS 9.999 - -#ifndef RANDOM_SEED -#define RANDOM_SEED 0 -#endif - -#ifndef BENCHMARK -#define BENCHMARK 0 -#endif - -using namespace std; - -// forces.cpp - -// Flop count: 2 ADDs, 6 MULTs -double w_ped(Vector2d e, Vector2d f) { - double f_norm_sq = f.dot(f); - double rhs = f_norm_sq * cos_phi_sq; - double e_dot_f = e.dot(f); - double e_dot_f_sq = e_dot_f * e_dot_f; - if (e_dot_f >= 0.0 || e_dot_f_sq <= rhs) - return dt_neg; - else - return c_neg_times_dt; -} - -double w_obstacle(Vector2d e, Vector2d f) { - double e_dot_f = (f.y == 0) ? 0 : (f.y < 0 ? -e.y : e.y); - if (e_dot_f >= cos_phi) - return dt_neg; - else - return c_neg_times_dt; -} - -// Flop count: 2 ADDs, 4 MULTs -Vector2d acceleration_term_dt(const Pedestrian &alpha) { - - // # 2 MULTs - Vector2d desired_velocity = alpha.e * alpha.desired_v0; - - Vector2d actual_velocity = alpha.curr_vel; - - // # 2 ADDs, 2 MULTs - Vector2d f_alpha_0_dt = (desired_velocity - actual_velocity) * dt_times_tau_alpha_inv; - - return f_alpha_0_dt; -} - -// Flop count: 12 ADDs, 24 MULTs, 1 DIV, 3 SQRTs, 1 EXP -Vector2d pedestrian_repulsion_dt(const Pedestrian &alpha, const Pedestrian &beta) { - - // # 2 ADDs - // $\vec{r}_{\alpha\beta}$ - Vector2d r_ab = alpha.curr_pos - beta.curr_pos; - - // # 1 ADD, 2 MULTs, 1 SQRT - // $\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert$ - double r_ab_norm = r_ab.norm(); - r_ab_norm = r_ab_norm > EPS ? r_ab_norm : EPS; // Avoid division by zero - - // $v_\beta \Delta t \vec{e}_\beta$ - Vector2d s_b_e = beta.s_b_e; - - // # 2 ADDs - // $\vec{r}_{\alpha\beta} - \vec{v}_{\beta}\varDelta_{t}\vec{e}_{\beta}$ - Vector2d r_ab2 = r_ab - s_b_e; - - // # 1 ADD, 2 MULTs, 1 SQRT - // $\Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert$ - double r_ab2_norm = r_ab2.norm(); - r_ab2_norm = r_ab2_norm > EPS ? r_ab2_norm : EPS; // Avoid division by zero - - // # 1 ADD - // $\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert$ - double r_ab_norm_sum = r_ab_norm + r_ab2_norm; - - // # 1 ADD, 1 MULTs, 1 SQRT - // - // $2b = \sqrt{(\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert)^2 - (v_\beta \Delta t)^2}$ - // - double two_b_sq = r_ab_norm_sum * r_ab_norm_sum - beta.s_b_sq; - two_b_sq = two_b_sq > EPS_SQ ? two_b_sq : EPS_SQ; // Avoid sqrt of negative and avoid division by zero - double two_b = sqrt(two_b_sq); - - // # 5 MULTs, 1 DIV, 1 EXP - // - // $\frac{-V^0_{\alpha\beta}}{2\sigma}\,e^{-b/\sigma}\,\frac{\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}{2b\,\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert\,\Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}$ - // - double scalar_mult_common = (r_ab_norm_sum / (two_b * r_ab_norm * r_ab2_norm)) * (NEG_V0_AB_div_TWO_SIGMA * exp(two_b * NEG_TWO_SIGMA_INV)); - - // # 1 MULT - // - // $\frac{-V^0_{\alpha\beta}}{2\sigma}\,e^{-b/\sigma}\,\frac{\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}{2b\,\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert}$ - // - double scalar_mult_r_ab = scalar_mult_common * r_ab2_norm; - - // # 1 MULT - // - // $\frac{-V^0_{\alpha\beta}}{2\sigma}\,e^{-b/\sigma}\,\frac{\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}{2b\,\Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}$ - // - double scalar_mult_r_ab2 = scalar_mult_common * r_ab_norm; - - // # 2 ADDs, 4 MULTs - // - // $\vec{-f}_{\alpha \beta}(\vec{r}_{\alpha \beta}) = \frac{-V^0_{\alpha \beta}}{2\sigma}e^{-b / \sigma} 2\grad_{\vec{r}_{\alpha \beta}}b $ - // - // - // $= \frac{-V^0_{\alpha\beta}}{2\sigma}\,e^{-b/\sigma}\,\frac{\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}{2b\,\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert}\,\vec{r}_{\alpha \beta} + \frac{-V^0_{\alpha\beta}}{2\sigma}\,e^{-b/\sigma}\,\frac{\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}{2b\,\Arrowvert \vec{r}_{\alpha \beta} - - // v_\beta \Delta t \vec{e}_\beta \Arrowvert}\,(\vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta)$ - // - // - // - // $2\grad_{\vec{r}_{\alpha \beta}}b = \frac{(\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert)(\frac{\vec{r}_{\alpha \beta}}{\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert} + \frac{\vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta}{\Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert})}{2b}$ - // - // - // $= \frac{\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta Delta t \vec{e}_\beta \Arrowvert}{2b\,\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert}\,\vec{r}_{\alpha \beta} + \frac{\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}{2b\,\Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}\,(\vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta)$ - // - // - Vector2d f_alpha_beta_neg = r_ab * scalar_mult_r_ab + r_ab2 * scalar_mult_r_ab2; - - // # 2 ADDs, 8 MULTs - // - // $\vec{F}_{\alpha\beta}(\vec{e}_{\alpha}, \vec{r}_{\alpha \beta}).dt = w(\vec{e}_{\alpha}, -\vec{f}_{\alpha\beta}(\vec{r}_{\alpha \beta})).\vec{f}_{\alpha\beta}(\vec{r}_{\alpha \beta}).dt$ - // - Vector2d f_alpha_beta_dt = f_alpha_beta_neg * w_ped(alpha.e, f_alpha_beta_neg); - - return f_alpha_beta_dt; -} - -// Flop count: 8 ADDs, 13 MULTs, 1 DIV, 1 SQRT, 1 EXP -Vector2d obstacle_repulsion_dt(const Pedestrian &alpha, const Obstacle &B) { - - // # 7 ADDs, 5 MULTs - // $\vec{r}_{\alpha B}$ - Vector2d r_aB = alpha.curr_pos - B.find_nearest_obstacle_point(alpha.curr_pos); - - // # 1 ADD, 2 MULTs, 1 SQRT - // $\Arrowvert \vec{r}_{\alpha B} \Arrowvert$ - double r_aB_norm = r_aB.norm(); - r_aB_norm = r_aB_norm > EPS ? r_aB_norm : EPS; // Avoid division by zero - - // # 2 MULTs, 1 DIV, 1 EXP - // - // $\frac{-U^0_{\alpha B}}{R} e^{- \Arrowvert \vec{r}_{\alpha B} \Arrowvert / R}$ - // - double scalar_mult = (NEG_U0_AB_div_R * exp(r_aB_norm * NEG_R_INV)) / r_aB_norm; - - // # 2 MULTs - // - // $-\vec{F}_{\alpha B}(\vec{r}_{\alpha B}) = \frac{-U^0_{\alpha B}}{R} e^{- \Arrowvert \vec{r}_{\alpha B} \Arrowvert / R} \frac{\vec{r}_{\alpha B}}{\Arrowvert \vec{r}_{\alpha B} \Arrowvert}$ - // - Vector2d f_alpha_B_neg = r_aB * scalar_mult; - - // # 2 MULTs - // - // $\vec{F}_{\alpha B}(\vec{e}_{\alpha}, \vec{r}_{\alpha B}).dt = w(\vec{e}_{\alpha}, -\vec{F}_{\alpha B}(\vec{r}_{\alpha B})).\vec{F}_{\alpha B}(\vec{r}_{\alpha B}).dt$ - // - Vector2d f_alpha_B_dt = f_alpha_B_neg * w_obstacle(alpha.e, f_alpha_B_neg); - - return f_alpha_B_dt; -} - -// social_force_baseline.cpp - -Pedestrian *init_pedestrians(int n_p) { -#if RANDOM_SEED - random_device rd; // Platform-provided random seed source. - mt19937 gen(rd()); // Non-reproducible initialization. -#else - mt19937 gen(SEED); // Reproducible initialization. -#endif - - normal_distribution velocity_gen(v0_mean, v0_std_dev); - uniform_real_distribution y_gen(0.3, BOARD_WIDTH - 0.3); - uniform_real_distribution x_gen_left(2.0, BOARD_LENGTH / 2.0); - uniform_real_distribution x_gen_right(BOARD_LENGTH / 2.0, BOARD_LENGTH - 2.0); - - Pedestrian *P = new Pedestrian[n_p]; - - int n_p_half = n_p / 2; - - // Group of pedestrians walking from left to right. - for (int i = 0; i < n_p_half; i++) { - double x = x_gen_left(gen); - double y = y_gen(gen); - - Vector2d start_pos = Vector2d(x, y); - Vector2d destination = Vector2d(BOARD_LENGTH + 30.0, y); - - Vector2d init_e = destination - start_pos; - double init_e_norm = init_e.norm(); - Vector2d init_e_unit = init_e / init_e_norm; - - // Reject near-zero desired speeds to avoid pedestrians that barely move, - // and reject negative speeds because they would reverse the intended - // walking direction. - double desired_v0; - do { - desired_v0 = velocity_gen(gen); - } while (desired_v0 <= 0.1); - - double v_max = 1.3 * desired_v0; - - P[i] = Pedestrian(WHITE, Vector2d(), start_pos, init_e_unit, v_max, desired_v0, destination); - } - - // Group of pedestrians walking from right to left. - for (int i = n_p_half; i < n_p; i++) { - double x = x_gen_right(gen); - double y = y_gen(gen); - - Vector2d start_pos = Vector2d(x, y); - Vector2d destination = Vector2d(-30.0, y); - - Vector2d init_e = destination - start_pos; - double init_e_norm = init_e.norm(); - Vector2d init_e_unit = init_e / init_e_norm; - - // Reject near-zero desired speeds to avoid pedestrians that barely move, - // and reject negative speeds because they would reverse the intended - // walking direction. - double desired_v0; - do { - desired_v0 = velocity_gen(gen); - } while (desired_v0 <= 0.1); - - double v_max = 1.3 * desired_v0; - - P[i] = Pedestrian(BLACK, Vector2d(), start_pos, init_e_unit, v_max, desired_v0, destination); - } - - return P; -} - -Board init_board_without_obstacles() { return Board(BOARD_LENGTH, BOARD_WIDTH, 0, nullptr); } - -int run_simulation(Pedestrian *pedestrians, Board &board, int n_p, int n_iter, const char *dump_path) { - -#if BENCHMARK - // Avoid compiler unused parameter warning - (void)dump_path; -#endif - -#if !BENCHMARK - const int dump_interval = 1; - - // Each sample stores: time step, pedestrian id, x, y, e.x, e.y, color. - vector> samples; - - if (dump_path) { - samples.reserve(static_cast(n_iter / dump_interval + 1) * n_p); - - // Store the initial state before the first simulation step. - for (int i = 0; i < n_p; i++) { - Pedestrian &p = pedestrians[i]; - samples.emplace_back(0, i, p.curr_pos.x, p.curr_pos.y, p.e.x, p.e.y, p.color); - } - } -#endif - - // # ADDs: n_iter * ( - // # n_p * 4 + - // # n_p * (n_p - 1) * 14 + - // # n_p * board.n_obstacles * 10 + - // # n_p * 7 - // # ) - - // # MULTs: n_iter * ( - // # n_p * 4 - // # n_p * (n_p - 1) * 24 + - // # n_p * board.n_obstacles * 13 + - // # n_p * 16 - // # ) - - // # DIVs: n_iter * ( - // # n_p * (n_p - 1) + - // # n_p * board.n_obstacles + - // # n_p * 2 - // # ) - - // # SQRTs: n_iter * ( - // # n_p * (n_p - 1) * 3 + - // # n_p * board.n_obstacles + - // # n_p * 3 - // # ) - - // # EXPs: n_iter * ( - // # n_p * (n_p - 1) + - // # n_p * board.n_obstacles - // # ) - - for (int t = 0; t < n_iter; t++) { - // Pass 1: Compute all forces from the current state. - for (int i = 0; i < n_p; i++) { - Pedestrian &p = pedestrians[i]; - - // # 2 ADDs, 4 MULTs - // - // $\vec{F}_\alpha(\vec{v}_\alpha, v^0_\alpha \vec{e}_\alpha) = \frac{1}{\tau_\alpha}(v^0_\alpha \vec{e}_\alpha - \vec{v}_\alpha)$ - // - Vector2d f_alpha_dt = acceleration_term_dt(p); - - // # (N-1) * (14 ADDs, 24 MULTs, 1 DIV, 3 SQRTs, 1 EXP) - // - // $\sum_{\beta} \vec{F}_{\alpha\beta}(\vec{e}_{\alpha},\vec{r}_{\alpha} - \vec{r}_{\beta})$ - // - for (int j = 0; j < n_p; j++) { - if (i == j) - continue; - - f_alpha_dt += pedestrian_repulsion_dt(p, pedestrians[j]); - } - - // # board.n_obstacles * (10 ADDs, 13 MULTs, 1 DIV, 1 SQRT, 1 EXP) - // - // $\sum_{B} \vec{F}_{\alpha B}(\vec{e}_{\alpha}, \vec{r}_{\alpha} - \vec{r}_{B})$ - // - for (int k = 0; k < board.n_obstacles; k++) { - f_alpha_dt += obstacle_repulsion_dt(p, board.obstacles[k]); - } - - // # 2 ADDs - // - // - // $\frac{d\vec{w}_\alpha}{dt} = \vec{F}_{\alpha}(t)$ - // - p.w += f_alpha_dt; - } - - // Pass 2: Update pedestrian states after all forces have been computed. - // # (7 ADDs, 16 MULTs, 2 DIVs, 3 SQRTs) * n_p - for (int i = 0; i < n_p; i++) { - Pedestrian &p = pedestrians[i]; - - // # 1 ADD, 2 MULTs, 1 SQRT - double w_norm = p.w.norm(); - - double w_norm2 = w_norm > EPS ? w_norm : EPS; // Avoid division by zero. - - // # 1 DIV - double correction_factor = std::min(1.0, p.v_max / w_norm2); - - // # 2 MULTs - p.curr_vel = p.w * correction_factor; - - // # 2 ADDs, 2 MULTs - // - // $\vec{r}_i^{t+1} = \vec{r}_i^t + \vec{v}_i^t dt$ - // - p.curr_pos = p.curr_pos + p.curr_vel * dt; - - // When a pedestrian reaches the end of the corridor, teleport them back - // to the opposite entrance while keeping the same y-coordinate. This - // maintains a constant pedestrian flow and makes lane formation easier to - // observe. - if (p.color == WHITE && p.curr_pos.x > (BOARD_LENGTH - 1.5)) { - p.curr_pos.x = 1.0; - p.destination.y = p.curr_pos.y; - p.w = Vector2d(p.desired_v0, 0.0); - p.curr_vel = p.w; - } - if (p.color == BLACK && p.curr_pos.x < 1.5) { - p.curr_pos.x = BOARD_LENGTH - 1.0; - p.destination.y = p.curr_pos.y; - p.w = Vector2d(-p.desired_v0, 0.0); - p.curr_vel = p.w; - } - - // Keep pedestrians inside the board boundaries. - p.curr_pos.y = clamp(p.curr_pos.y, WALL_EPS, BOARD_WIDTH_MINUS_WALL_EPS); - - // # 2 ADDs - Vector2d e = p.destination - p.curr_pos; - - // # 1 ADD, 2 MULTs, 1 DIV, 1 SQRT - double e_norm = e.norm(); - double e_norm2 = e_norm > EPS ? e_norm : EPS; // Avoid division by zero. - double e_norm_inv = 1.0 / e_norm2; - - // # 2 MULTs - p.e = e * e_norm_inv; - - // # 1 ADD, 2 MULTs - double velocity_norm_sq = p.curr_vel.dot(p.curr_vel); - - // # 1 SQRT - double velocity_norm = std::sqrt(velocity_norm_sq); - - // # 1 MULT - p.s_b_sq = velocity_norm_sq * delta_t_sq; - - // # 3 MULTs - p.s_b_e = p.e * (velocity_norm * delta_t); - -#if !BENCHMARK - // Store the current pedestrian state for visualization. - if (dump_path && (t % dump_interval == 0)) { - samples.emplace_back(t + 1, i, p.curr_pos.x, p.curr_pos.y, p.e.x, p.e.y, p.color); - } -#endif - } - } - -#if !BENCHMARK - // Export recorded simulation states to CSV. - if (dump_path) { - ofstream out(dump_path); - out << "time,id,x,y,p.e.x,p.e.y,color\n"; - - for (auto &[t, id, x, y, e_x, e_y, c] : samples) { - out << t << ',' << id << ',' << x << ',' << y << ',' << e_x << ',' << e_y << ',' << c << '\n'; - } - } -#endif - - return 0; -} \ No newline at end of file diff --git a/variants/variant_float_baseline_flopcount.cpp b/variants/variant_float_baseline_flopcount.cpp deleted file mode 100644 index 820f769..0000000 --- a/variants/variant_float_baseline_flopcount.cpp +++ /dev/null @@ -1,438 +0,0 @@ -// # Changes: Added per-instruction flop count -// # Predecessor: variant_float_baseline.cpp - -#include "../headers/board_float.hpp" -#include "../headers/constants_float.hpp" -#include "../headers/pedestrian_float.hpp" -#include "../headers/vector2d_float.hpp" - -#include -#include -#include -#include -#include // keep — forces code uses unqualified sqrt/exp -#include -#include -#include - -#define SEED 42 -#define BOARD_WIDTH 10.0f -#define BOARD_LENGTH 50.0f -#define WALL_EPS 1e-3f -#define BOARD_WIDTH_MINUS_WALL_EPS 9.999f - -#ifndef RANDOM_SEED -#define RANDOM_SEED 0 -#endif - -#ifndef BENCHMARK -#define BENCHMARK 0 -#endif - -using namespace std; - -// forces.cpp - -// Flop count: 2 ADDs, 6 MULTs -float w_ped(Vector2d e, Vector2d f) { - float f_norm_sq = f.dot(f); - float rhs = f_norm_sq * cos_phi_sq; - float e_dot_f = e.dot(f); - float e_dot_f_sq = e_dot_f * e_dot_f; - if (e_dot_f >= 0.0 || e_dot_f_sq <= rhs) - return dt_neg; - else - return c_neg_times_dt; -} - -float w_obstacle(Vector2d e, Vector2d f) { - float e_dot_f = (f.y == 0) ? 0 : (f.y < 0 ? -e.y : e.y); - if (e_dot_f >= cos_phi) - return dt_neg; - else - return c_neg_times_dt; -} - -// Flop count: 2 ADDs, 4 MULTs -Vector2d acceleration_term_dt(const Pedestrian &alpha) { - - // # 2 MULTs - Vector2d desired_velocity = alpha.e * alpha.desired_v0; - - Vector2d actual_velocity = alpha.curr_vel; - - // # 2 ADDs, 2 MULTs - Vector2d f_alpha_0_dt = (desired_velocity - actual_velocity) * dt_times_tau_alpha_inv; - - return f_alpha_0_dt; -} - -// Flop count: 12 ADDs, 24 MULTs, 1 DIV, 3 SQRTs, 1 EXP -Vector2d pedestrian_repulsion_dt(const Pedestrian &alpha, const Pedestrian &beta) { - - // # 2 ADDs - // $\vec{r}_{\alpha\beta}$ - Vector2d r_ab = alpha.curr_pos - beta.curr_pos; - - // # 1 ADD, 2 MULTs, 1 SQRT - // $\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert$ - float r_ab_norm = r_ab.norm(); - r_ab_norm = r_ab_norm > EPS ? r_ab_norm : EPS; // Avoid division by zero - - // $v_\beta \Delta t \vec{e}_\beta$ - Vector2d s_b_e = beta.s_b_e; - - // # 2 ADDs - // $\vec{r}_{\alpha\beta} - \vec{v}_{\beta}\varDelta_{t}\vec{e}_{\beta}$ - Vector2d r_ab2 = r_ab - s_b_e; - - // # 1 ADD, 2 MULTs, 1 SQRT - // $\Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert$ - float r_ab2_norm = r_ab2.norm(); - r_ab2_norm = r_ab2_norm > EPS ? r_ab2_norm : EPS; // Avoid division by zero - - // # 1 ADD - // $\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert$ - float r_ab_norm_sum = r_ab_norm + r_ab2_norm; - - // # 1 ADD, 1 MULTs, 1 SQRT - // - // $2b = \sqrt{(\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert)^2 - (v_\beta \Delta t)^2}$ - // - float two_b_sq = r_ab_norm_sum * r_ab_norm_sum - beta.s_b_sq; - two_b_sq = two_b_sq > EPS_SQ ? two_b_sq : EPS_SQ; // Avoid sqrt of negative and avoid division by zero - float two_b = sqrt(two_b_sq); - - // # 5 MULTs, 1 DIV, 1 EXP - // - // $\frac{-V^0_{\alpha\beta}}{2\sigma}\,e^{-b/\sigma}\,\frac{\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}{2b\,\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert\,\Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}$ - // - float scalar_mult_common = (r_ab_norm_sum / (two_b * r_ab_norm * r_ab2_norm)) * (NEG_V0_AB_div_TWO_SIGMA * exp(two_b * NEG_TWO_SIGMA_INV)); - - // # 1 MULT - // - // $\frac{-V^0_{\alpha\beta}}{2\sigma}\,e^{-b/\sigma}\,\frac{\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}{2b\,\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert}$ - // - float scalar_mult_r_ab = scalar_mult_common * r_ab2_norm; - - // # 1 MULT - // - // $\frac{-V^0_{\alpha\beta}}{2\sigma}\,e^{-b/\sigma}\,\frac{\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}{2b\,\Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}$ - // - float scalar_mult_r_ab2 = scalar_mult_common * r_ab_norm; - - // # 2 ADDs, 4 MULTs - // - // $\vec{-f}_{\alpha \beta}(\vec{r}_{\alpha \beta}) = \frac{-V^0_{\alpha \beta}}{2\sigma}e^{-b / \sigma} 2\grad_{\vec{r}_{\alpha \beta}}b $ - // - // - // $= \frac{-V^0_{\alpha\beta}}{2\sigma}\,e^{-b/\sigma}\,\frac{\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}{2b\,\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert}\,\vec{r}_{\alpha \beta} + \frac{-V^0_{\alpha\beta}}{2\sigma}\,e^{-b/\sigma}\,\frac{\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}{2b\,\Arrowvert \vec{r}_{\alpha \beta} - - // v_\beta \Delta t \vec{e}_\beta \Arrowvert}\,(\vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta)$ - // - // - // - // $2\grad_{\vec{r}_{\alpha \beta}}b = \frac{(\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert)(\frac{\vec{r}_{\alpha \beta}}{\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert} + \frac{\vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta}{\Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert})}{2b}$ - // - // - // $= \frac{\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta Delta t \vec{e}_\beta \Arrowvert}{2b\,\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert}\,\vec{r}_{\alpha \beta} + \frac{\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}{2b\,\Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}\,(\vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta)$ - // - // - Vector2d f_alpha_beta_neg = r_ab * scalar_mult_r_ab + r_ab2 * scalar_mult_r_ab2; - - // # 2 ADDs, 8 MULTs - // - // $\vec{F}_{\alpha\beta}(\vec{e}_{\alpha}, \vec{r}_{\alpha \beta}).dt = w(\vec{e}_{\alpha}, -\vec{f}_{\alpha\beta}(\vec{r}_{\alpha \beta})).\vec{f}_{\alpha\beta}(\vec{r}_{\alpha \beta}).dt$ - // - Vector2d f_alpha_beta_dt = f_alpha_beta_neg * w_ped(alpha.e, f_alpha_beta_neg); - - return f_alpha_beta_dt; -} - -// Flop count: 8 ADDs, 13 MULTs, 1 DIV, 1 SQRT, 1 EXP -Vector2d obstacle_repulsion_dt(const Pedestrian &alpha, const Obstacle &B) { - - // # 7 ADDs, 5 MULTs - // $\vec{r}_{\alpha B}$ - Vector2d r_aB = alpha.curr_pos - B.find_nearest_obstacle_point(alpha.curr_pos); - - // # 1 ADD, 2 MULTs, 1 SQRT - // $\Arrowvert \vec{r}_{\alpha B} \Arrowvert$ - float r_aB_norm = r_aB.norm(); - r_aB_norm = r_aB_norm > EPS ? r_aB_norm : EPS; // Avoid division by zero - - // # 2 MULTs, 1 DIV, 1 EXP - // - // $\frac{-U^0_{\alpha B}}{R} e^{- \Arrowvert \vec{r}_{\alpha B} \Arrowvert / R}$ - // - float scalar_mult = (NEG_U0_AB_div_R * exp(r_aB_norm * NEG_R_INV)) / r_aB_norm; - - // # 2 MULTs - // - // $-\vec{F}_{\alpha B}(\vec{r}_{\alpha B}) = \frac{-U^0_{\alpha B}}{R} e^{- \Arrowvert \vec{r}_{\alpha B} \Arrowvert / R} \frac{\vec{r}_{\alpha B}}{\Arrowvert \vec{r}_{\alpha B} \Arrowvert}$ - // - Vector2d f_alpha_B_neg = r_aB * scalar_mult; - - // # 2 MULTs - // - // $\vec{F}_{\alpha B}(\vec{e}_{\alpha}, \vec{r}_{\alpha B}).dt = w(\vec{e}_{\alpha}, -\vec{F}_{\alpha B}(\vec{r}_{\alpha B})).\vec{F}_{\alpha B}(\vec{r}_{\alpha B}).dt$ - // - Vector2d f_alpha_B_dt = f_alpha_B_neg * w_obstacle(alpha.e, f_alpha_B_neg); - - return f_alpha_B_dt; -} - -// social_force_baseline.cpp - -Pedestrian *init_pedestrians(int n_p) { -#if RANDOM_SEED - random_device rd; // Platform-provided random seed source. - mt19937 gen(rd()); // Non-reproducible initialization. -#else - mt19937 gen(SEED); // Reproducible initialization. -#endif - - normal_distribution velocity_gen(v0_mean, v0_std_dev); - uniform_real_distribution y_gen(0.3, BOARD_WIDTH - 0.3); - uniform_real_distribution x_gen_left(2.0, BOARD_LENGTH / 2.0); - uniform_real_distribution x_gen_right(BOARD_LENGTH / 2.0, BOARD_LENGTH - 2.0); - - Pedestrian *P = new Pedestrian[n_p]; - - int n_p_half = n_p / 2; - - // Group of pedestrians walking from left to right. - for (int i = 0; i < n_p_half; i++) { - float x = x_gen_left(gen); - float y = y_gen(gen); - - Vector2d start_pos = Vector2d(x, y); - Vector2d destination = Vector2d(BOARD_LENGTH + 30.0, y); - - Vector2d init_e = destination - start_pos; - float init_e_norm = init_e.norm(); - Vector2d init_e_unit = init_e / init_e_norm; - - // Reject near-zero desired speeds to avoid pedestrians that barely move, - // and reject negative speeds because they would reverse the intended - // walking direction. - float desired_v0; - do { - desired_v0 = velocity_gen(gen); - } while (desired_v0 <= 0.1); - - float v_max = 1.3 * desired_v0; - - P[i] = Pedestrian(WHITE, Vector2d(), start_pos, init_e_unit, v_max, desired_v0, destination); - } - - // Group of pedestrians walking from right to left. - for (int i = n_p_half; i < n_p; i++) { - float x = x_gen_right(gen); - float y = y_gen(gen); - - Vector2d start_pos = Vector2d(x, y); - Vector2d destination = Vector2d(-30.0, y); - - Vector2d init_e = destination - start_pos; - float init_e_norm = init_e.norm(); - Vector2d init_e_unit = init_e / init_e_norm; - - // Reject near-zero desired speeds to avoid pedestrians that barely move, - // and reject negative speeds because they would reverse the intended - // walking direction. - float desired_v0; - do { - desired_v0 = velocity_gen(gen); - } while (desired_v0 <= 0.1); - - float v_max = 1.3 * desired_v0; - - P[i] = Pedestrian(BLACK, Vector2d(), start_pos, init_e_unit, v_max, desired_v0, destination); - } - - return P; -} - -Board init_board_without_obstacles() { return Board(BOARD_LENGTH, BOARD_WIDTH, 0, nullptr); } - -int run_simulation(Pedestrian *pedestrians, Board &board, int n_p, int n_iter, const char *dump_path) { - -#if BENCHMARK - // Avoid compiler unused parameter warning - (void)dump_path; -#endif - -#if !BENCHMARK - const int dump_interval = 1; - - // Each sample stores: time step, pedestrian id, x, y, e.x, e.y, color. - vector> samples; - - if (dump_path) { - samples.reserve(static_cast(n_iter / dump_interval + 1) * n_p); - - // Store the initial state before the first simulation step. - for (int i = 0; i < n_p; i++) { - Pedestrian &p = pedestrians[i]; - samples.emplace_back(0, i, p.curr_pos.x, p.curr_pos.y, p.e.x, p.e.y, p.color); - } - } -#endif - - // # ADDs: n_iter * ( - // # n_p * 4 + - // # n_p * (n_p - 1) * 14 + - // # n_p * board.n_obstacles * 10 + - // # n_p * 7 - // # ) - - // # MULTs: n_iter * ( - // # n_p * 4 - // # n_p * (n_p - 1) * 24 + - // # n_p * board.n_obstacles * 13 + - // # n_p * 16 - // # ) - - // # DIVs: n_iter * ( - // # n_p * (n_p - 1) + - // # n_p * board.n_obstacles + - // # n_p * 2 - // # ) - - // # SQRTs: n_iter * ( - // # n_p * (n_p - 1) * 3 + - // # n_p * board.n_obstacles + - // # n_p * 3 - // # ) - - // # EXPs: n_iter * ( - // # n_p * (n_p - 1) + - // # n_p * board.n_obstacles - // # ) - - for (int t = 0; t < n_iter; t++) { - // Pass 1: Compute all forces from the current state. - for (int i = 0; i < n_p; i++) { - Pedestrian &p = pedestrians[i]; - - // # 2 ADDs, 4 MULTs - // - // $\vec{F}_\alpha(\vec{v}_\alpha, v^0_\alpha \vec{e}_\alpha) = \frac{1}{\tau_\alpha}(v^0_\alpha \vec{e}_\alpha - \vec{v}_\alpha)$ - // - Vector2d f_alpha_dt = acceleration_term_dt(p); - - // # (N-1) * (14 ADDs, 24 MULTs, 1 DIV, 3 SQRTs, 1 EXP) - // - // $\sum_{\beta} \vec{F}_{\alpha\beta}(\vec{e}_{\alpha},\vec{r}_{\alpha} - \vec{r}_{\beta})$ - // - for (int j = 0; j < n_p; j++) { - if (i == j) - continue; - - f_alpha_dt += pedestrian_repulsion_dt(p, pedestrians[j]); - } - - // # board.n_obstacles * (10 ADDs, 13 MULTs, 1 DIV, 1 SQRT, 1 EXP) - // - // $\sum_{B} \vec{F}_{\alpha B}(\vec{e}_{\alpha}, \vec{r}_{\alpha} - \vec{r}_{B})$ - // - for (int k = 0; k < board.n_obstacles; k++) { - f_alpha_dt += obstacle_repulsion_dt(p, board.obstacles[k]); - } - - // # 2 ADDs - // - // - // $\frac{d\vec{w}_\alpha}{dt} = \vec{F}_{\alpha}(t)$ - // - p.w += f_alpha_dt; - } - - // Pass 2: Update pedestrian states after all forces have been computed. - // # (7 ADDs, 16 MULTs, 2 DIVs, 3 SQRTs) * n_p - for (int i = 0; i < n_p; i++) { - Pedestrian &p = pedestrians[i]; - - // # 1 ADD, 2 MULTs, 1 SQRT - float w_norm = p.w.norm(); - - float w_norm2 = w_norm > EPS ? w_norm : EPS; // Avoid division by zero. - - // # 1 DIV - float correction_factor = std::min(1.0f, p.v_max / w_norm2); - - // # 2 MULTs - p.curr_vel = p.w * correction_factor; - - // # 2 ADDs, 2 MULTs - // - // $\vec{r}_i^{t+1} = \vec{r}_i^t + \vec{v}_i^t dt$ - // - p.curr_pos = p.curr_pos + p.curr_vel * dt; - - // When a pedestrian reaches the end of the corridor, teleport them back - // to the opposite entrance while keeping the same y-coordinate. This - // maintains a constant pedestrian flow and makes lane formation easier to - // observe. - if (p.color == WHITE && p.curr_pos.x > (BOARD_LENGTH - 1.5)) { - p.curr_pos.x = 1.0; - p.destination.y = p.curr_pos.y; - p.w = Vector2d(p.desired_v0, 0.0); - p.curr_vel = p.w; - } - if (p.color == BLACK && p.curr_pos.x < 1.5) { - p.curr_pos.x = BOARD_LENGTH - 1.0; - p.destination.y = p.curr_pos.y; - p.w = Vector2d(-p.desired_v0, 0.0); - p.curr_vel = p.w; - } - - // Keep pedestrians inside the board boundaries. - p.curr_pos.y = clamp(p.curr_pos.y, WALL_EPS, BOARD_WIDTH_MINUS_WALL_EPS); - - // # 2 ADDs - Vector2d e = p.destination - p.curr_pos; - - // # 1 ADD, 2 MULTs, 1 DIV, 1 SQRT - float e_norm = e.norm(); - float e_norm2 = e_norm > EPS ? e_norm : EPS; // Avoid division by zero. - float e_norm_inv = 1.0 / e_norm2; - - // # 2 MULTs - p.e = e * e_norm_inv; - - // # 1 ADD, 2 MULTs - float velocity_norm_sq = p.curr_vel.dot(p.curr_vel); - - // # 1 SQRT - float velocity_norm = std::sqrt(velocity_norm_sq); - - // # 1 MULT - p.s_b_sq = velocity_norm_sq * delta_t_sq; - - // # 3 MULTs - p.s_b_e = p.e * (velocity_norm * delta_t); - -#if !BENCHMARK - // Store the current pedestrian state for visualization. - if (dump_path && (t % dump_interval == 0)) { - samples.emplace_back(t + 1, i, p.curr_pos.x, p.curr_pos.y, p.e.x, p.e.y, p.color); - } -#endif - } - } - -#if !BENCHMARK - // Export recorded simulation states to CSV. - if (dump_path) { - ofstream out(dump_path); - out << "time,id,x,y,p.e.x,p.e.y,color\n"; - - for (auto &[t, id, x, y, e_x, e_y, c] : samples) { - out << t << ',' << id << ',' << x << ',' << y << ',' << e_x << ',' << e_y << ',' << c << '\n'; - } - } -#endif - - return 0; -} \ No newline at end of file diff --git a/variants/variant_float_inlining_SoA_NEON_polyexp_morevect_ilp_2_rsqrt_morersqrt.cpp b/variants/variant_float_inlining_SoA_NEON_polyexp_morevect_ilp_2_rsqrt_morersqrt.cpp index 14ce633..76eaed7 100644 --- a/variants/variant_float_inlining_SoA_NEON_polyexp_morevect_ilp_2_rsqrt_morersqrt.cpp +++ b/variants/variant_float_inlining_SoA_NEON_polyexp_morevect_ilp_2_rsqrt_morersqrt.cpp @@ -187,6 +187,7 @@ static inline float32x4_t _neon_exp_ps(float32x4_t x) { // # 4 MULTs, 4 ADDs float32x4_t r = vfmsq_f32(x, nf, ln2); // r = x - nf*ln2 + const float32x4_t c13 = vdupq_n_f32(c13_scalar); const float32x4_t c12 = vdupq_n_f32(c12_scalar); const float32x4_t c11 = vdupq_n_f32(c11_scalar); const float32x4_t c10 = vdupq_n_f32(c10_scalar); @@ -201,7 +202,7 @@ static inline float32x4_t _neon_exp_ps(float32x4_t x) { const float32x4_t c1 = vdupq_n_f32(c1_scalar); const float32x4_t c0 = vdupq_n_f32(c0_scalar); - float32x4_t p = c12; + float32x4_t p = c13; // # 13 * (4 MULTs, 4 ADDs) // # = 52 MULTs, 52 ADDs diff --git a/variants/variant_float_inlining_SoA_flopcount.cpp b/variants/variant_float_inlining_SoA_flopcount.cpp deleted file mode 100644 index 1842150..0000000 --- a/variants/variant_float_inlining_SoA_flopcount.cpp +++ /dev/null @@ -1,380 +0,0 @@ -// # Changes: SoA layout for pedestrians (file-scope arrays instead of struct fields) -// # Predecessor: variant_float_inlining.cpp - -#include "../headers/board_reduced_float.hpp" -#include "../headers/constants_float.hpp" -#include "../headers/pedestrian_reduced_float.hpp" - -#include -#include -#include -#include -#include -#include -#include -#include - -#define SEED 42 -#define BOARD_WIDTH 10.0f -#define BOARD_LENGTH 50.0f -#define WALL_EPS 1e-3f -#define BOARD_WIDTH_MINUS_WALL_EPS 9.999f - -#ifndef RANDOM_SEED -#define RANDOM_SEED 0 -#endif - -#ifndef BENCHMARK -#define BENCHMARK 0 -#endif - -using namespace std; - -// we wanna use arrays so we can implement SoA, one array per pedestrian field -vector pos_x, pos_y; -vector sb_sq; -vector sb_e_x, sb_e_y; -vector e_x, e_y; -vector curr_vel_x, curr_vel_y; -vector w_x, w_y; -vector desired_v0; -vector v_max; -vector destination_x, destination_y; -vector color; - -Pedestrian *init_pedestrians(int n_p) { -#if RANDOM_SEED - random_device rd; - mt19937 gen(rd()); -#else - mt19937 gen(SEED); -#endif - - normal_distribution velocity_gen(v0_mean, v0_std_dev); - uniform_real_distribution y_gen(0.3, BOARD_WIDTH - 0.3); - uniform_real_distribution x_gen_left(2.0, BOARD_LENGTH / 2.0); - uniform_real_distribution x_gen_right(BOARD_LENGTH / 2.0, BOARD_LENGTH - 2.0); - - // resize all arrays for this run - pos_x.resize(n_p); - pos_y.resize(n_p); - sb_sq.resize(n_p); - sb_e_x.resize(n_p); - sb_e_y.resize(n_p); - e_x.resize(n_p); - e_y.resize(n_p); - curr_vel_x.resize(n_p); - curr_vel_y.resize(n_p); - w_x.resize(n_p); - w_y.resize(n_p); - desired_v0.resize(n_p); - v_max.resize(n_p); - destination_x.resize(n_p); - destination_y.resize(n_p); - color.resize(n_p); - - int n_p_half = n_p / 2; - - for (int i = 0; i < n_p_half; i++) { - float x = x_gen_left(gen); - float y = y_gen(gen); - - float dest_x = BOARD_LENGTH + 30.0; // renamed to avoid shadowing destination_x[] - float dest_y = y; - - float init_e_x = dest_x - x; - float init_e_y = dest_y - y; - float init_e_norm = sqrt(init_e_x * init_e_x + init_e_y * init_e_y); - float init_e_unit_x = init_e_x / init_e_norm; - float init_e_unit_y = init_e_y / init_e_norm; - - float d_v0; - do { - d_v0 = velocity_gen(gen); - } while (d_v0 <= 0.1); - - float v_m = 1.3 * d_v0; - - color[i] = 0; - pos_x[i] = x; - pos_y[i] = y; - e_x[i] = init_e_unit_x; - e_y[i] = init_e_unit_y; - v_max[i] = v_m; - desired_v0[i] = d_v0; - destination_x[i] = dest_x; - destination_y[i] = dest_y; - - // Fields the Pedestrian constructor would have zero-initialized (c_v_x=c_v_y=0): - w_x[i] = 0.0f; - w_y[i] = 0.0f; - curr_vel_x[i] = 0.0f; - curr_vel_y[i] = 0.0f; - sb_sq[i] = 0.0f; - sb_e_x[i] = 0.0f; - sb_e_y[i] = 0.0f; - } - - for (int i = n_p_half; i < n_p; i++) { - float x = x_gen_right(gen); - float y = y_gen(gen); - - float dest_x = -30; - float dest_y = y; - - float init_e_x = dest_x - x; - float init_e_y = dest_y - y; - float init_e_norm = sqrt(init_e_x * init_e_x + init_e_y * init_e_y); - float init_e_unit_x = init_e_x / init_e_norm; - float init_e_unit_y = init_e_y / init_e_norm; - - float d_v0; - do { - d_v0 = velocity_gen(gen); - } while (d_v0 <= 0.1); - - float v_m = 1.3 * d_v0; - - color[i] = 1; - pos_x[i] = x; - pos_y[i] = y; - e_x[i] = init_e_unit_x; - e_y[i] = init_e_unit_y; - v_max[i] = v_m; - desired_v0[i] = d_v0; - destination_x[i] = dest_x; - destination_y[i] = dest_y; - - w_x[i] = 0.0f; - w_y[i] = 0.0f; - curr_vel_x[i] = 0.0f; - curr_vel_y[i] = 0.0f; - sb_sq[i] = 0.0f; - sb_e_x[i] = 0.0f; - sb_e_y[i] = 0.0f; - } - - return nullptr; -} - -Board init_board_without_obstacles() { return Board(BOARD_LENGTH, BOARD_WIDTH, 0, nullptr); } - -int run_simulation(Pedestrian *pedestrians, Board &board, int n_p, int n_iter, const char *dump_path) { - (void)pedestrians; -#if BENCHMARK - (void)dump_path; -#endif - -#if !BENCHMARK - const int dump_interval = 1; - - vector> samples; - - if (dump_path) { - samples.reserve(static_cast(n_iter / dump_interval + 1) * n_p); - - for (int i = 0; i < n_p; i++) { - samples.emplace_back(0, i, pos_x[i], pos_y[i], e_x[i], e_y[i], color[i]); - } - } -#endif - - // # ADDs: n_iter * ( - // # n_p * 4 + - // # n_p * (n_p - 1) * 14 + - // # n_p * board.n_obstacles * 10 + - // # n_p * 7 - // # ) - - // # MULTs: n_iter * ( - // # n_p * 4 - // # n_p * (n_p - 1) * 24 + - // # n_p * board.n_obstacles * 13 + - // # n_p * 17 - // # ) - - // # DIVs: n_iter * ( - // # n_p * (n_p - 1) + - // # n_p * board.n_obstacles + - // # n_p * 2 - // # ) - - // # SQRTs: n_iter * ( - // # n_p * (n_p - 1) * 3 + - // # n_p * board.n_obstacles + - // # n_p * 3 - // # ) - - // # EXPs: n_iter * ( - // # n_p * (n_p - 1) + - // # n_p * board.n_obstacles - // # ) - - for (int t = 0; t < n_iter; t++) { - for (int i = 0; i < n_p; i++) { - - float desired_velocity_x = e_x[i] * desired_v0[i]; - float desired_velocity_y = e_y[i] * desired_v0[i]; - float f_alpha_dt_x = (desired_velocity_x - curr_vel_x[i]) * dt_times_tau_alpha_inv; - float f_alpha_dt_y = (desired_velocity_y - curr_vel_y[i]) * dt_times_tau_alpha_inv; - - for (int j = 0; j < n_p; j++) { - if (i == j) - continue; - - float r_ab_x = pos_x[i] - pos_x[j]; - float r_ab_y = pos_y[i] - pos_y[j]; - - float r_ab_norm = sqrt(r_ab_x * r_ab_x + r_ab_y * r_ab_y); - r_ab_norm = r_ab_norm > EPS ? r_ab_norm : EPS; - - float r_ab2_x = r_ab_x - sb_e_x[j]; - float r_ab2_y = r_ab_y - sb_e_y[j]; - - float r_ab2_norm = sqrt(r_ab2_x * r_ab2_x + r_ab2_y * r_ab2_y); - r_ab2_norm = r_ab2_norm > EPS ? r_ab2_norm : EPS; - - float r_ab_norm_sum = r_ab_norm + r_ab2_norm; - - float two_b_sq = r_ab_norm_sum * r_ab_norm_sum - sb_sq[j]; - two_b_sq = two_b_sq > EPS_SQ ? two_b_sq : EPS_SQ; - float two_b = sqrt(two_b_sq); - float scalar_mult_common = (r_ab_norm_sum / (two_b * r_ab_norm * r_ab2_norm)) * (NEG_V0_AB_div_TWO_SIGMA * exp(two_b * NEG_TWO_SIGMA_INV)); - float scalar_mult_r_ab = scalar_mult_common * r_ab2_norm; - float scalar_mult_r_ab2 = scalar_mult_common * r_ab_norm; - - float f_alpha_beta_neg_x = r_ab_x * scalar_mult_r_ab + r_ab2_x * scalar_mult_r_ab2; - float f_alpha_beta_neg_y = r_ab_y * scalar_mult_r_ab + r_ab2_y * scalar_mult_r_ab2; - - float w_val; - - float f_norm_sq = f_alpha_beta_neg_x * f_alpha_beta_neg_x + f_alpha_beta_neg_y * f_alpha_beta_neg_y; - float rhs = f_norm_sq * cos_phi_sq; - float e_dot_f = e_x[i] * f_alpha_beta_neg_x + e_y[i] * f_alpha_beta_neg_y; - float e_dot_f_sq = e_dot_f * e_dot_f; - if (e_dot_f >= 0.0 || e_dot_f_sq <= rhs) - w_val = dt_neg; - else - w_val = c_neg_times_dt; - - f_alpha_dt_x = f_alpha_dt_x + f_alpha_beta_neg_x * w_val; - f_alpha_dt_y = f_alpha_dt_y + f_alpha_beta_neg_y * w_val; - } - - for (int k = 0; k < board.n_obstacles; k++) { - - float v1_x = pos_x[i] - board.obstacles[k].start_x; - float v1_y = pos_y[i] - board.obstacles[k].start_y; - float v2_x = board.obstacles[k].obstacle_v_x; - float v2_y = board.obstacles[k].obstacle_v_y; - - float projection = (v1_x * v2_x + v1_y * v2_y) * board.obstacles[k].length_sq_inv; - - if (projection < 0.0) - projection = 0.0; - else if (projection > 1.0) - projection = 1.0; - - float r_aB_x = pos_x[i] - (board.obstacles[k].start_x + board.obstacles[k].obstacle_v_x * projection); - float r_aB_y = pos_y[i] - (board.obstacles[k].start_y + board.obstacles[k].obstacle_v_y * projection); - - float r_aB_norm = sqrt(r_aB_x * r_aB_x + r_aB_y * r_aB_y); - r_aB_norm = r_aB_norm > EPS ? r_aB_norm : EPS; - - float scalar_mult = (NEG_U0_AB_div_R * exp(r_aB_norm * NEG_R_INV)) / r_aB_norm; - float f_alpha_B_neg_x = r_aB_x * scalar_mult; - float f_alpha_B_neg_y = r_aB_y * scalar_mult; - - float w_val; - - float e_dot_f = (f_alpha_B_neg_y == 0) ? 0 : (f_alpha_B_neg_y < 0 ? -e_y[i] : e_y[i]); - - if (e_dot_f >= cos_phi) - w_val = dt_neg; - else - w_val = c_neg_times_dt; - - float f_alpha_B_dt_x = f_alpha_B_neg_x * w_val; - float f_alpha_B_dt_y = f_alpha_B_neg_y * w_val; - - f_alpha_dt_x += f_alpha_B_dt_x; - f_alpha_dt_y += f_alpha_B_dt_y; - } - - w_x[i] += f_alpha_dt_x; - w_y[i] += f_alpha_dt_y; - } - - for (int i = 0; i < n_p; i++) { - - float w_norm = sqrt(w_x[i] * w_x[i] + w_y[i] * w_y[i]); - - float w_norm2 = w_norm > EPS ? w_norm : EPS; - - float correction_factor = min(1.0f, v_max[i] / w_norm2); - - curr_vel_x[i] = w_x[i] * correction_factor; - curr_vel_y[i] = w_y[i] * correction_factor; - - pos_x[i] = pos_x[i] + curr_vel_x[i] * dt; - pos_y[i] = pos_y[i] + curr_vel_y[i] * dt; - - if (color[i] == 0 && pos_x[i] > (BOARD_LENGTH - 1.5)) { - pos_x[i] = 1.0; - destination_y[i] = pos_y[i]; - w_x[i] = desired_v0[i]; - w_y[i] = 0.0; - curr_vel_x[i] = w_x[i]; - curr_vel_y[i] = w_y[i]; - } - if (color[i] == 1 && pos_x[i] < 1.5) { - pos_x[i] = BOARD_LENGTH - 1.0; - destination_y[i] = pos_y[i]; - w_x[i] = -desired_v0[i]; - w_y[i] = 0.0; - curr_vel_x[i] = w_x[i]; - curr_vel_y[i] = w_y[i]; - } - - pos_y[i] = clamp(pos_y[i], WALL_EPS, BOARD_WIDTH_MINUS_WALL_EPS); - - float ex = destination_x[i] - pos_x[i]; - float ey = destination_y[i] - pos_y[i]; - - float e_norm = sqrt(ex * ex + ey * ey); - float e_norm2 = e_norm > EPS ? e_norm : EPS; - float e_norm_inv = 1.0 / e_norm2; - - e_x[i] = ex * e_norm_inv; - e_y[i] = ey * e_norm_inv; - - float velocity_norm_sq = curr_vel_x[i] * curr_vel_x[i] + curr_vel_y[i] * curr_vel_y[i]; - - float velocity_norm = sqrt(velocity_norm_sq); - - sb_sq[i] = velocity_norm_sq * delta_t_sq; - - sb_e_x[i] = e_x[i] * (velocity_norm * delta_t); - sb_e_y[i] = e_y[i] * (velocity_norm * delta_t); - -#if !BENCHMARK - if (dump_path && (t % dump_interval == 0)) { - samples.emplace_back(t + 1, i, pos_x[i], pos_y[i], e_x[i], e_y[i], color[i]); - } -#endif - } - } - -#if !BENCHMARK - if (dump_path) { - ofstream out(dump_path); - out << "time,id,x,y,p.e.x,p.e.y,color\n"; - - for (auto &[t, id, x, y, e_x, e_y, c] : samples) { - out << t << ',' << id << ',' << x << ',' << y << ',' << e_x << ',' << e_y << ',' << c << '\n'; - } - } -#endif - - return 0; -} diff --git a/variants/variant_float_inlining_flopcount.cpp b/variants/variant_float_inlining_flopcount.cpp deleted file mode 100644 index 482b915..0000000 --- a/variants/variant_float_inlining_flopcount.cpp +++ /dev/null @@ -1,366 +0,0 @@ -// # Changes: Replaced «double» with «float» in variable declarations -// # Predecessor: variant_float_baseline.cpp - -#include "../headers/board_reduced_float.hpp" -#include "../headers/constants_float.hpp" -#include "../headers/pedestrian_reduced_float.hpp" - -#include -#include -#include -#include -#include -#include -#include -#include - -#define SEED 42 -#define BOARD_WIDTH 10.0f -#define BOARD_LENGTH 50.0f -#define WALL_EPS 1e-3f -#define BOARD_WIDTH_MINUS_WALL_EPS 9.999f - -#ifndef RANDOM_SEED -#define RANDOM_SEED 0 -#endif - -#ifndef BENCHMARK -#define BENCHMARK 0 -#endif - -using namespace std; - -Pedestrian *init_pedestrians(int n_p) { -#if RANDOM_SEED - random_device rd; - mt19937 gen(rd()); -#else - mt19937 gen(SEED); -#endif - - normal_distribution velocity_gen(v0_mean, v0_std_dev); - uniform_real_distribution y_gen(0.3, BOARD_WIDTH - 0.3); - uniform_real_distribution x_gen_left(2.0, BOARD_LENGTH / 2.0); - uniform_real_distribution x_gen_right(BOARD_LENGTH / 2.0, BOARD_LENGTH - 2.0); - - Pedestrian *P = new Pedestrian[n_p]; - - int n_p_half = n_p / 2; - - for (int i = 0; i < n_p_half; i++) { - float x = x_gen_left(gen); - float y = y_gen(gen); - - float start_pos_x = x; - float start_pos_y = y; - - float destination_x = BOARD_LENGTH + 30.0; - float destination_y = y; - - float init_e_x = destination_x - start_pos_x; - float init_e_y = destination_y - start_pos_y; - float init_e_norm = sqrt(init_e_x * init_e_x + init_e_y * init_e_y); - float init_e_unit_x = init_e_x / init_e_norm; - float init_e_unit_y = init_e_y / init_e_norm; - - float desired_v0; - do { - desired_v0 = velocity_gen(gen); - } while (desired_v0 <= 0.1); - - float v_max = 1.3 * desired_v0; - - P[i] = Pedestrian(0, 0, 0, start_pos_x, start_pos_y, init_e_unit_x, init_e_unit_y, v_max, desired_v0, destination_x, destination_y); - } - - for (int i = n_p_half; i < n_p; i++) { - float x = x_gen_right(gen); - float y = y_gen(gen); - - float start_pos_x = x; - float start_pos_y = y; - float destination_x = -30; - float destination_y = y; - - float init_e_x = destination_x - start_pos_x; - float init_e_y = destination_y - start_pos_y; - float init_e_norm = sqrt(init_e_x * init_e_x + init_e_y * init_e_y); - float init_e_unit_x = init_e_x / init_e_norm; - float init_e_unit_y = init_e_y / init_e_norm; - - float desired_v0; - do { - desired_v0 = velocity_gen(gen); - } while (desired_v0 <= 0.1); - - float v_max = 1.3 * desired_v0; - - P[i] = Pedestrian(1, 0, 0, start_pos_x, start_pos_y, init_e_unit_x, init_e_unit_y, v_max, desired_v0, destination_x, destination_y); - } - - return P; -} - -Board init_board_without_obstacles() { return Board(BOARD_LENGTH, BOARD_WIDTH, 0, nullptr); } - -int run_simulation(Pedestrian *pedestrians, Board &board, int n_p, int n_iter, const char *dump_path) { - -#if BENCHMARK - (void)dump_path; -#endif - -#if !BENCHMARK - const int dump_interval = 1; - - vector> samples; - - if (dump_path) { - samples.reserve(static_cast(n_iter / dump_interval + 1) * n_p); - - for (int i = 0; i < n_p; i++) { - Pedestrian &p = pedestrians[i]; - samples.emplace_back(0, i, p.curr_pos_x, p.curr_pos_y, p.e_x, p.e_y, p.color); - } - } -#endif - - // # ADDs: n_iter * ( - // # n_p * 4 + - // # n_p * (n_p - 1) * 14 + - // # n_p * board.n_obstacles * 10 + - // # n_p * 7 - // # ) - - // # MULTs: n_iter * ( - // # n_p * 4 - // # n_p * (n_p - 1) * 24 + - // # n_p * board.n_obstacles * 13 + - // # n_p * 17 - // # ) - - // # DIVs: n_iter * ( - // # n_p * (n_p - 1) + - // # n_p * board.n_obstacles + - // # n_p * 2 - // # ) - - // # SQRTs: n_iter * ( - // # n_p * (n_p - 1) * 3 + - // # n_p * board.n_obstacles + - // # n_p * 3 - // # ) - - // # EXPs: n_iter * ( - // # n_p * (n_p - 1) + - // # n_p * board.n_obstacles - // # ) - - for (int t = 0; t < n_iter; t++) { - for (int i = 0; i < n_p; i++) { - Pedestrian &p = pedestrians[i]; - - // # 2 MULTs - float desired_velocity_x = p.e_x * p.desired_v0; - float desired_velocity_y = p.e_y * p.desired_v0; - - // # 2 MULTs, 2 ADDs - float f_alpha_dt_x = (desired_velocity_x - p.curr_vel_x) * dt_times_tau_alpha_inv; - float f_alpha_dt_y = (desired_velocity_y - p.curr_vel_y) * dt_times_tau_alpha_inv; - - for (int j = 0; j < n_p; j++) { - if (i == j) - continue; - - // # 2 ADDs - float r_ab_x = p.curr_pos_x - pedestrians[j].curr_pos_x; - float r_ab_y = p.curr_pos_y - pedestrians[j].curr_pos_y; - - // # 1 SQRT, 2 MULTs, 1 ADD - float r_ab_norm = sqrt(r_ab_x * r_ab_x + r_ab_y * r_ab_y); - r_ab_norm = r_ab_norm > EPS ? r_ab_norm : EPS; - - // # 2 ADDs - float r_ab2_x = r_ab_x - pedestrians[j].s_b_e_x; - float r_ab2_y = r_ab_y - pedestrians[j].s_b_e_y; - - // # 1 SQRT, 2 MULTs, 1 ADD - float r_ab2_norm = sqrt(r_ab2_x * r_ab2_x + r_ab2_y * r_ab2_y); - r_ab2_norm = r_ab2_norm > EPS ? r_ab2_norm : EPS; - - // # 1 ADD - float r_ab_norm_sum = r_ab_norm + r_ab2_norm; - - // # 1 MULT 1 ADD - float two_b_sq = r_ab_norm_sum * r_ab_norm_sum - pedestrians[j].s_b_sq; - two_b_sq = two_b_sq > EPS_SQ ? two_b_sq : EPS_SQ; - - // # 1 SQRT - float two_b = sqrt(two_b_sq); - - // # 5 MULTs, 1 DIV, 1 EXP - float scalar_mult_common = (r_ab_norm_sum / (two_b * r_ab_norm * r_ab2_norm)) * (NEG_V0_AB_div_TWO_SIGMA * exp(two_b * NEG_TWO_SIGMA_INV)); - - // # 2 MULTs - float scalar_mult_r_ab = scalar_mult_common * r_ab2_norm; - float scalar_mult_r_ab2 = scalar_mult_common * r_ab_norm; - - // # 4 MULTs, 2 ADDs - float f_alpha_beta_neg_x = r_ab_x * scalar_mult_r_ab + r_ab2_x * scalar_mult_r_ab2; - float f_alpha_beta_neg_y = r_ab_y * scalar_mult_r_ab + r_ab2_y * scalar_mult_r_ab2; - - float w_val; - - // # 6 MULTs, 2 ADDs - float f_norm_sq = f_alpha_beta_neg_x * f_alpha_beta_neg_x + f_alpha_beta_neg_y * f_alpha_beta_neg_y; - float rhs = f_norm_sq * cos_phi_sq; - float e_dot_f = p.e_x * f_alpha_beta_neg_x + p.e_y * f_alpha_beta_neg_y; - float e_dot_f_sq = e_dot_f * e_dot_f; - - if (e_dot_f >= 0.0 || e_dot_f_sq <= rhs) - w_val = dt_neg; - else - w_val = c_neg_times_dt; - - // # 2 MULTs, 2 ADDs - f_alpha_dt_x = f_alpha_dt_x + f_alpha_beta_neg_x * w_val; - f_alpha_dt_y = f_alpha_dt_y + f_alpha_beta_neg_y * w_val; - } - - for (int k = 0; k < board.n_obstacles; k++) { - - // # 2 ADDs - float v1_x = p.curr_pos_x - board.obstacles[k].start_x; - float v1_y = p.curr_pos_y - board.obstacles[k].start_y; - float v2_x = board.obstacles[k].obstacle_v_x; - float v2_y = board.obstacles[k].obstacle_v_y; - - // # 3 MULTs, 1 ADD - float projection = (v1_x * v2_x + v1_y * v2_y) * board.obstacles[k].length_sq_inv; - - if (projection < 0.0) - projection = 0.0; - else if (projection > 1.0) - projection = 1.0; - - // # 4 ADDs, 2 MULTs - float r_aB_x = p.curr_pos_x - (board.obstacles[k].start_x + board.obstacles[k].obstacle_v_x * projection); - float r_aB_y = p.curr_pos_y - (board.obstacles[k].start_y + board.obstacles[k].obstacle_v_y * projection); - - // # 1 SQRT, 2 MULTs, 1 ADD - float r_aB_norm = sqrt(r_aB_x * r_aB_x + r_aB_y * r_aB_y); - r_aB_norm = r_aB_norm > EPS ? r_aB_norm : EPS; - - // # 4 MULTs, 1 EXP, 1 DIV - float scalar_mult = (NEG_U0_AB_div_R * exp(r_aB_norm * NEG_R_INV)) / r_aB_norm; - float f_alpha_B_neg_x = r_aB_x * scalar_mult; - float f_alpha_B_neg_y = r_aB_y * scalar_mult; - - float w_val; - - float e_dot_f = (f_alpha_B_neg_y == 0) ? 0 : (f_alpha_B_neg_y < 0 ? -p.e_y : p.e_y); - - if (e_dot_f >= cos_phi) - w_val = dt_neg; - else - w_val = c_neg_times_dt; - - // # 2 MULTS - float f_alpha_B_dt_x = f_alpha_B_neg_x * w_val; - float f_alpha_B_dt_y = f_alpha_B_neg_y * w_val; - - // # 2 ADDs - f_alpha_dt_x += f_alpha_B_dt_x; - f_alpha_dt_y += f_alpha_B_dt_y; - } - - // # 2 ADDs - p.w_x += f_alpha_dt_x; - p.w_y += f_alpha_dt_y; - } - - for (int i = 0; i < n_p; i++) { - Pedestrian &p = pedestrians[i]; - - // # 1 SQRT, 2 MULTs, 1 ADD - float w_norm = sqrt(p.w_x * p.w_x + p.w_y * p.w_y); - - float w_norm2 = w_norm > EPS ? w_norm : EPS; - - // # 1 DIV - float correction_factor = min(1.0f, p.v_max / w_norm2); - - // # 2 MULTs - p.curr_vel_x = p.w_x * correction_factor; - p.curr_vel_y = p.w_y * correction_factor; - - // # 2 MULTs, 2 ADDs - p.curr_pos_x = p.curr_pos_x + p.curr_vel_x * dt; - p.curr_pos_y = p.curr_pos_y + p.curr_vel_y * dt; - - if (p.color == 0 && p.curr_pos_x > (BOARD_LENGTH - 1.5)) { - p.curr_pos_x = 1.0; - p.destination_y = p.curr_pos_y; - p.w_x = p.desired_v0; - p.w_y = 0.0; - p.curr_vel_x = p.w_x; - p.curr_vel_y = p.w_y; - } - if (p.color == 1 && p.curr_pos_x < 1.5) { - p.curr_pos_x = BOARD_LENGTH - 1.0; - p.destination_y = p.curr_pos_y; - p.w_x = -p.desired_v0; - p.w_y = 0.0; - p.curr_vel_x = p.w_x; - p.curr_vel_y = p.w_y; - } - - p.curr_pos_y = clamp(p.curr_pos_y, WALL_EPS, BOARD_WIDTH_MINUS_WALL_EPS); - - // # 2 ADDs - float e_x = p.destination_x - p.curr_pos_x; - float e_y = p.destination_y - p.curr_pos_y; - - // # 1 SQRT, 2 MULTs, 1 ADD - float e_norm = sqrt(e_x * e_x + e_y * e_y); - - float e_norm2 = e_norm > EPS ? e_norm : EPS; - - // # 1 DIV - float e_norm_inv = 1.0 / e_norm2; - - // # 2 MULTs - p.e_x = e_x * e_norm_inv; - p.e_y = e_y * e_norm_inv; - - // # 2 MULTs, 1 ADD - float velocity_norm_sq = p.curr_vel_x * p.curr_vel_x + p.curr_vel_y * p.curr_vel_y; - - // # 1 SQRT - float velocity_norm = sqrt(velocity_norm_sq); - - // # 5 MULTs - p.s_b_sq = velocity_norm_sq * delta_t_sq; - p.s_b_e_x = p.e_x * (velocity_norm * delta_t); - p.s_b_e_y = p.e_y * (velocity_norm * delta_t); - -#if !BENCHMARK - if (dump_path && (t % dump_interval == 0)) { - samples.emplace_back(t + 1, i, p.curr_pos_x, p.curr_pos_y, p.e_x, p.e_y, p.color); - } -#endif - } - } - -#if !BENCHMARK - if (dump_path) { - ofstream out(dump_path); - out << "time,id,x,y,p.e.x,p.e.y,color\n"; - - for (auto &[t, id, x, y, e_x, e_y, c] : samples) { - out << t << ',' << id << ',' << x << ',' << y << ',' << e_x << ',' << e_y << ',' << c << '\n'; - } - } -#endif - - return 0; -} diff --git a/variants/variant_inlining_7_SoA_flopcount.cpp b/variants/variant_inlining_7_SoA_flopcount.cpp deleted file mode 100644 index 8c10a6a..0000000 --- a/variants/variant_inlining_7_SoA_flopcount.cpp +++ /dev/null @@ -1,362 +0,0 @@ -// # This file uses avx2 -// # This is the first variant and should serve as a base for further avx2 optimizations - -// Okay our first step is to use a SoA instead of an AoS -// We only do SoA for the pedestrian struct, there we wanna optimize such that we access for the ped-ped forces adjacent pedestrians (in the arrays) - -#include "../headers/board.hpp" -#include "../headers/constants.hpp" -#include "../headers/pedestrian_reduced.hpp" - -#include -#include -#include -#include -#include -#include -#include -#include - -#define SEED 42 -#define BOARD_WIDTH 10.0 -#define BOARD_LENGTH 50.0 -#define WALL_EPS 1e-3 -#define BOARD_WIDTH_MINUS_WALL_EPS 9.999 - -#ifndef RANDOM_SEED -#define RANDOM_SEED 0 -#endif - -#ifndef BENCHMARK -#define BENCHMARK 0 -#endif - -using namespace std; - -// we wanna use arrays so we can implement SOA, we create arrays for all the data we need and access -vector pos_x, pos_y; -vector sb_sq; -vector sb_e_x, sb_e_y; -vector e_x, e_y; -vector curr_vel_x, curr_vel_y; -vector w_x, w_y; -vector desired_v0; -vector v_max; -vector destination_x, destination_y; -vector color; - -Pedestrian *init_pedestrians(int n_p) { - -#if RANDOM_SEED - random_device rd; - mt19937 gen(rd()); -#else - mt19937 gen(SEED); -#endif - - normal_distribution velocity_gen(v0_mean, v0_std_dev); - uniform_real_distribution y_gen(0.3, BOARD_WIDTH - 0.3); - uniform_real_distribution x_gen_left(2.0, BOARD_LENGTH / 2.0); - uniform_real_distribution x_gen_right(BOARD_LENGTH / 2.0, BOARD_LENGTH - 2.0); - - // resize all our arrays for this run - pos_x.resize(n_p); - pos_y.resize(n_p); - sb_sq.resize(n_p); - sb_e_x.resize(n_p); - sb_e_y.resize(n_p); - e_x.resize(n_p); - e_y.resize(n_p); - curr_vel_x.resize(n_p); - curr_vel_y.resize(n_p); - w_x.resize(n_p); - w_y.resize(n_p); - desired_v0.resize(n_p); - v_max.resize(n_p); - destination_x.resize(n_p); - destination_y.resize(n_p); - color.resize(n_p); - - int n_p_half = n_p / 2; - - for (int i = 0; i < n_p_half; i++) { - double x = x_gen_left(gen); - double y = y_gen(gen); - - double dest_x = BOARD_LENGTH + 30.0; // renamed to avoid shadowing destination_x[] - double dest_y = y; // same idea for destination_y[] - - double init_e_x = dest_x - x; - double init_e_y = dest_y - y; - double init_e_norm = sqrt(init_e_x * init_e_x + init_e_y * init_e_y); - double init_e_unit_x = init_e_x / init_e_norm; - double init_e_unit_y = init_e_y / init_e_norm; - - double d_v0; - do { - d_v0 = velocity_gen(gen); - } while (d_v0 <= 0.1); - - double v_m = 1.3 * d_v0; - - color[i] = 0; - pos_x[i] = x; - pos_y[i] = y; - e_x[i] = init_e_unit_x; - e_y[i] = init_e_unit_y; - v_max[i] = v_m; - desired_v0[i] = d_v0; - destination_x[i] = dest_x; - destination_y[i] = dest_y; - - // Fields the Pedestrian constructor would have zero-initialized: - w_x[i] = 0.0; - w_y[i] = 0.0; - curr_vel_x[i] = 0.0; - curr_vel_y[i] = 0.0; - sb_sq[i] = 0.0; - sb_e_x[i] = 0.0; - sb_e_y[i] = 0.0; - } - - for (int i = n_p_half; i < n_p; i++) { - double x = x_gen_right(gen); - double y = y_gen(gen); - - double dest_x = -30; - double dest_y = y; - - double init_e_x = dest_x - x; - double init_e_y = dest_y - y; - double init_e_norm = sqrt(init_e_x * init_e_x + init_e_y * init_e_y); - double init_e_unit_x = init_e_x / init_e_norm; - double init_e_unit_y = init_e_y / init_e_norm; - - double d_v0; - do { - d_v0 = velocity_gen(gen); - } while (d_v0 <= 0.1); - - double v_m = 1.3 * d_v0; - - color[i] = 1; - pos_x[i] = x; - pos_y[i] = y; - e_x[i] = init_e_unit_x; - e_y[i] = init_e_unit_y; - v_max[i] = v_m; - desired_v0[i] = d_v0; - destination_x[i] = dest_x; - destination_y[i] = dest_y; - - // Fields the Pedestrian constructor would have zero-initialized: - w_x[i] = 0.0; - w_y[i] = 0.0; - curr_vel_x[i] = 0.0; - curr_vel_y[i] = 0.0; - sb_sq[i] = 0.0; - sb_e_x[i] = 0.0; - sb_e_y[i] = 0.0; - } - - return nullptr; -} - -Board init_board_without_obstacles() { return Board(BOARD_LENGTH, BOARD_WIDTH, 0, nullptr); } - -int run_simulation(Pedestrian *pedestrians, Board &board, int n_p, int n_iter, const char *dump_path) { - -#if BENCHMARK - (void)dump_path; -#endif - -#if !BENCHMARK - const int dump_interval = 1; - - vector> samples; - - if (dump_path) { - samples.reserve(static_cast(n_iter / dump_interval + 1) * n_p); - - for (int i = 0; i < n_p; i++) { - // Pedestrian &p = pedestrians[i]; - samples.emplace_back(0, i, pos_x[i], pos_y[i], e_x[i], e_y[i], color[i]); - } - } -#endif - - // # ADDs: n_iter * n_p * 11 + - // # n_iter * n_p * (n_p - 1) * 14 + - // # n_iter * n_p * board.n_obstacles * 12 - - // # MULTs: n_iter * n_p * 21 + - // # n_iter * n_p * (n_p - 1) * 24 + - // # n_iter * n_p * board.n_obstacles * 19 - - // # DIVs: n_iter * n_p * 2 + - // # n_iter * n_p * (n_p - 1) + - // # n_iter * n_p * board.n_obstacles - - // # SQRTs: n_iter * n_p * 3 + - // # n_iter * n_p * (n_p - 1) * 3 + - // # n_iter * n_p * board.n_obstacles - - // # EXPs: n_iter * n_p * (n_p - 1) + - // # n_iter * n_p * board.n_obstacles - - for (int t = 0; t < n_iter; t++) { - - for (int i = 0; i < n_p; i++) { - - double desired_velocity_x = e_x[i] * desired_v0[i]; - double desired_velocity_y = e_y[i] * desired_v0[i]; - - double f_alpha_dt_x = (desired_velocity_x - curr_vel_x[i]) * dt_times_tau_alpha_inv; - double f_alpha_dt_y = (desired_velocity_y - curr_vel_y[i]) * dt_times_tau_alpha_inv; - - for (int j = 0; j < n_p; j++) { - - if (i == j) - continue; - - double r_ab_x = pos_x[i] - pos_x[j]; - double r_ab_y = pos_y[i] - pos_y[j]; - - double r_ab_norm = sqrt(r_ab_x * r_ab_x + r_ab_y * r_ab_y); - - double r_ab2_x = r_ab_x - sb_e_x[j]; - double r_ab2_y = r_ab_y - sb_e_y[j]; - - double r_ab2_norm = sqrt(r_ab2_x * r_ab2_x + r_ab2_y * r_ab2_y); - double r_ab_norm_sum = r_ab_norm + r_ab2_norm; - double two_b_sq = r_ab_norm_sum * r_ab_norm_sum - sb_sq[j]; - double two_b = sqrt(two_b_sq); - double scalar_mult_common = (r_ab_norm_sum / (two_b * r_ab_norm * r_ab2_norm)) * (NEG_V0_AB_div_TWO_SIGMA * exp(two_b * NEG_TWO_SIGMA_INV)); - double scalar_mult_r_ab = scalar_mult_common * r_ab2_norm; - double scalar_mult_r_ab2 = scalar_mult_common * r_ab_norm; - - double f_alpha_beta_neg_x = r_ab_x * scalar_mult_r_ab + r_ab2_x * scalar_mult_r_ab2; - double f_alpha_beta_neg_y = r_ab_y * scalar_mult_r_ab + r_ab2_y * scalar_mult_r_ab2; - - double w_val; - double a = (e_x[i] * f_alpha_beta_neg_x + e_y[i] * f_alpha_beta_neg_y) * cos_phi_inv; - double b_sq = f_alpha_beta_neg_x * f_alpha_beta_neg_x + f_alpha_beta_neg_y * f_alpha_beta_neg_y; - - w_val = dt_neg; - if (a >= 0 && (a * a > b_sq)) { - w_val = c_neg_times_dt; - } - - f_alpha_dt_x = f_alpha_dt_x + f_alpha_beta_neg_x * w_val; - f_alpha_dt_y = f_alpha_dt_y + f_alpha_beta_neg_y * w_val; - } - - for (int k = 0; k < board.n_obstacles; k++) { - - double v1_x = pos_x[i] - board.obstacles[k].start.x; - double v1_y = pos_y[i] - board.obstacles[k].start.y; - - double projection = (v1_x * board.obstacles[k].obstacle_v.x + v1_y * board.obstacles[k].obstacle_v.y) * board.obstacles[k].length_sq_inv; - - if (projection < 0.0) - projection = 0.0; - else if (projection > 1.0) - projection = 1.0; - - double r_aB_x = pos_x[i] - (board.obstacles[k].start.x + board.obstacles[k].obstacle_v.x * projection); - double r_aB_y = pos_y[i] - (board.obstacles[k].start.y + board.obstacles[k].obstacle_v.y * projection); - - double r_aB_norm = sqrt(r_aB_x * r_aB_x + r_aB_y * r_aB_y); - double scalar_mult = (NEG_U0_AB_div_R * exp(r_aB_norm * NEG_R_INV)) / r_aB_norm; - double f_alpha_B_neg_x = r_aB_x * scalar_mult; - double f_alpha_B_neg_y = r_aB_y * scalar_mult; - - double w_val; - double a = (e_x[i] * f_alpha_B_neg_x + e_y[i] * f_alpha_B_neg_y) * cos_phi_inv; - double b_sq = f_alpha_B_neg_x * f_alpha_B_neg_x + f_alpha_B_neg_y * f_alpha_B_neg_y; - - w_val = dt_neg; - if (a >= 0 && (a * a > b_sq)) { - w_val = c_neg_times_dt; - } - - double f_alpha_B_dt_x = f_alpha_B_neg_x * w_val; - double f_alpha_B_dt_y = f_alpha_B_neg_y * w_val; - - f_alpha_dt_x += f_alpha_B_dt_x; - f_alpha_dt_y += f_alpha_B_dt_y; - } - - w_x[i] += f_alpha_dt_x; - w_y[i] += f_alpha_dt_y; - } - - for (int i = 0; i < n_p; i++) { - - double w_norm = sqrt(w_x[i] * w_x[i] + w_y[i] * w_y[i]); - - double correction_factor = min(1.0, v_max[i] / w_norm); - - curr_vel_x[i] = w_x[i] * correction_factor; - curr_vel_y[i] = w_y[i] * correction_factor; - - pos_x[i] = pos_x[i] + curr_vel_x[i] * dt; - pos_y[i] = pos_y[i] + curr_vel_y[i] * dt; - - if (color[i] == 0 && pos_x[i] > (BOARD_LENGTH - 1.5)) { - pos_x[i] = 1.0; - destination_y[i] = pos_y[i]; - w_x[i] = desired_v0[i]; - w_y[i] = 0.0; - curr_vel_x[i] = w_x[i]; - curr_vel_y[i] = w_y[i]; - } - if (color[i] == 1 && pos_x[i] < 1.5) { - pos_x[i] = BOARD_LENGTH - 1.0; - destination_y[i] = pos_y[i]; - w_x[i] = -desired_v0[i]; - w_y[i] = 0.0; - curr_vel_x[i] = w_x[i]; - curr_vel_y[i] = w_y[i]; - } - - pos_y[i] = clamp(pos_y[i], WALL_EPS, BOARD_WIDTH_MINUS_WALL_EPS); - - double ex = destination_x[i] - pos_x[i]; - double ey = destination_y[i] - pos_y[i]; - - double e_norm_inv = 1.0 / sqrt(ex * ex + ey * ey); - - e_x[i] = ex * e_norm_inv; - e_y[i] = ey * e_norm_inv; - - double velocity_norm_sq = curr_vel_x[i] * curr_vel_x[i] + curr_vel_y[i] * curr_vel_y[i]; - - double velocity_norm = sqrt(velocity_norm_sq); - - sb_sq[i] = velocity_norm_sq * delta_t_sq; - - sb_e_x[i] = e_x[i] * (velocity_norm * delta_t); - sb_e_y[i] = e_y[i] * (velocity_norm * delta_t); - -#if !BENCHMARK - if (dump_path && (t % dump_interval == 0)) { - samples.emplace_back(t + 1, i, pos_x[i], pos_y[i], e_x[i], e_y[i], color[i]); - } -#endif - } - } - -#if !BENCHMARK - if (dump_path) { - ofstream out(dump_path); - out << "time,id,x,y,p.e.x,p.e.y,color\n"; - - for (auto &[t, id, x, y, e_x, e_y, c] : samples) { - out << t << ',' << id << ',' << x << ',' << y << ',' << e_x << ',' << e_y << ',' << c << '\n'; - } - } -#endif - - return 0; -} diff --git a/variants/variant_inlining_7_flopcount.cpp b/variants/variant_inlining_7_flopcount.cpp deleted file mode 100644 index 9e5b28a..0000000 --- a/variants/variant_inlining_7_flopcount.cpp +++ /dev/null @@ -1,352 +0,0 @@ -// # This file performs in-lining for run_simulation. -// # Changes: Removed 2 sqrt() calls from the loop bodies -// # Predecessor: variant_inlining_6.cpp - -#include "../headers/board.hpp" -#include "../headers/constants.hpp" -#include "../headers/pedestrian_reduced.hpp" - -#include -#include -#include -#include -#include -#include -#include -#include - -#define SEED 42 -#define BOARD_WIDTH 10.0 -#define BOARD_LENGTH 50.0 -#define BOARD_LENGTH_MINUS_ONE BOARD_LENGTH - 1.0 -#define BOARD_LENGTH_MINUS_ONE_POINT_FIVE BOARD_LENGTH - 1.5 -#define WALL_EPS 1e-3 -#define BOARD_WIDTH_MINUS_WALL_EPS 9.999 - -#ifndef RANDOM_SEED -#define RANDOM_SEED 0 -#endif - -#ifndef BENCHMARK -#define BENCHMARK 0 -#endif - -using namespace std; - -Pedestrian *init_pedestrians(int n_p) { -#if RANDOM_SEED - random_device rd; - mt19937 gen(rd()); -#else - mt19937 gen(SEED); -#endif - - normal_distribution velocity_gen(v0_mean, v0_std_dev); - uniform_real_distribution y_gen(0.3, BOARD_WIDTH - 0.3); - uniform_real_distribution x_gen_left(2.0, BOARD_LENGTH / 2.0); - uniform_real_distribution x_gen_right(BOARD_LENGTH / 2.0, BOARD_LENGTH - 2.0); - - Pedestrian *P = new Pedestrian[n_p]; - - int n_p_half = n_p / 2; - - for (int i = 0; i < n_p_half; i++) { - double x = x_gen_left(gen); - double y = y_gen(gen); - - double start_pos_x = x; - double start_pos_y = y; - - double destination_x = BOARD_LENGTH + 30.0; - double destination_y = y; - - double init_e_x = destination_x - start_pos_x; - double init_e_y = destination_y - start_pos_y; - double init_e_norm = sqrt(init_e_x * init_e_x + init_e_y * init_e_y); - double init_e_unit_x = init_e_x / init_e_norm; - double init_e_unit_y = init_e_y / init_e_norm; - - double desired_v0; - do { - desired_v0 = velocity_gen(gen); - } while (desired_v0 <= 0.1); - - double v_max = 1.3 * desired_v0; - - P[i] = Pedestrian(0, 0, 0, start_pos_x, start_pos_y, init_e_unit_x, init_e_unit_y, v_max, desired_v0, destination_x, destination_y); - } - - for (int i = n_p_half; i < n_p; i++) { - double x = x_gen_right(gen); - double y = y_gen(gen); - - double start_pos_x = x; - double start_pos_y = y; - double destination_x = -30; - double destination_y = y; - - double init_e_x = destination_x - start_pos_x; - double init_e_y = destination_y - start_pos_y; - double init_e_norm = sqrt(init_e_x * init_e_x + init_e_y * init_e_y); - double init_e_unit_x = init_e_x / init_e_norm; - double init_e_unit_y = init_e_y / init_e_norm; - - double desired_v0; - do { - desired_v0 = velocity_gen(gen); - } while (desired_v0 <= 0.1); - - double v_max = 1.3 * desired_v0; - - P[i] = Pedestrian(1, 0, 0, start_pos_x, start_pos_y, init_e_unit_x, init_e_unit_y, v_max, desired_v0, destination_x, destination_y); - } - - return P; -} - -Board init_board_without_obstacles() { return Board(BOARD_LENGTH, BOARD_WIDTH, 0, nullptr); } - -int run_simulation(Pedestrian *pedestrians, Board &board, int n_p, int n_iter, const char *dump_path) { - -#if BENCHMARK - (void)dump_path; -#endif - -#if !BENCHMARK - const int dump_interval = 1; - - vector> samples; - - if (dump_path) { - samples.reserve(static_cast(n_iter / dump_interval + 1) * n_p); - - for (int i = 0; i < n_p; i++) { - Pedestrian &p = pedestrians[i]; - samples.emplace_back(0, i, p.curr_pos_x, p.curr_pos_y, p.e_x, p.e_y, p.color); - } - } -#endif - - // # ADDs: n_iter * n_p * 11 + - // # n_iter * n_p * (n_p - 1) * 14 + - // # n_iter * n_p * board.n_obstacles * 12 - - // # MULTs: n_iter * n_p * 21 + - // # n_iter * n_p * (n_p - 1) * 24 + - // # n_iter * n_p * board.n_obstacles * 19 - - // # DIVs: n_iter * n_p * 2 + - // # n_iter * n_p * (n_p - 1) + - // # n_iter * n_p * board.n_obstacles - - // # SQRTs: n_iter * n_p * 3 + - // # n_iter * n_p * (n_p - 1) * 3 + - // # n_iter * n_p * board.n_obstacles - - // # EXPs: n_iter * n_p * (n_p - 1) + - // # n_iter * n_p * board.n_obstacles - - for (int t = 0; t < n_iter; t++) { - for (int i = 0; i < n_p; i++) { - Pedestrian &p = pedestrians[i]; - - // # 2 MULTs - double desired_velocity_x = p.e_x * p.desired_v0; - double desired_velocity_y = p.e_y * p.desired_v0; - - // # 2 MULTs, 2 ADDs - double f_alpha_dt_x = (desired_velocity_x - p.curr_vel_x) * dt_times_tau_alpha_inv; - double f_alpha_dt_y = (desired_velocity_y - p.curr_vel_y) * dt_times_tau_alpha_inv; - - for (int j = 0; j < n_p; j++) { - if (i == j) - continue; - - // # 2 ADDs - double r_ab_x = p.curr_pos_x - pedestrians[j].curr_pos_x; - double r_ab_y = p.curr_pos_y - pedestrians[j].curr_pos_y; - - // # 1 SQRT, 2 MULTs, 1 ADD - double r_ab_norm = sqrt(r_ab_x * r_ab_x + r_ab_y * r_ab_y); - - // # 2 ADDs - double r_ab2_x = r_ab_x - pedestrians[j].s_b_e_x; - double r_ab2_y = r_ab_y - pedestrians[j].s_b_e_y; - - // # 1 SQRT, 2 MULTs, 1 ADD - double r_ab2_norm = sqrt(r_ab2_x * r_ab2_x + r_ab2_y * r_ab2_y); - - // # 1ADD - double r_ab_norm_sum = r_ab_norm + r_ab2_norm; - - // # 1 MULT, 1 ADD - double two_b_sq = r_ab_norm_sum * r_ab_norm_sum - pedestrians[j].s_b_sq; - - // # 1 SQRT - double two_b = sqrt(two_b_sq); - - // # 5 MULTs, 1 EXP, 1 DIV - double scalar_mult_common = (r_ab_norm_sum / (two_b * r_ab_norm * r_ab2_norm)) * (NEG_V0_AB_div_TWO_SIGMA * exp(two_b * NEG_TWO_SIGMA_INV)); - - // # 2 MULTs - double scalar_mult_r_ab = scalar_mult_common * r_ab2_norm; - double scalar_mult_r_ab2 = scalar_mult_common * r_ab_norm; - - // # 4 MULTs, 2 ADDs - double f_alpha_beta_neg_x = r_ab_x * scalar_mult_r_ab + r_ab2_x * scalar_mult_r_ab2; - double f_alpha_beta_neg_y = r_ab_y * scalar_mult_r_ab + r_ab2_y * scalar_mult_r_ab2; - - double w_val; - - // # 3 MULTs, 1 ADD - double a = (p.e_x * f_alpha_beta_neg_x + p.e_y * f_alpha_beta_neg_y) * cos_phi_inv; - - // # 2 MULTs, 1 ADD - double b_sq = f_alpha_beta_neg_x * f_alpha_beta_neg_x + f_alpha_beta_neg_y * f_alpha_beta_neg_y; - - w_val = dt_neg; - - // # 1 MULT - if (a >= 0 && (a * a > b_sq)) { - w_val = c_neg_times_dt; - } - - // # 2 MULTs, 2 ADDs - f_alpha_dt_x = f_alpha_dt_x + f_alpha_beta_neg_x * w_val; - f_alpha_dt_y = f_alpha_dt_y + f_alpha_beta_neg_y * w_val; - } - - for (int k = 0; k < board.n_obstacles; k++) { - - // # 2 ADDs - double v1_x = p.curr_pos_x - board.obstacles[k].start.x; - double v1_y = p.curr_pos_y - board.obstacles[k].start.y; - - // # 3 MULTs, 1 ADD - double projection = (v1_x * board.obstacles[k].obstacle_v.x + v1_y * board.obstacles[k].obstacle_v.y) * board.obstacles[k].length_sq_inv; - - if (projection < 0.0) - projection = 0.0; - else if (projection > 1.0) - projection = 1.0; - - // # 4 ADDs, 2 MULTs - double r_aB_x = p.curr_pos_x - (board.obstacles[k].start.x + board.obstacles[k].obstacle_v.x * projection); - double r_aB_y = p.curr_pos_y - (board.obstacles[k].start.y + board.obstacles[k].obstacle_v.y * projection); - - // # 1 SQRT, 2 MULTs, 1 ADD - double r_aB_norm = sqrt(r_aB_x * r_aB_x + r_aB_y * r_aB_y); - - // # 4 MULTs, 1 EXP, 1 DIV - double scalar_mult = (NEG_U0_AB_div_R * exp(r_aB_norm * NEG_R_INV)) / r_aB_norm; - double f_alpha_B_neg_x = r_aB_x * scalar_mult; - double f_alpha_B_neg_y = r_aB_y * scalar_mult; - - double w_val; - - // # 3 MULTs, 1 ADD - double a = (p.e_x * f_alpha_B_neg_x + p.e_y * f_alpha_B_neg_y) * cos_phi_inv; - - // # 2 MULTS, 1 ADD - double b_sq = f_alpha_B_neg_x * f_alpha_B_neg_x + f_alpha_B_neg_y * f_alpha_B_neg_y; - - w_val = dt_neg; - - // # 1 MULT - if (a >= 0 && (a * a > b_sq)) { - w_val = c_neg_times_dt; - } - - // # 2 MULTs - double f_alpha_B_dt_x = f_alpha_B_neg_x * w_val; - double f_alpha_B_dt_y = f_alpha_B_neg_y * w_val; - - // # 2 ADDs - f_alpha_dt_x += f_alpha_B_dt_x; - f_alpha_dt_y += f_alpha_B_dt_y; - } - - // # 2 ADDs - p.w_x += f_alpha_dt_x; - p.w_y += f_alpha_dt_y; - } - - for (int i = 0; i < n_p; i++) { - Pedestrian &p = pedestrians[i]; - - // # 1 SQRT, 2 MULTs, 1 ADD - double w_norm = sqrt(p.w_x * p.w_x + p.w_y * p.w_y); - - // # 1 DIV - double correction_factor = min(1.0, p.v_max / w_norm); - - // # 2 MULTs - p.curr_vel_x = p.w_x * correction_factor; - p.curr_vel_y = p.w_y * correction_factor; - - // # 2 MULTs, 2 ADDs - p.curr_pos_x = p.curr_pos_x + p.curr_vel_x * dt; - p.curr_pos_y = p.curr_pos_y + p.curr_vel_y * dt; - - if (p.color == 0 && p.curr_pos_x > BOARD_LENGTH_MINUS_ONE_POINT_FIVE) { - p.curr_pos_x = 1.0; - p.destination_y = p.curr_pos_y; - p.w_x = p.desired_v0; - p.w_y = 0.0; - p.curr_vel_x = p.w_x; - p.curr_vel_y = p.w_y; - } - if (p.color == 1 && p.curr_pos_x < 1.5) { - p.curr_pos_x = BOARD_LENGTH_MINUS_ONE; - p.destination_y = p.curr_pos_y; - p.w_x = -p.desired_v0; - p.w_y = 0.0; - p.curr_vel_x = p.w_x; - p.curr_vel_y = p.w_y; - } - - p.curr_pos_y = clamp(p.curr_pos_y, WALL_EPS, BOARD_WIDTH_MINUS_WALL_EPS); - - // # 2 ADDs - double e_x = p.destination_x - p.curr_pos_x; - double e_y = p.destination_y - p.curr_pos_y; - - // # 1 DIV, 1 SQRT, 2 MULTs, 1 ADD - double e_norm_inv = 1.0 / sqrt(e_x * e_x + e_y * e_y); - - // # 2 MULTs - p.e_x = e_x * e_norm_inv; - p.e_y = e_y * e_norm_inv; - - // # 2 MULTs, 1 ADD - double velocity_norm_sq = p.curr_vel_x * p.curr_vel_x + p.curr_vel_y * p.curr_vel_y; - - // # 1 SQRT - double velocity_norm = sqrt(velocity_norm_sq); - - // # 5 MULTs - p.s_b_sq = velocity_norm_sq * delta_t_sq; - p.s_b_e_x = p.e_x * (velocity_norm * delta_t); - p.s_b_e_y = p.e_y * (velocity_norm * delta_t); - -#if !BENCHMARK - if (dump_path && (t % dump_interval == 0)) { - samples.emplace_back(t + 1, i, p.curr_pos_x, p.curr_pos_y, p.e_x, p.e_y, p.color); - } -#endif - } - } - -#if !BENCHMARK - if (dump_path) { - ofstream out(dump_path); - out << "time,id,x,y,p.e.x,p.e.y,color\n"; - - for (auto &[t, id, x, y, e_x, e_y, c] : samples) { - out << t << ',' << id << ',' << x << ',' << y << ',' << e_x << ',' << e_y << ',' << c << '\n'; - } - } -#endif - - return 0; -} diff --git a/variants/variant_inlining_8_flopcount.cpp b/variants/variant_inlining_8_flopcount.cpp deleted file mode 100644 index 45ffc24..0000000 --- a/variants/variant_inlining_8_flopcount.cpp +++ /dev/null @@ -1,343 +0,0 @@ -// # This file performs in-lining for run_simulation. -// # Changes: Replaced references to Vector2d structs inside each Obstacle of Board.obstacles with references to doubles -// # Predecessor: variant_inlining_7.cpp - -#include "../headers/board_reduced.hpp" -#include "../headers/constants.hpp" -#include "../headers/pedestrian_reduced.hpp" - -#include -#include -#include -#include -#include -#include -#include -#include - -#define SEED 42 -#define BOARD_WIDTH 10.0 -#define BOARD_LENGTH 50.0 -#define BOARD_LENGTH_MINUS_ONE BOARD_LENGTH - 1.0 -#define BOARD_LENGTH_MINUS_ONE_POINT_FIVE BOARD_LENGTH - 1.5 -#define WALL_EPS 1e-3 -#define BOARD_WIDTH_MINUS_WALL_EPS 9.999 - -#ifndef RANDOM_SEED -#define RANDOM_SEED 0 -#endif - -#ifndef BENCHMARK -#define BENCHMARK 0 -#endif - -using namespace std; - -Pedestrian *init_pedestrians(int n_p) { -#if RANDOM_SEED - random_device rd; - mt19937 gen(rd()); -#else - mt19937 gen(SEED); -#endif - - normal_distribution velocity_gen(v0_mean, v0_std_dev); - uniform_real_distribution y_gen(0.3, BOARD_WIDTH - 0.3); - uniform_real_distribution x_gen_left(2.0, BOARD_LENGTH / 2.0); - uniform_real_distribution x_gen_right(BOARD_LENGTH / 2.0, BOARD_LENGTH - 2.0); - - Pedestrian *P = new Pedestrian[n_p]; - - int n_p_half = n_p / 2; - - for (int i = 0; i < n_p_half; i++) { - double x = x_gen_left(gen); - double y = y_gen(gen); - - double start_pos_x = x; - double start_pos_y = y; - - double destination_x = BOARD_LENGTH + 30.0; - double destination_y = y; - - double init_e_x = destination_x - start_pos_x; - double init_e_y = destination_y - start_pos_y; - double init_e_norm = sqrt(init_e_x * init_e_x + init_e_y * init_e_y); - double init_e_unit_x = init_e_x / init_e_norm; - double init_e_unit_y = init_e_y / init_e_norm; - - double desired_v0; - do { - desired_v0 = velocity_gen(gen); - } while (desired_v0 <= 0.1); - - double v_max = 1.3 * desired_v0; - - P[i] = Pedestrian(0, 0, 0, start_pos_x, start_pos_y, init_e_unit_x, init_e_unit_y, v_max, desired_v0, destination_x, destination_y); - } - - for (int i = n_p_half; i < n_p; i++) { - double x = x_gen_right(gen); - double y = y_gen(gen); - - double start_pos_x = x; - double start_pos_y = y; - double destination_x = -30; - double destination_y = y; - - double init_e_x = destination_x - start_pos_x; - double init_e_y = destination_y - start_pos_y; - double init_e_norm = sqrt(init_e_x * init_e_x + init_e_y * init_e_y); - double init_e_unit_x = init_e_x / init_e_norm; - double init_e_unit_y = init_e_y / init_e_norm; - - double desired_v0; - do { - desired_v0 = velocity_gen(gen); - } while (desired_v0 <= 0.1); - - double v_max = 1.3 * desired_v0; - - P[i] = Pedestrian(1, 0, 0, start_pos_x, start_pos_y, init_e_unit_x, init_e_unit_y, v_max, desired_v0, destination_x, destination_y); - } - - return P; -} - -Board init_board_without_obstacles() { return Board(BOARD_LENGTH, BOARD_WIDTH, 0, nullptr); } - -int run_simulation(Pedestrian *pedestrians, Board &board, int n_p, int n_iter, const char *dump_path) { - -#if BENCHMARK - (void)dump_path; -#endif - -#if !BENCHMARK - const int dump_interval = 1; - - vector> samples; - - if (dump_path) { - samples.reserve(static_cast(n_iter / dump_interval + 1) * n_p); - - for (int i = 0; i < n_p; i++) { - Pedestrian &p = pedestrians[i]; - samples.emplace_back(0, i, p.curr_pos_x, p.curr_pos_y, p.e_x, p.e_y, p.color); - } - } -#endif - - // # ADDs: n_iter * n_p * 11 + - // # n_iter * n_p * (n_p - 1) * 14 + - // # n_iter * n_p * board.n_obstacles * 10 - - // # MULTs: n_iter * n_p * 21 + - // # n_iter * n_p * (n_p - 1) * 24 + - // # n_iter * n_p * board.n_obstacles * 13 - - // # DIVs: n_iter * n_p * 2 + - // # n_iter * n_p * (n_p - 1) + - // # n_iter * n_p * board.n_obstacles - - // # SQRTs: n_iter * n_p * 3 + - // # n_iter * n_p * (n_p - 1) * 3 + - // # n_iter * n_p * board.n_obstacles - - // # EXPs: n_iter * n_p * (n_p - 1) + - // # n_iter * n_p * board.n_obstacles - - for (int t = 0; t < n_iter; t++) { - for (int i = 0; i < n_p; i++) { - Pedestrian &p = pedestrians[i]; - - // # 2 MULTs - double desired_velocity_x = p.e_x * p.desired_v0; - double desired_velocity_y = p.e_y * p.desired_v0; - - // # 2 MULTs, 2 ADDs - double f_alpha_dt_x = (desired_velocity_x - p.curr_vel_x) * dt_times_tau_alpha_inv; - double f_alpha_dt_y = (desired_velocity_y - p.curr_vel_y) * dt_times_tau_alpha_inv; - - for (int j = 0; j < n_p; j++) { - if (i == j) - continue; - - // # 2 ADDs - double r_ab_x = p.curr_pos_x - pedestrians[j].curr_pos_x; - double r_ab_y = p.curr_pos_y - pedestrians[j].curr_pos_y; - - // # 1 SQRT, 2 MULTs, 1 ADD - double r_ab_norm = sqrt(r_ab_x * r_ab_x + r_ab_y * r_ab_y); - - // # 2 ADDs - double r_ab2_x = r_ab_x - pedestrians[j].s_b_e_x; - double r_ab2_y = r_ab_y - pedestrians[j].s_b_e_y; - - // # 1 SQRT, 2 MULTs, 1 ADD - double r_ab2_norm = sqrt(r_ab2_x * r_ab2_x + r_ab2_y * r_ab2_y); - - // # 1ADD - double r_ab_norm_sum = r_ab_norm + r_ab2_norm; - - // # 1 MULT, 1 ADD - double two_b_sq = r_ab_norm_sum * r_ab_norm_sum - pedestrians[j].s_b_sq; - - // # 1 SQRT - double two_b = sqrt(two_b_sq); - - // # 5 MULTs, 1 EXP, 1 DIV - double scalar_mult_common = (r_ab_norm_sum / (two_b * r_ab_norm * r_ab2_norm)) * (NEG_V0_AB_div_TWO_SIGMA * exp(two_b * NEG_TWO_SIGMA_INV)); - - // # 2 MULTs - double scalar_mult_r_ab = scalar_mult_common * r_ab2_norm; - double scalar_mult_r_ab2 = scalar_mult_common * r_ab_norm; - - // # 4 MULTs, 2 ADDs - double f_alpha_beta_neg_x = r_ab_x * scalar_mult_r_ab + r_ab2_x * scalar_mult_r_ab2; - double f_alpha_beta_neg_y = r_ab_y * scalar_mult_r_ab + r_ab2_y * scalar_mult_r_ab2; - - double w_val; - - // # 3 MULTs, 1 ADD - double a = (p.e_x * f_alpha_beta_neg_x + p.e_y * f_alpha_beta_neg_y) * cos_phi_inv; - - // # 2 MULTs, 1 ADD - double b_sq = f_alpha_beta_neg_x * f_alpha_beta_neg_x + f_alpha_beta_neg_y * f_alpha_beta_neg_y; - - w_val = dt_neg; - - // # 1 MULT - if (a >= 0 && (a * a > b_sq)) { - w_val = c_neg_times_dt; - } - - // # 2 MULTs, 2 ADDs - f_alpha_dt_x = f_alpha_dt_x + f_alpha_beta_neg_x * w_val; - f_alpha_dt_y = f_alpha_dt_y + f_alpha_beta_neg_y * w_val; - } - - for (int k = 0; k < board.n_obstacles; k++) { - - // # 2 ADDs - double v1_x = p.curr_pos_x - board.obstacles[k].start_x; - double v1_y = p.curr_pos_y - board.obstacles[k].start_y; - - // # 3 MULTs, 1 ADD - double projection = (v1_x * board.obstacles[k].obstacle_v_x + v1_y * board.obstacles[k].obstacle_v_y) * board.obstacles[k].length_sq_inv; - - if (projection < 0.0) - projection = 0.0; - else if (projection > 1.0) - projection = 1.0; - - // # 4 ADDs, 2 MULTs - double r_aB_x = p.curr_pos_x - (board.obstacles[k].start_x + board.obstacles[k].obstacle_v_x * projection); - double r_aB_y = p.curr_pos_y - (board.obstacles[k].start_y + board.obstacles[k].obstacle_v_y * projection); - - // # 1 SQRT, 2 MULTs, 1 ADD - double r_aB_norm = sqrt(r_aB_x * r_aB_x + r_aB_y * r_aB_y); - // # 4 MULTs, 1 EXP, 1 DIV - double scalar_mult = (NEG_U0_AB_div_R * exp(r_aB_norm * NEG_R_INV)) / r_aB_norm; - double f_alpha_B_neg_x = r_aB_x * scalar_mult; - double f_alpha_B_neg_y = r_aB_y * scalar_mult; - - double w_val; - double e_dot_f = (f_alpha_B_neg_y == 0) ? 0 : (f_alpha_B_neg_y < 0 ? -p.e_y : p.e_y); - if (e_dot_f >= cos_phi) - w_val = dt_neg; - else - w_val = c_neg_times_dt; - - // # 2 MULTs - double f_alpha_B_dt_x = f_alpha_B_neg_x * w_val; - double f_alpha_B_dt_y = f_alpha_B_neg_y * w_val; - - // # 2 ADDs - f_alpha_dt_x += f_alpha_B_dt_x; - f_alpha_dt_y += f_alpha_B_dt_y; - } - - // # 2 ADDs - p.w_x += f_alpha_dt_x; - p.w_y += f_alpha_dt_y; - } - - for (int i = 0; i < n_p; i++) { - Pedestrian &p = pedestrians[i]; - - // # 1 SQRT, 2 MULTs, 1 ADD - double w_norm = sqrt(p.w_x * p.w_x + p.w_y * p.w_y); - - // # 1 DIV - double correction_factor = min(1.0, p.v_max / w_norm); - - // # 2 MULTs - p.curr_vel_x = p.w_x * correction_factor; - p.curr_vel_y = p.w_y * correction_factor; - - // # 2 MULTs, 2 ADDs - p.curr_pos_x = p.curr_pos_x + p.curr_vel_x * dt; - p.curr_pos_y = p.curr_pos_y + p.curr_vel_y * dt; - - if (p.color == 0 && p.curr_pos_x > BOARD_LENGTH_MINUS_ONE_POINT_FIVE) { - p.curr_pos_x = 1.0; - p.destination_y = p.curr_pos_y; - p.w_x = p.desired_v0; - p.w_y = 0.0; - p.curr_vel_x = p.w_x; - p.curr_vel_y = p.w_y; - } - if (p.color == 1 && p.curr_pos_x < 1.5) { - p.curr_pos_x = BOARD_LENGTH_MINUS_ONE; - p.destination_y = p.curr_pos_y; - p.w_x = -p.desired_v0; - p.w_y = 0.0; - p.curr_vel_x = p.w_x; - p.curr_vel_y = p.w_y; - } - - p.curr_pos_y = clamp(p.curr_pos_y, WALL_EPS, BOARD_WIDTH_MINUS_WALL_EPS); - - // # 2 ADDs - double e_x = p.destination_x - p.curr_pos_x; - double e_y = p.destination_y - p.curr_pos_y; - - // # 1 DIV, 1 SQRT, 2 MULTs, 1 ADD - double e_norm_inv = 1.0 / sqrt(e_x * e_x + e_y * e_y); - - // # 2 MULTs - p.e_x = e_x * e_norm_inv; - p.e_y = e_y * e_norm_inv; - - // # 2 MULTs, 1 ADD - double velocity_norm_sq = p.curr_vel_x * p.curr_vel_x + p.curr_vel_y * p.curr_vel_y; - - // # 1 SQRT - double velocity_norm = sqrt(velocity_norm_sq); - - // # 5 MULTs - p.s_b_sq = velocity_norm_sq * delta_t_sq; - p.s_b_e_x = p.e_x * (velocity_norm * delta_t); - p.s_b_e_y = p.e_y * (velocity_norm * delta_t); - -#if !BENCHMARK - if (dump_path && (t % dump_interval == 0)) { - samples.emplace_back(t + 1, i, p.curr_pos_x, p.curr_pos_y, p.e_x, p.e_y, p.color); - } -#endif - } - } - -#if !BENCHMARK - if (dump_path) { - ofstream out(dump_path); - out << "time,id,x,y,p.e.x,p.e.y,color\n"; - - for (auto &[t, id, x, y, e_x, e_y, c] : samples) { - out << t << ',' << id << ',' << x << ',' << y << ',' << e_x << ',' << e_y << ',' << c << '\n'; - } - } -#endif - - return 0; -}