From 402e72a2484e24b64b7b82130a6d323fa85b905b Mon Sep 17 00:00:00 2001 From: Mark Date: Thu, 5 Sep 2019 20:34:57 -0400 Subject: [PATCH 1/8] finish part 1 --- src/kernel.cu | 67 ++++++++++++++++++++++++++++++++++++++++++++++++--- src/main.cpp | 2 +- 2 files changed, 64 insertions(+), 5 deletions(-) diff --git a/src/kernel.cu b/src/kernel.cu index 74dffcb..caa7f1a 100644 --- a/src/kernel.cu +++ b/src/kernel.cu @@ -137,6 +137,7 @@ __global__ void kernGenerateRandomPosArray(int time, int N, glm::vec3 * arr, flo * Initialize memory, update some globals */ void Boids::initSimulation(int N) { + //initialize numObjects here numObjects = N; dim3 fullBlocksPerGrid((N + blockSize - 1) / blockSize); @@ -231,9 +232,44 @@ void Boids::copyBoidsToVBO(float *vbodptr_positions, float *vbodptr_velocities) */ __device__ glm::vec3 computeVelocityChange(int N, int iSelf, const glm::vec3 *pos, const glm::vec3 *vel) { // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves - // Rule 2: boids try to stay a distance d away from each other - // Rule 3: boids try to match the speed of surrounding boids - return glm::vec3(0.0f, 0.0f, 0.0f); + glm::vec3 perceived_center; + glm::vec3 c; + glm::vec3 perceived_velocity; + glm::vec3 return_vel; + float neighbor_count_rule1 = 0; + float neighbor_count_rule3 = 0; + for (int idx = 0; idx < N; ++idx) + { + //if b = boid skip the rest actions + if (idx == iSelf) continue; + + // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves + if (glm::distance(pos[idx], pos[iSelf]) < rule1Distance) + { + perceived_center += pos[idx]; + neighbor_count_rule1++; + } + + // Rule 2: boids try to stay a distance d away from each other + if (glm::distance(pos[idx], pos[iSelf]) < rule2Distance) + { + c += (pos[idx] - pos[iSelf]); + } + + // Rule 3: boids try to match the speed of surrounding boids + if (glm::distance(pos[idx], pos[iSelf]) < rule3Distance) + { + perceived_velocity += vel[idx]; + neighbor_count_rule3++; + } + } + + perceived_center /= neighbor_count_rule1; + perceived_velocity /= neighbor_count_rule3; + + return_vel += ((perceived_center - pos[iSelf]) * rule1Scale) + (c * rule2Scale) + (perceived_center * rule3Scale); + + return return_vel; } /** @@ -243,8 +279,20 @@ __device__ glm::vec3 computeVelocityChange(int N, int iSelf, const glm::vec3 *po __global__ void kernUpdateVelocityBruteForce(int N, glm::vec3 *pos, glm::vec3 *vel1, glm::vec3 *vel2) { // Compute a new velocity based on pos and vel1 + //Compute the index of current thread + int index = threadIdx.x + (blockIdx.x * blockDim.x); + glm::vec3 curr_velocity = computeVelocityChange(N, index, pos, vel1); // Clamp the speed - // Record the new velocity into vel2. Question: why NOT vel1? + glm::vec3 zero_vec; + float curr_speed = glm::distance(curr_velocity, zero_vec); + //if the total speed of vel is larger than maxSpeed, we normalize the vel and apply the maxSpeed we allow -- do we need to care negative speed? + if (curr_speed > maxSpeed) + { + curr_velocity = glm::normalize(curr_velocity) * maxSpeed; + } + + // Record the new velocity into vel2. Question: why NOT vel1? --- because other boids might need that + vel2[index] = curr_velocity; } /** @@ -348,7 +396,18 @@ __global__ void kernUpdateVelNeighborSearchCoherent( */ void Boids::stepSimulationNaive(float dt) { // TODO-1.2 - use the kernels you wrote to step the simulation forward in time. + int gridSize = numObjects / blockSize; + dim3 blocksPerGrid(gridSize); + //first compute the new velocity + kernUpdateVelocityBruteForce <<< blocksPerGrid, threadsPerBlock >>> (numObjects, dev_pos, dev_vel1, dev_vel2); + + //Then update the pos + kernUpdatePos << < blocksPerGrid, threadsPerBlock >> > (numObjects, dt, dev_pos, dev_vel2); + // TODO-1.2 ping-pong the velocity buffers + glm::vec3 *dev_temp = dev_vel1; + dev_vel1 = dev_vel2; + dev_vel2 = dev_temp; } void Boids::stepSimulationScatteredGrid(float dt) { diff --git a/src/main.cpp b/src/main.cpp index b82c8c6..a7b8586 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -18,7 +18,7 @@ #define COHERENT_GRID 0 // LOOK-1.2 - change this to adjust particle count in the simulation -const int N_FOR_VIS = 5000; +const int N_FOR_VIS = 10000; const float DT = 0.2f; /** From dbd869938ddb8a03fed4f93a59b310dd66fffffc Mon Sep 17 00:00:00 2001 From: Mark Date: Thu, 5 Sep 2019 21:47:48 -0400 Subject: [PATCH 2/8] fixed bugs in part 1 --- src/kernel.cu | 44 +++++++++++++++++++++++++++----------------- 1 file changed, 27 insertions(+), 17 deletions(-) diff --git a/src/kernel.cu b/src/kernel.cu index caa7f1a..a2339dc 100644 --- a/src/kernel.cu +++ b/src/kernel.cu @@ -232,42 +232,55 @@ void Boids::copyBoidsToVBO(float *vbodptr_positions, float *vbodptr_velocities) */ __device__ glm::vec3 computeVelocityChange(int N, int iSelf, const glm::vec3 *pos, const glm::vec3 *vel) { // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves - glm::vec3 perceived_center; - glm::vec3 c; - glm::vec3 perceived_velocity; - glm::vec3 return_vel; + glm::vec3 perceived_center(0.f,0.f,0.f); + glm::vec3 avoidance_velocity(0.f, 0.f, 0.f); + glm::vec3 perceived_velocity(0.f, 0.f, 0.f); + glm::vec3 return_vel(0.f, 0.f, 0.f); float neighbor_count_rule1 = 0; float neighbor_count_rule3 = 0; + + //pre load all needed data + glm::vec3 curr_boid_pos = pos[iSelf]; for (int idx = 0; idx < N; ++idx) { //if b = boid skip the rest actions if (idx == iSelf) continue; + //load current boid pos + glm::vec3 idx_boid_pos = pos[idx]; + float dist = glm::distance(idx_boid_pos, curr_boid_pos); // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves - if (glm::distance(pos[idx], pos[iSelf]) < rule1Distance) + if (dist < rule1Distance) { - perceived_center += pos[idx]; + perceived_center += idx_boid_pos; neighbor_count_rule1++; } // Rule 2: boids try to stay a distance d away from each other - if (glm::distance(pos[idx], pos[iSelf]) < rule2Distance) + if (dist < rule2Distance) { - c += (pos[idx] - pos[iSelf]); + avoidance_velocity -= (idx_boid_pos - curr_boid_pos); } // Rule 3: boids try to match the speed of surrounding boids - if (glm::distance(pos[idx], pos[iSelf]) < rule3Distance) + if (dist < rule3Distance) { perceived_velocity += vel[idx]; neighbor_count_rule3++; } } + //if we use N-1, the particles will shrink to the center of cube -- safer + //glm::vec3 rule1_component = neighbor_count_rule1 > 0 ? (perceived_center / neighbor_count_rule1 - curr_boid_pos) * rule1Scale : glm::vec3(0.f, 0.f, 0.f); + //glm::vec3 rule2_component = avoidance_velocity * rule2Scale; + //glm::vec3 rule3_component = neighbor_count_rule3 > 0 ? (perceived_velocity / neighbor_count_rule3) * rule3Scale : glm::vec3(0.f, 0.f, 0.f); + + //if we use N-1, the particles will shrink to the center of cube -- similar if we set perceived_center/velocity = 0 when neighbor_count is 0 perceived_center /= neighbor_count_rule1; perceived_velocity /= neighbor_count_rule3; - return_vel += ((perceived_center - pos[iSelf]) * rule1Scale) + (c * rule2Scale) + (perceived_center * rule3Scale); + //helped by Hanna ReadMe Rule part sum all rules' and current velocity + return_vel += vel[iSelf] + ((perceived_center - curr_boid_pos) * rule1Scale) + (avoidance_velocity * rule2Scale) + (perceived_velocity * rule3Scale); return return_vel; } @@ -283,8 +296,7 @@ __global__ void kernUpdateVelocityBruteForce(int N, glm::vec3 *pos, int index = threadIdx.x + (blockIdx.x * blockDim.x); glm::vec3 curr_velocity = computeVelocityChange(N, index, pos, vel1); // Clamp the speed - glm::vec3 zero_vec; - float curr_speed = glm::distance(curr_velocity, zero_vec); + float curr_speed = glm::length(curr_velocity); //if the total speed of vel is larger than maxSpeed, we normalize the vel and apply the maxSpeed we allow -- do we need to care negative speed? if (curr_speed > maxSpeed) { @@ -396,7 +408,7 @@ __global__ void kernUpdateVelNeighborSearchCoherent( */ void Boids::stepSimulationNaive(float dt) { // TODO-1.2 - use the kernels you wrote to step the simulation forward in time. - int gridSize = numObjects / blockSize; + int gridSize = (numObjects + blockSize - 1) / blockSize; //helped by Gangzheng Tong dim3 blocksPerGrid(gridSize); //first compute the new velocity kernUpdateVelocityBruteForce <<< blocksPerGrid, threadsPerBlock >>> (numObjects, dev_pos, dev_vel1, dev_vel2); @@ -404,10 +416,8 @@ void Boids::stepSimulationNaive(float dt) { //Then update the pos kernUpdatePos << < blocksPerGrid, threadsPerBlock >> > (numObjects, dt, dev_pos, dev_vel2); - // TODO-1.2 ping-pong the velocity buffers - glm::vec3 *dev_temp = dev_vel1; - dev_vel1 = dev_vel2; - dev_vel2 = dev_temp; + // TODO-1.2 ping-pong the velocity buffers -- swap content + cudaMemcpy(dev_vel1, dev_vel2, sizeof(glm::vec3) * numObjects, cudaMemcpyDeviceToDevice); } void Boids::stepSimulationScatteredGrid(float dt) { From 213de2915c5cfeec36e45c64455bf40cdbbc27bd Mon Sep 17 00:00:00 2001 From: Mark Date: Fri, 6 Sep 2019 15:32:14 -0400 Subject: [PATCH 3/8] finish part 2.1 --- src/kernel.cu | 387 ++++++++++++++++++++++++++++++++++++++++---------- src/main.cpp | 4 +- 2 files changed, 311 insertions(+), 80 deletions(-) diff --git a/src/kernel.cu b/src/kernel.cu index a2339dc..ed9a521 100644 --- a/src/kernel.cu +++ b/src/kernel.cu @@ -17,6 +17,8 @@ #define checkCUDAErrorWithLine(msg) checkCUDAError(msg, __LINE__) +#define CELL_WIDTH_FACTOR 2.0f + /** * Check for CUDA errors; print and exit if there was a problem. */ @@ -158,7 +160,7 @@ void Boids::initSimulation(int N) { checkCUDAErrorWithLine("kernGenerateRandomPosArray failed!"); // LOOK-2.1 computing grid params - gridCellWidth = 2.0f * std::max(std::max(rule1Distance, rule2Distance), rule3Distance); + gridCellWidth = CELL_WIDTH_FACTOR * std::max(std::max(rule1Distance, rule2Distance), rule3Distance); int halfSideCount = (int)(scene_scale / gridCellWidth) + 1; gridSideCount = 2 * halfSideCount; @@ -170,6 +172,22 @@ void Boids::initSimulation(int N) { gridMinimum.z -= halfGridWidth; // TODO-2.1 TODO-2.3 - Allocate additional buffers here. + cudaMalloc((void**)&dev_particleArrayIndices, N * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc dev_particleArrayIndices failed!"); + + cudaMalloc((void**)&dev_particleGridIndices, N * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc dev_particleGridIndices failed!"); + + //these two arrays only use to store cell information, so no need to be as many as numObjects + cudaMalloc((void**)&dev_gridCellStartIndices, gridCellCount * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc dev_gridCellStartIndices failed!"); + + cudaMalloc((void**)&dev_gridCellEndIndices, gridCellCount * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc dev_gridCellEndIndices failed!"); + + dev_thrust_particleArrayIndices = thrust::device_ptr(dev_particleArrayIndices); + dev_thrust_particleGridIndices = thrust::device_ptr(dev_particleGridIndices); + cudaDeviceSynchronize(); } @@ -270,17 +288,22 @@ __device__ glm::vec3 computeVelocityChange(int N, int iSelf, const glm::vec3 *po } } - //if we use N-1, the particles will shrink to the center of cube -- safer - //glm::vec3 rule1_component = neighbor_count_rule1 > 0 ? (perceived_center / neighbor_count_rule1 - curr_boid_pos) * rule1Scale : glm::vec3(0.f, 0.f, 0.f); - //glm::vec3 rule2_component = avoidance_velocity * rule2Scale; - //glm::vec3 rule3_component = neighbor_count_rule3 > 0 ? (perceived_velocity / neighbor_count_rule3) * rule3Scale : glm::vec3(0.f, 0.f, 0.f); - - //if we use N-1, the particles will shrink to the center of cube -- similar if we set perceived_center/velocity = 0 when neighbor_count is 0 - perceived_center /= neighbor_count_rule1; - perceived_velocity /= neighbor_count_rule3; + //if we use N-1, the particles will shrink to the center of cube -- helped by Hannar + glm::vec3 rule1_component = glm::vec3(0.f, 0.f, 0.f); + glm::vec3 rule2_component = glm::vec3(0.f, 0.f, 0.f); + glm::vec3 rule3_component = glm::vec3(0.f, 0.f, 0.f); + if (neighbor_count_rule1 > 0) + { + rule1_component = (perceived_center / neighbor_count_rule1 - curr_boid_pos) * rule1Scale; + } + rule2_component = avoidance_velocity * rule2Scale; + if (neighbor_count_rule3 > 0) + { + rule3_component = (perceived_velocity / neighbor_count_rule3) * rule3Scale; + } //helped by Hanna ReadMe Rule part sum all rules' and current velocity - return_vel += vel[iSelf] + ((perceived_center - curr_boid_pos) * rule1Scale) + (avoidance_velocity * rule2Scale) + (perceived_velocity * rule3Scale); + return_vel += vel[iSelf] + rule1_component + rule2_component + rule3_component; return return_vel; } @@ -294,17 +317,21 @@ __global__ void kernUpdateVelocityBruteForce(int N, glm::vec3 *pos, // Compute a new velocity based on pos and vel1 //Compute the index of current thread int index = threadIdx.x + (blockIdx.x * blockDim.x); - glm::vec3 curr_velocity = computeVelocityChange(N, index, pos, vel1); + if (index > N) + { + return; + } + glm::vec3 new_velocity = computeVelocityChange(N, index, pos, vel1); // Clamp the speed - float curr_speed = glm::length(curr_velocity); + float curr_speed = glm::length(new_velocity); //if the total speed of vel is larger than maxSpeed, we normalize the vel and apply the maxSpeed we allow -- do we need to care negative speed? if (curr_speed > maxSpeed) { - curr_velocity = glm::normalize(curr_velocity) * maxSpeed; + new_velocity = glm::normalize(new_velocity) * maxSpeed; } // Record the new velocity into vel2. Question: why NOT vel1? --- because other boids might need that - vel2[index] = curr_velocity; + vel2[index] = new_velocity; } /** @@ -347,8 +374,25 @@ __global__ void kernComputeIndices(int N, int gridResolution, glm::vec3 *pos, int *indices, int *gridIndices) { // TODO-2.1 // - Label each boid with the index of its grid cell. - // - Set up a parallel array of integer indices as pointers to the actual - // boid data in pos and vel1/vel2 + //compute the correspond index in x, y, z axis and use gridIndex3Dto1D to store the actual 1D index -- why we need inverseCellWidth -- mult is faster than divide + int idx = threadIdx.x + (blockIdx.x * blockDim.x); + if (idx >= N) { + return; + } + //pre store info + glm::vec3 curr_pos = pos[idx]; + int idx_x = (curr_pos.x - gridMin.x) * inverseCellWidth; + int idx_y = (curr_pos.y - gridMin.y) * inverseCellWidth; + int idx_z = (curr_pos.z - gridMin.z) * inverseCellWidth; + + //combine to get the 1D index + int gridIndex = gridIndex3Dto1D(idx_x, idx_y, idx_z, gridResolution); + + //store to indices and gridIndices correspondingly + // - Set up a parallel array of integer indices as pointers to the actual + // boid data in pos and vel1/vel2 + indices[idx] = idx; + gridIndices[idx] = gridIndex; } // LOOK-2.1 Consider how this could be useful for indicating that a cell @@ -366,8 +410,53 @@ __global__ void kernIdentifyCellStartEnd(int N, int *particleGridIndices, // Identify the start point of each cell in the gridIndices array. // This is basically a parallel unrolling of a loop that goes // "this index doesn't match the one before it, must be a new cell!" + //what happen for those we don't include any particle? is there an identifier to show its identity? + + //so why do we need th Particles array? we don't even pass in it -- the arranged start and end indices are for us to access dev_particleArrayIndices, which has been sorted + //may have better way to do + int idx = threadIdx.x + (blockIdx.x * blockDim.x); + if (idx >= N) { + return; + } + int target_grid_index = particleGridIndices[idx]; + //look at the one before and the one after, if diff, store + //head must be a start + if (idx == 0) + { + gridCellStartIndices[target_grid_index] = idx; + } + + if (idx == N - 1) + { + gridCellEndIndices[target_grid_index] = idx; + return; + } + + + //check one before and one after + if (target_grid_index != 0 && target_grid_index != particleGridIndices[idx - 1]) + { + //start of a cell + gridCellStartIndices[target_grid_index] = idx; + + } + + if (target_grid_index != N - 1 && target_grid_index != particleGridIndices[idx + 1]) + { + gridCellEndIndices[target_grid_index] = idx; + + } + + //int next_grid_index = particleGridIndices[idx + 1]; + //if (target_grid_index != next_grid_index) + //{ + // gridCellStartIndices[next_grid_index] = idx + 1; + // gridCellEndIndices[target_grid_index] = idx; + //} + } +//very similar to photon mapping __global__ void kernUpdateVelNeighborSearchScattered( int N, int gridResolution, glm::vec3 gridMin, float inverseCellWidth, float cellWidth, @@ -377,13 +466,124 @@ __global__ void kernUpdateVelNeighborSearchScattered( // TODO-2.1 - Update a boid's velocity using the uniform grid to reduce // the number of boids that need to be checked. // - Identify the grid cell that this particle is in + int idx = threadIdx.x + (blockIdx.x * blockDim.x); + if (idx >= N) { + return; + } + glm::vec3 curr_boid_pos = pos[idx]; + int idx_x = (curr_boid_pos.x - gridMin.x) * inverseCellWidth; + int idx_y = (curr_boid_pos.y - gridMin.y) * inverseCellWidth; + int idx_z = (curr_boid_pos.z - gridMin.z) * inverseCellWidth; + + //combine to get the 1D index + int gridIndex = gridIndex3Dto1D(idx_x, idx_y, idx_z, gridResolution); // - Identify which cells may contain neighbors. This isn't always 8. + //by calculating those cells that interact with the sphere with the neighbor_radius(std::max(std::max(rule1Distance, rule2Distance), rule3Distance)) -not allowed + //float neighbor_radius = cellWidth / CELL_WIDTH_FACTOR; + const float neighbor_radius = glm::max(glm::max(rule1Distance, rule2Distance), rule3Distance); + int max_x = (curr_boid_pos.x + neighbor_radius - gridMin.x) * inverseCellWidth; + max_x = max_x > gridResolution - 1 ? gridResolution - 1 : max_x; + int min_x = (curr_boid_pos.x - neighbor_radius - gridMin.x) * inverseCellWidth; + min_x = min_x < 0 ? 0 : min_x; + int max_y = (curr_boid_pos.y + neighbor_radius - gridMin.y) * inverseCellWidth; + max_y = max_y > gridResolution - 1 ? gridResolution - 1 : max_y; + int min_y = (curr_boid_pos.y + neighbor_radius - gridMin.y) * inverseCellWidth; + min_y = min_y < 0 ? 0 : min_y; + int max_z = (curr_boid_pos.z + neighbor_radius - gridMin.z) * inverseCellWidth; + max_z = max_z > gridResolution - 1 ? gridResolution - 1 : max_z; + int min_z = (curr_boid_pos.z + neighbor_radius - gridMin.z) * inverseCellWidth; + min_z = min_z < 0 ? 0 : min_z; + // - For each cell, read the start/end indices in the boid pointer array. - // - Access each boid in the cell and compute velocity change from - // the boids rules, if this boid is within the neighborhood distance. + glm::vec3 perceived_center(0.f, 0.f, 0.f); + glm::vec3 avoidance_velocity(0.f, 0.f, 0.f); + glm::vec3 perceived_velocity(0.f, 0.f, 0.f); + glm::vec3 new_velocity(0.f, 0.f, 0.f); + float neighbor_count_rule1 = 0; + float neighbor_count_rule3 = 0; + for (int x_cord_idx = min_x; x_cord_idx <= max_x; ++x_cord_idx) + { + for (int y_cord_idx = min_y; y_cord_idx <= max_y; ++y_cord_idx) + { + for (int z_cord_idx = min_z; z_cord_idx <= max_z; ++z_cord_idx) + { + int curr_gridIndex = gridIndex3Dto1D(x_cord_idx, y_cord_idx, z_cord_idx, gridResolution); + //read the start/end indices + int start = gridCellStartIndices[curr_gridIndex]; + int end = gridCellEndIndices[curr_gridIndex]; + if (start == -1 || end == -1) { + continue; //no boid in this cell + } + else + { + // - Access each boid in the cell and compute velocity change from + // the boids rules, if this boid is within the neighborhood distance. + for (int boid_array_idx = start; boid_array_idx <= end; ++boid_array_idx) + { + //the boid_array_idx is only the index in particleArray, need to load it out -- buggy for only few particles moving + int boid_idx = particleArrayIndices[boid_array_idx]; + + if (boid_idx == idx) continue; + + //pre load temp boid pos + glm::vec3 idx_boid_pos = pos[boid_idx]; + float dist = glm::distance(idx_boid_pos, curr_boid_pos); + // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves + if (dist < rule1Distance) + { + perceived_center += idx_boid_pos; + neighbor_count_rule1++; + } + + // Rule 2: boids try to stay a distance d away from each other + if (dist < rule2Distance) + { + avoidance_velocity -= (idx_boid_pos - curr_boid_pos); + } + + // Rule 3: boids try to match the speed of surrounding boids + if (dist < rule3Distance) + { + perceived_velocity += vel1[idx]; + neighbor_count_rule3++; + } + + } + + } + } + + } + + } + + //compute the new velocity + glm::vec3 rule1_component = glm::vec3(0.f, 0.f, 0.f); + glm::vec3 rule2_component = glm::vec3(0.f, 0.f, 0.f); + glm::vec3 rule3_component = glm::vec3(0.f, 0.f, 0.f); + if (neighbor_count_rule1 > 0) + { + rule1_component = (perceived_center / neighbor_count_rule1 - curr_boid_pos) * rule1Scale; + } + rule2_component = avoidance_velocity * rule2Scale; + if (neighbor_count_rule3 > 0) + { + rule3_component = (perceived_velocity / neighbor_count_rule3) * rule3Scale; + } + new_velocity += vel1[idx] + rule1_component + rule2_component + rule3_component; // - Clamp the speed change before putting the new speed in vel2 + float curr_speed = glm::length(new_velocity); + //if the total speed of vel is larger than maxSpeed, we normalize the vel and apply the maxSpeed we allow + if (curr_speed > maxSpeed) + { + new_velocity = glm::normalize(new_velocity) * maxSpeed; + } + + vel2[idx] = new_velocity; } + + __global__ void kernUpdateVelNeighborSearchCoherent( int N, int gridResolution, glm::vec3 gridMin, float inverseCellWidth, float cellWidth, @@ -424,15 +624,41 @@ void Boids::stepSimulationScatteredGrid(float dt) { // TODO-2.1 // Uniform Grid Neighbor search using Thrust sort. // In Parallel: + // - label each particle with its array index as well as its grid index. // Use 2x width grids. + + //set up blocks + int gridSizeParticle = (numObjects + blockSize - 1) / blockSize; //helped by Gangzheng Tong + dim3 blocksPerGridParticle(gridSizeParticle); + int gridSizeGridCell = (gridCellCount + blockSize - 1) / blockSize; + dim3 blocksPerGridGridCell(gridSizeGridCell); + //call kernel + kernComputeIndices << < blocksPerGridParticle, threadsPerBlock >> > (numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, dev_pos, dev_particleArrayIndices, dev_particleGridIndices); + // - Unstable key sort using Thrust. A stable sort isn't necessary, but you // are welcome to do a performance comparison. + thrust::sort_by_key(dev_thrust_particleGridIndices, dev_thrust_particleGridIndices + numObjects, dev_thrust_particleArrayIndices); //sort by grid index + // - Naively unroll the loop for finding the start and end indices of each // cell's data pointers in the array of boid indices + //first initialize the two indices list to be -1 + kernResetIntBuffer <<< blocksPerGridGridCell, threadsPerBlock >>> (gridCellCount, dev_gridCellStartIndices, -1); + kernResetIntBuffer <<< blocksPerGridGridCell, threadsPerBlock >>> (gridCellCount, dev_gridCellEndIndices, -1); + + //then call the kernel to compute those who contain boids + kernIdentifyCellStartEnd <<< blocksPerGridParticle, threadsPerBlock >>> (numObjects, dev_particleGridIndices, dev_gridCellStartIndices, dev_gridCellEndIndices); // - Perform velocity updates using neighbor search + kernUpdateVelNeighborSearchScattered << < blocksPerGridParticle, threadsPerBlock >> > ( + numObjects, gridSideCount, gridMinimum, + gridInverseCellWidth, gridCellWidth, + dev_gridCellStartIndices, dev_gridCellEndIndices, + dev_particleArrayIndices, + dev_pos, dev_vel1, dev_vel2); // - Update positions + kernUpdatePos << < blocksPerGridParticle, threadsPerBlock >> > (numObjects, dt, dev_pos, dev_vel2); // - Ping-pong buffers as needed + cudaMemcpy(dev_vel1, dev_vel2, sizeof(glm::vec3) * numObjects, cudaMemcpyDeviceToDevice); } void Boids::stepSimulationCoherentGrid(float dt) { @@ -459,68 +685,73 @@ void Boids::endSimulation() { cudaFree(dev_pos); // TODO-2.1 TODO-2.3 - Free any additional buffers here. + //2.1 free + cudaFree(dev_particleArrayIndices); + cudaFree(dev_particleGridIndices); + cudaFree(dev_gridCellStartIndices); + cudaFree(dev_gridCellEndIndices); } void Boids::unitTest() { - // LOOK-1.2 Feel free to write additional tests here. - - // test unstable sort - int *dev_intKeys; - int *dev_intValues; - int N = 10; - - std::unique_ptrintKeys{ new int[N] }; - std::unique_ptrintValues{ new int[N] }; - - intKeys[0] = 0; intValues[0] = 0; - intKeys[1] = 1; intValues[1] = 1; - intKeys[2] = 0; intValues[2] = 2; - intKeys[3] = 3; intValues[3] = 3; - intKeys[4] = 0; intValues[4] = 4; - intKeys[5] = 2; intValues[5] = 5; - intKeys[6] = 2; intValues[6] = 6; - intKeys[7] = 0; intValues[7] = 7; - intKeys[8] = 5; intValues[8] = 8; - intKeys[9] = 6; intValues[9] = 9; - - cudaMalloc((void**)&dev_intKeys, N * sizeof(int)); - checkCUDAErrorWithLine("cudaMalloc dev_intKeys failed!"); - - cudaMalloc((void**)&dev_intValues, N * sizeof(int)); - checkCUDAErrorWithLine("cudaMalloc dev_intValues failed!"); - - dim3 fullBlocksPerGrid((N + blockSize - 1) / blockSize); - - std::cout << "before unstable sort: " << std::endl; - for (int i = 0; i < N; i++) { - std::cout << " key: " << intKeys[i]; - std::cout << " value: " << intValues[i] << std::endl; - } + // LOOK-1.2 Feel free to write additional tests here. + + // test unstable sort + int *dev_intKeys; + int *dev_intValues; + int N = 10; + + std::unique_ptrintKeys{ new int[N] }; + std::unique_ptrintValues{ new int[N] }; + + intKeys[0] = 0; intValues[0] = 0; + intKeys[1] = 1; intValues[1] = 1; + intKeys[2] = 0; intValues[2] = 2; + intKeys[3] = 3; intValues[3] = 3; + intKeys[4] = 0; intValues[4] = 4; + intKeys[5] = 2; intValues[5] = 5; + intKeys[6] = 2; intValues[6] = 6; + intKeys[7] = 0; intValues[7] = 7; + intKeys[8] = 5; intValues[8] = 8; + intKeys[9] = 6; intValues[9] = 9; + + cudaMalloc((void**)&dev_intKeys, N * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc dev_intKeys failed!"); + + cudaMalloc((void**)&dev_intValues, N * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc dev_intValues failed!"); + + dim3 fullBlocksPerGrid((N + blockSize - 1) / blockSize); + + std::cout << "before unstable sort: " << std::endl; + for (int i = 0; i < N; i++) { + std::cout << " key: " << intKeys[i]; + std::cout << " value: " << intValues[i] << std::endl; + } - // How to copy data to the GPU - cudaMemcpy(dev_intKeys, intKeys.get(), sizeof(int) * N, cudaMemcpyHostToDevice); - cudaMemcpy(dev_intValues, intValues.get(), sizeof(int) * N, cudaMemcpyHostToDevice); - - // Wrap device vectors in thrust iterators for use with thrust. - thrust::device_ptr dev_thrust_keys(dev_intKeys); - thrust::device_ptr dev_thrust_values(dev_intValues); - // LOOK-2.1 Example for using thrust::sort_by_key - thrust::sort_by_key(dev_thrust_keys, dev_thrust_keys + N, dev_thrust_values); - - // How to copy data back to the CPU side from the GPU - cudaMemcpy(intKeys.get(), dev_intKeys, sizeof(int) * N, cudaMemcpyDeviceToHost); - cudaMemcpy(intValues.get(), dev_intValues, sizeof(int) * N, cudaMemcpyDeviceToHost); - checkCUDAErrorWithLine("memcpy back failed!"); - - std::cout << "after unstable sort: " << std::endl; - for (int i = 0; i < N; i++) { - std::cout << " key: " << intKeys[i]; - std::cout << " value: " << intValues[i] << std::endl; - } + // How to copy data to the GPU + cudaMemcpy(dev_intKeys, intKeys.get(), sizeof(int) * N, cudaMemcpyHostToDevice); + cudaMemcpy(dev_intValues, intValues.get(), sizeof(int) * N, cudaMemcpyHostToDevice); + + // Wrap device vectors in thrust iterators for use with thrust. + thrust::device_ptr dev_thrust_keys(dev_intKeys); + thrust::device_ptr dev_thrust_values(dev_intValues); + // LOOK-2.1 Example for using thrust::sort_by_key + thrust::sort_by_key(dev_thrust_keys, dev_thrust_keys + N, dev_thrust_values); + + // How to copy data back to the CPU side from the GPU + cudaMemcpy(intKeys.get(), dev_intKeys, sizeof(int) * N, cudaMemcpyDeviceToHost); + cudaMemcpy(intValues.get(), dev_intValues, sizeof(int) * N, cudaMemcpyDeviceToHost); + checkCUDAErrorWithLine("memcpy back failed!"); + + std::cout << "after unstable sort: " << std::endl; + for (int i = 0; i < N; i++) { + std::cout << " key: " << intKeys[i]; + std::cout << " value: " << intValues[i] << std::endl; + } - // cleanup - cudaFree(dev_intKeys); - cudaFree(dev_intValues); - checkCUDAErrorWithLine("cudaFree failed!"); - return; -} + // cleanup + cudaFree(dev_intKeys); + cudaFree(dev_intValues); + checkCUDAErrorWithLine("cudaFree failed!"); + return; +} \ No newline at end of file diff --git a/src/main.cpp b/src/main.cpp index a7b8586..e3b6f63 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -14,11 +14,11 @@ // LOOK-2.1 LOOK-2.3 - toggles for UNIFORM_GRID and COHERENT_GRID #define VISUALIZE 1 -#define UNIFORM_GRID 0 +#define UNIFORM_GRID 1 #define COHERENT_GRID 0 // LOOK-1.2 - change this to adjust particle count in the simulation -const int N_FOR_VIS = 10000; +const int N_FOR_VIS = 50000; const float DT = 0.2f; /** From 09aa3656749ea5f2307599b565aaa0246cb65f5d Mon Sep 17 00:00:00 2001 From: Mark Date: Fri, 6 Sep 2019 18:32:48 -0400 Subject: [PATCH 4/8] finish part 2.3, but the particle movement doesn't seems right --- src/kernel.cu | 247 ++++++++++++++++++++++++++++++++++++++++++++------ src/main.cpp | 6 +- 2 files changed, 224 insertions(+), 29 deletions(-) diff --git a/src/kernel.cu b/src/kernel.cu index ed9a521..9c7a6b8 100644 --- a/src/kernel.cu +++ b/src/kernel.cu @@ -87,6 +87,9 @@ int *dev_gridCellEndIndices; // to this cell? // TODO-2.3 - consider what additional buffers you might need to reshuffle // the position and velocity data to be coherent within cells. +glm::vec3 *dev_pos_buffer; +glm::vec3 *dev_vel1_buffer; +glm::vec3 *dev_vel2_buffer; // LOOK-2.1 - Grid parameters based on simulation parameters. // These are automatically computed for you in Boids::initSimulation @@ -188,6 +191,16 @@ void Boids::initSimulation(int N) { dev_thrust_particleArrayIndices = thrust::device_ptr(dev_particleArrayIndices); dev_thrust_particleGridIndices = thrust::device_ptr(dev_particleGridIndices); + //2.3 additional buffer to help reshuffle + cudaMalloc((void**)&dev_pos_buffer, N * sizeof(glm::vec3)); + checkCUDAErrorWithLine("cudaMalloc dev_pos_buffer failed!"); + + cudaMalloc((void**)&dev_vel1_buffer, N * sizeof(glm::vec3)); + checkCUDAErrorWithLine("cudaMalloc dev_vel1_buffer failed!"); + + cudaMalloc((void**)&dev_vel2_buffer, N * sizeof(glm::vec3)); + checkCUDAErrorWithLine("cudaMalloc dev_vel2_buffer failed!"); + cudaDeviceSynchronize(); } @@ -362,7 +375,7 @@ __global__ void kernUpdatePos(int N, float dt, glm::vec3 *pos, glm::vec3 *vel) { // LOOK-2.1 Consider this method of computing a 1D index from a 3D grid index. // LOOK-2.3 Looking at this method, what would be the most memory efficient // order for iterating over neighboring grid cells? -// for(x) +// for(x) tick // for(y) // for(z)? Or some other order? __device__ int gridIndex3Dto1D(int x, int y, int z, int gridResolution) { @@ -447,13 +460,6 @@ __global__ void kernIdentifyCellStartEnd(int N, int *particleGridIndices, } - //int next_grid_index = particleGridIndices[idx + 1]; - //if (target_grid_index != next_grid_index) - //{ - // gridCellStartIndices[next_grid_index] = idx + 1; - // gridCellEndIndices[target_grid_index] = idx; - //} - } //very similar to photon mapping @@ -501,6 +507,7 @@ __global__ void kernUpdateVelNeighborSearchScattered( glm::vec3 new_velocity(0.f, 0.f, 0.f); float neighbor_count_rule1 = 0; float neighbor_count_rule3 = 0; + //is this a grid-looping optimization? for (int x_cord_idx = min_x; x_cord_idx <= max_x; ++x_cord_idx) { for (int y_cord_idx = min_y; y_cord_idx <= max_y; ++y_cord_idx) @@ -570,7 +577,7 @@ __global__ void kernUpdateVelNeighborSearchScattered( { rule3_component = (perceived_velocity / neighbor_count_rule3) * rule3Scale; } - new_velocity += vel1[idx] + rule1_component + rule2_component + rule3_component; + new_velocity = vel1[idx] + rule1_component + rule2_component + rule3_component; // - Clamp the speed change before putting the new speed in vel2 float curr_speed = glm::length(new_velocity); //if the total speed of vel is larger than maxSpeed, we normalize the vel and apply the maxSpeed we allow @@ -594,13 +601,154 @@ __global__ void kernUpdateVelNeighborSearchCoherent( // This should expect gridCellStartIndices and gridCellEndIndices to refer // directly to pos and vel1. // - Identify the grid cell that this particle is in - // - Identify which cells may contain neighbors. This isn't always 8. - // - For each cell, read the start/end indices in the boid pointer array. // DIFFERENCE: For best results, consider what order the cells should be // checked in to maximize the memory benefits of reordering the boids data. // - Access each boid in the cell and compute velocity change from // the boids rules, if this boid is within the neighborhood distance. // - Clamp the speed change before putting the new speed in vel2 + + + + int idx = threadIdx.x + (blockIdx.x * blockDim.x); + if (idx >= N) { + return; + } + glm::vec3 curr_boid_pos = pos[idx]; + int idx_x = (curr_boid_pos.x - gridMin.x) * inverseCellWidth; + int idx_y = (curr_boid_pos.y - gridMin.y) * inverseCellWidth; + int idx_z = (curr_boid_pos.z - gridMin.z) * inverseCellWidth; + + //combine to get the 1D index + int gridIndex = gridIndex3Dto1D(idx_x, idx_y, idx_z, gridResolution); + // - Identify which cells may contain neighbors. This isn't always 8. + //by calculating those cells that interact with the sphere with the neighbor_radius(std::max(std::max(rule1Distance, rule2Distance), rule3Distance)) + //float neighbor_radius = cellWidth / CELL_WIDTH_FACTOR; + const float neighbor_radius = glm::max(glm::max(rule1Distance, rule2Distance), rule3Distance); + int max_x = (curr_boid_pos.x + neighbor_radius - gridMin.x) * inverseCellWidth; + max_x = max_x > gridResolution - 1 ? gridResolution - 1 : max_x; + int min_x = (curr_boid_pos.x - neighbor_radius - gridMin.x) * inverseCellWidth; + min_x = min_x < 0 ? 0 : min_x; + int max_y = (curr_boid_pos.y + neighbor_radius - gridMin.y) * inverseCellWidth; + max_y = max_y > gridResolution - 1 ? gridResolution - 1 : max_y; + int min_y = (curr_boid_pos.y + neighbor_radius - gridMin.y) * inverseCellWidth; + min_y = min_y < 0 ? 0 : min_y; + int max_z = (curr_boid_pos.z + neighbor_radius - gridMin.z) * inverseCellWidth; + max_z = max_z > gridResolution - 1 ? gridResolution - 1 : max_z; + int min_z = (curr_boid_pos.z + neighbor_radius - gridMin.z) * inverseCellWidth; + min_z = min_z < 0 ? 0 : min_z; + + // - For each cell, read the start/end indices in the boid pointer array. + glm::vec3 perceived_center(0.f, 0.f, 0.f); + glm::vec3 avoidance_velocity(0.f, 0.f, 0.f); + glm::vec3 perceived_velocity(0.f, 0.f, 0.f); + glm::vec3 new_velocity(0.f, 0.f, 0.f); + float neighbor_count_rule1 = 0; + float neighbor_count_rule3 = 0; + //is this a grid-looping optimization? + // DIFFERENCE: For best results, consider what order the cells should be + // checked in to maximize the memory benefits of reordering the boids data. + for (int z_cord_idx = min_z; z_cord_idx <= max_z; ++z_cord_idx) + { + for (int y_cord_idx = min_y; y_cord_idx <= max_y; ++y_cord_idx) + { + for (int x_cord_idx = min_x; x_cord_idx <= max_x; ++x_cord_idx) + { + int curr_gridIndex = gridIndex3Dto1D(x_cord_idx, y_cord_idx, z_cord_idx, gridResolution); + //read the start/end indices + int start = gridCellStartIndices[curr_gridIndex]; + int end = gridCellEndIndices[curr_gridIndex]; + if (start == -1 || end == -1) { + continue; //no boid in this cell + } + else + { + // - Access each boid in the cell and compute velocity change from + // the boids rules, if this boid is within the neighborhood distance. + for (int boid_array_idx = start; boid_array_idx <= end; ++boid_array_idx) + { + //the boid_array_idx is only the index in particleArray, need to load it out -- buggy for only few particles moving + int boid_idx = boid_array_idx; + + if (boid_idx == idx) continue; + + //pre load temp boid pos + glm::vec3 idx_boid_pos = pos[boid_idx]; + float dist = glm::distance(idx_boid_pos, curr_boid_pos); + // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves + if (dist < rule1Distance) + { + perceived_center += idx_boid_pos; + neighbor_count_rule1++; + } + + // Rule 2: boids try to stay a distance d away from each other + if (dist < rule2Distance) + { + avoidance_velocity -= (idx_boid_pos - curr_boid_pos); + } + + // Rule 3: boids try to match the speed of surrounding boids + if (dist < rule3Distance) + { + perceived_velocity += vel1[idx]; + neighbor_count_rule3++; + } + + } + + } + } + + } + + } + + //compute the new velocity + glm::vec3 rule1_component = glm::vec3(0.f, 0.f, 0.f); + glm::vec3 rule2_component = glm::vec3(0.f, 0.f, 0.f); + glm::vec3 rule3_component = glm::vec3(0.f, 0.f, 0.f); + if (neighbor_count_rule1 > 0) + { + rule1_component = (perceived_center / neighbor_count_rule1 - curr_boid_pos) * rule1Scale; + } + rule2_component = avoidance_velocity * rule2Scale; + if (neighbor_count_rule3 > 0) + { + rule3_component = (perceived_velocity / neighbor_count_rule3) * rule3Scale; + } + new_velocity = vel1[idx] + rule1_component + rule2_component + rule3_component; + // - Clamp the speed change before putting the new speed in vel2 + float curr_speed = glm::length(new_velocity); + //if the total speed of vel is larger than maxSpeed, we normalize the vel and apply the maxSpeed we allow + if (curr_speed > maxSpeed) + { + new_velocity = glm::normalize(new_velocity) * maxSpeed; + } + + vel2[idx] = new_velocity; +} + + +__global__ void kernReshuffleArray(int N, int *particle_array_indices, glm::vec3 *source, glm::vec3 *destination) +{ + int idx = threadIdx.x + (blockIdx.x * blockDim.x); + if (idx >= N) { + return; + } + + int particle_idx = particle_array_indices[idx]; + destination[idx] = source[particle_idx]; +} + +__global__ void kernReverseArray(int N, int *particle_array_indices, glm::vec3 *source, glm::vec3 *destination) +{ + int idx = threadIdx.x + (blockIdx.x * blockDim.x); + if (idx >= N) { + return; + } + + int particle_idx = particle_array_indices[idx]; + destination[particle_idx] = source[idx]; } /** @@ -617,7 +765,9 @@ void Boids::stepSimulationNaive(float dt) { kernUpdatePos << < blocksPerGrid, threadsPerBlock >> > (numObjects, dt, dev_pos, dev_vel2); // TODO-1.2 ping-pong the velocity buffers -- swap content - cudaMemcpy(dev_vel1, dev_vel2, sizeof(glm::vec3) * numObjects, cudaMemcpyDeviceToDevice); + glm::vec3 *temp_vel = dev_vel1; + dev_vel1 = dev_vel2; + dev_vel2 = temp_vel; } void Boids::stepSimulationScatteredGrid(float dt) { @@ -658,25 +808,66 @@ void Boids::stepSimulationScatteredGrid(float dt) { // - Update positions kernUpdatePos << < blocksPerGridParticle, threadsPerBlock >> > (numObjects, dt, dev_pos, dev_vel2); // - Ping-pong buffers as needed - cudaMemcpy(dev_vel1, dev_vel2, sizeof(glm::vec3) * numObjects, cudaMemcpyDeviceToDevice); + glm::vec3 *temp_vel = dev_vel1; + dev_vel1 = dev_vel2; + dev_vel2 = temp_vel; } void Boids::stepSimulationCoherentGrid(float dt) { - // TODO-2.3 - start by copying Boids::stepSimulationNaiveGrid - // Uniform Grid Neighbor search using Thrust sort on cell-coherent data. - // In Parallel: - // - Label each particle with its array index as well as its grid index. - // Use 2x width grids - // - Unstable key sort using Thrust. A stable sort isn't necessary, but you - // are welcome to do a performance comparison. + // TODO-2.3 - start by copying Boids::stepSimulationNaiveGrid + // Uniform Grid Neighbor search using Thrust sort on cell-coherent data. + // In Parallel: + // - Label each particle with its array index as well as its grid index. + // Use 2x width grids + int gridSizeParticle = (numObjects + blockSize - 1) / blockSize; //helped by Gangzheng Tong + dim3 blocksPerGridParticle(gridSizeParticle); + int gridSizeGridCell = (gridCellCount + blockSize - 1) / blockSize; + dim3 blocksPerGridGridCell(gridSizeGridCell); + + kernComputeIndices << < blocksPerGridParticle, threadsPerBlock >> > (numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, dev_pos, dev_particleArrayIndices, dev_particleGridIndices); + // - Unstable key sort using Thrust. A stable sort isn't necessary, but you + // are welcome to do a performance comparison. + thrust::sort_by_key(dev_thrust_particleGridIndices, dev_thrust_particleGridIndices + numObjects, dev_thrust_particleArrayIndices); //sort by grid index + + // - Naively unroll the loop for finding the start and end indices of each // cell's data pointers in the array of boid indices - // - BIG DIFFERENCE: use the rearranged array index buffer to reshuffle all - // the particle data in the simulation array. - // CONSIDER WHAT ADDITIONAL BUFFERS YOU NEED - // - Perform velocity updates using neighbor search - // - Update positions - // - Ping-pong buffers as needed. THIS MAY BE DIFFERENT FROM BEFORE. + //first initialize the two indices list to be -1 + kernResetIntBuffer << < blocksPerGridGridCell, threadsPerBlock >> > (gridCellCount, dev_gridCellStartIndices, -1); + kernResetIntBuffer << < blocksPerGridGridCell, threadsPerBlock >> > (gridCellCount, dev_gridCellEndIndices, -1); + + kernIdentifyCellStartEnd << < blocksPerGridParticle, threadsPerBlock >> > (numObjects, dev_particleGridIndices, dev_gridCellStartIndices, dev_gridCellEndIndices); + // - BIG DIFFERENCE: use the rearranged array index buffer to reshuffle all + // the particle data in the simulation array. + // CONSIDER WHAT ADDITIONAL BUFFERS YOU NEED + + kernReshuffleArray << < blocksPerGridParticle, threadsPerBlock >> > (numObjects, dev_particleArrayIndices, dev_pos, dev_pos_buffer); + kernReshuffleArray << < blocksPerGridParticle, threadsPerBlock >> > (numObjects, dev_particleArrayIndices, dev_vel1, dev_vel1_buffer); + //reassign the pos and vel buffer back to pos and vel1 + cudaMemcpy(dev_pos, dev_pos_buffer, sizeof(glm::vec3) * numObjects, cudaMemcpyDeviceToDevice); + cudaMemcpy(dev_vel1, dev_vel1_buffer, sizeof(glm::vec3) * numObjects, cudaMemcpyDeviceToDevice); + + // - Perform velocity updates using neighbor search + kernUpdateVelNeighborSearchCoherent << < blocksPerGridParticle, threadsPerBlock >> > ( + numObjects, gridSideCount, gridMinimum, + gridInverseCellWidth, gridCellWidth, + dev_gridCellStartIndices, dev_gridCellEndIndices, + dev_pos, dev_vel1, dev_vel2); + + //dev_pos has been reshuffled -- dev_vel2 also need to bu reshuffled + //kernReshuffleArray << < blocksPerGridParticle, threadsPerBlock >> > (numObjects, dev_particleArrayIndices, dev_vel2, dev_vel2_buffer); + //cudaMemcpy(dev_vel2, dev_vel2_buffer, sizeof(glm::vec3) * numObjects, cudaMemcpyDeviceToDevice); + // - Update positions + kernUpdatePos << < blocksPerGridParticle, threadsPerBlock >> > (numObjects, dt, dev_pos, dev_vel2); + // - Ping-pong buffers as needed. THIS MAY BE DIFFERENT FROM BEFORE. -- shuffle back and store back + //kernReverseArray <<< blocksPerGridParticle, threadsPerBlock >> > (numObjects, dev_particleArrayIndices, dev_vel2, dev_vel2_buffer); + //kernReverseArray << < blocksPerGridParticle, threadsPerBlock >> > (numObjects, dev_particleArrayIndices, dev_pos, dev_pos_buffer); + //cudaMemcpy(dev_vel1, dev_vel2_buffer, sizeof(glm::vec3) * numObjects, cudaMemcpyDeviceToDevice); + //cudaMemcpy(dev_pos, dev_pos_buffer, sizeof(glm::vec3) * numObjects, cudaMemcpyDeviceToDevice); + + glm::vec3 *temp_vel = dev_vel1; + dev_vel1 = dev_vel2; + dev_vel2 = temp_vel; } void Boids::endSimulation() { @@ -690,6 +881,10 @@ void Boids::endSimulation() { cudaFree(dev_particleGridIndices); cudaFree(dev_gridCellStartIndices); cudaFree(dev_gridCellEndIndices); + + cudaFree(dev_vel1_buffer); + cudaFree(dev_vel2_buffer); + cudaFree(dev_pos_buffer); } void Boids::unitTest() { diff --git a/src/main.cpp b/src/main.cpp index e3b6f63..6641a62 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -13,12 +13,12 @@ // ================ // LOOK-2.1 LOOK-2.3 - toggles for UNIFORM_GRID and COHERENT_GRID -#define VISUALIZE 1 +#define VISUALIZE 0 #define UNIFORM_GRID 1 -#define COHERENT_GRID 0 +#define COHERENT_GRID 1 // LOOK-1.2 - change this to adjust particle count in the simulation -const int N_FOR_VIS = 50000; +const int N_FOR_VIS = 5000; const float DT = 0.2f; /** From 6239c49fa4e2b2be7b221dbf954a75e460d28f54 Mon Sep 17 00:00:00 2001 From: Mark Date: Sat, 7 Sep 2019 20:49:48 -0400 Subject: [PATCH 5/8] finish part 2 --- src/kernel.cu | 79 ++++++++++++++++++++++++++++++++------------------- 1 file changed, 49 insertions(+), 30 deletions(-) diff --git a/src/kernel.cu b/src/kernel.cu index 9c7a6b8..7a23618 100644 --- a/src/kernel.cu +++ b/src/kernel.cu @@ -394,6 +394,13 @@ __global__ void kernComputeIndices(int N, int gridResolution, } //pre store info glm::vec3 curr_pos = pos[idx]; + + //test + //indices[idx] = idx; + //const glm::vec3 gridIndex = glm::floor((curr_pos - gridMin) * inverseCellWidth); + //gridIndices[idx] = gridIndex3Dto1D(gridIndex.x, gridIndex.y, gridIndex.z, gridResolution); + + //original int idx_x = (curr_pos.x - gridMin.x) * inverseCellWidth; int idx_y = (curr_pos.y - gridMin.y) * inverseCellWidth; int idx_z = (curr_pos.z - gridMin.z) * inverseCellWidth; @@ -445,7 +452,16 @@ __global__ void kernIdentifyCellStartEnd(int N, int *particleGridIndices, return; } + //test + //const int nextGrid = particleGridIndices[idx + 1]; + + //if (target_grid_index != nextGrid) + //{ + // gridCellEndIndices[target_grid_index] = idx; + // gridCellStartIndices[nextGrid] = idx + 1; + //} + //original //check one before and one after if (target_grid_index != 0 && target_grid_index != particleGridIndices[idx - 1]) { @@ -469,7 +485,7 @@ __global__ void kernUpdateVelNeighborSearchScattered( int *gridCellStartIndices, int *gridCellEndIndices, int *particleArrayIndices, glm::vec3 *pos, glm::vec3 *vel1, glm::vec3 *vel2) { - // TODO-2.1 - Update a boid's velocity using the uniform grid to reduce + // TODO-2.1 - Update a boid's velocity using the uniform grid to reduce // the number of boids that need to be checked. // - Identify the grid cell that this particle is in int idx = threadIdx.x + (blockIdx.x * blockDim.x); @@ -477,36 +493,37 @@ __global__ void kernUpdateVelNeighborSearchScattered( return; } glm::vec3 curr_boid_pos = pos[idx]; + int idx_x = (curr_boid_pos.x - gridMin.x) * inverseCellWidth; int idx_y = (curr_boid_pos.y - gridMin.y) * inverseCellWidth; int idx_z = (curr_boid_pos.z - gridMin.z) * inverseCellWidth; //combine to get the 1D index int gridIndex = gridIndex3Dto1D(idx_x, idx_y, idx_z, gridResolution); - // - Identify which cells may contain neighbors. This isn't always 8. - //by calculating those cells that interact with the sphere with the neighbor_radius(std::max(std::max(rule1Distance, rule2Distance), rule3Distance)) -not allowed - //float neighbor_radius = cellWidth / CELL_WIDTH_FACTOR; + // - Identify which cells may contain neighbors. This isn't always 8. + + //by calculating those cells that interact with the sphere with the neighbor_radius(std::max(std::max(rule1Distance, rule2Distance), rule3Distance)) -not allowed unless it is const float const float neighbor_radius = glm::max(glm::max(rule1Distance, rule2Distance), rule3Distance); - int max_x = (curr_boid_pos.x + neighbor_radius - gridMin.x) * inverseCellWidth; + int max_x = glm::floor((curr_boid_pos.x + neighbor_radius - gridMin.x) * inverseCellWidth); max_x = max_x > gridResolution - 1 ? gridResolution - 1 : max_x; - int min_x = (curr_boid_pos.x - neighbor_radius - gridMin.x) * inverseCellWidth; + int min_x = glm::floor((curr_boid_pos.x - neighbor_radius - gridMin.x) * inverseCellWidth); min_x = min_x < 0 ? 0 : min_x; - int max_y = (curr_boid_pos.y + neighbor_radius - gridMin.y) * inverseCellWidth; + int max_y = glm::floor((curr_boid_pos.y + neighbor_radius - gridMin.y) * inverseCellWidth); max_y = max_y > gridResolution - 1 ? gridResolution - 1 : max_y; - int min_y = (curr_boid_pos.y + neighbor_radius - gridMin.y) * inverseCellWidth; + int min_y = glm::floor((curr_boid_pos.y + neighbor_radius - gridMin.y) * inverseCellWidth); min_y = min_y < 0 ? 0 : min_y; - int max_z = (curr_boid_pos.z + neighbor_radius - gridMin.z) * inverseCellWidth; + int max_z = glm::floor((curr_boid_pos.z + neighbor_radius - gridMin.z) * inverseCellWidth); max_z = max_z > gridResolution - 1 ? gridResolution - 1 : max_z; - int min_z = (curr_boid_pos.z + neighbor_radius - gridMin.z) * inverseCellWidth; + int min_z = glm::floor((curr_boid_pos.z + neighbor_radius - gridMin.z) * inverseCellWidth); min_z = min_z < 0 ? 0 : min_z; - // - For each cell, read the start/end indices in the boid pointer array. + // - For each cell, read the start/end indices in the boid pointer array. glm::vec3 perceived_center(0.f, 0.f, 0.f); glm::vec3 avoidance_velocity(0.f, 0.f, 0.f); glm::vec3 perceived_velocity(0.f, 0.f, 0.f); glm::vec3 new_velocity(0.f, 0.f, 0.f); - float neighbor_count_rule1 = 0; - float neighbor_count_rule3 = 0; + int neighbor_count_rule1 = 0; + int neighbor_count_rule3 = 0; //is this a grid-looping optimization? for (int x_cord_idx = min_x; x_cord_idx <= max_x; ++x_cord_idx) { @@ -534,34 +551,36 @@ __global__ void kernUpdateVelNeighborSearchScattered( //pre load temp boid pos glm::vec3 idx_boid_pos = pos[boid_idx]; - float dist = glm::distance(idx_boid_pos, curr_boid_pos); + float dist = glm::distance(curr_boid_pos, idx_boid_pos); // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves - if (dist < rule1Distance) + if (dist <= rule1Distance) { perceived_center += idx_boid_pos; neighbor_count_rule1++; } // Rule 2: boids try to stay a distance d away from each other - if (dist < rule2Distance) + if (dist <= rule2Distance) { avoidance_velocity -= (idx_boid_pos - curr_boid_pos); } // Rule 3: boids try to match the speed of surrounding boids - if (dist < rule3Distance) + if (dist <= rule3Distance) { - perceived_velocity += vel1[idx]; + //use the boid_idx to access its speed, not the origianl speed .... + //perceived_velocity += vel1[idx]; + perceived_velocity += vel1[boid_idx]; neighbor_count_rule3++; } - + } - + } } - + } - + } //compute the new velocity @@ -570,15 +589,15 @@ __global__ void kernUpdateVelNeighborSearchScattered( glm::vec3 rule3_component = glm::vec3(0.f, 0.f, 0.f); if (neighbor_count_rule1 > 0) { - rule1_component = (perceived_center / neighbor_count_rule1 - curr_boid_pos) * rule1Scale; + rule1_component = (perceived_center / (float)neighbor_count_rule1 - curr_boid_pos) * rule1Scale; } rule2_component = avoidance_velocity * rule2Scale; if (neighbor_count_rule3 > 0) { - rule3_component = (perceived_velocity / neighbor_count_rule3) * rule3Scale; + rule3_component = (perceived_velocity / (float)neighbor_count_rule3) * rule3Scale; } new_velocity = vel1[idx] + rule1_component + rule2_component + rule3_component; - // - Clamp the speed change before putting the new speed in vel2 + // - Clamp the speed change before putting the new speed in vel2 float curr_speed = glm::length(new_velocity); //if the total speed of vel is larger than maxSpeed, we normalize the vel and apply the maxSpeed we allow if (curr_speed > maxSpeed) @@ -642,8 +661,8 @@ __global__ void kernUpdateVelNeighborSearchCoherent( glm::vec3 avoidance_velocity(0.f, 0.f, 0.f); glm::vec3 perceived_velocity(0.f, 0.f, 0.f); glm::vec3 new_velocity(0.f, 0.f, 0.f); - float neighbor_count_rule1 = 0; - float neighbor_count_rule3 = 0; + int neighbor_count_rule1 = 0; + int neighbor_count_rule3 = 0; //is this a grid-looping optimization? // DIFFERENCE: For best results, consider what order the cells should be // checked in to maximize the memory benefits of reordering the boids data. @@ -690,7 +709,7 @@ __global__ void kernUpdateVelNeighborSearchCoherent( // Rule 3: boids try to match the speed of surrounding boids if (dist < rule3Distance) { - perceived_velocity += vel1[idx]; + perceived_velocity += vel1[boid_idx]; neighbor_count_rule3++; } @@ -709,12 +728,12 @@ __global__ void kernUpdateVelNeighborSearchCoherent( glm::vec3 rule3_component = glm::vec3(0.f, 0.f, 0.f); if (neighbor_count_rule1 > 0) { - rule1_component = (perceived_center / neighbor_count_rule1 - curr_boid_pos) * rule1Scale; + rule1_component = (perceived_center / (float)neighbor_count_rule1 - curr_boid_pos) * rule1Scale; } rule2_component = avoidance_velocity * rule2Scale; if (neighbor_count_rule3 > 0) { - rule3_component = (perceived_velocity / neighbor_count_rule3) * rule3Scale; + rule3_component = (perceived_velocity / (float)neighbor_count_rule3) * rule3Scale; } new_velocity = vel1[idx] + rule1_component + rule2_component + rule3_component; // - Clamp the speed change before putting the new speed in vel2 From c444b16cffadc441ef05a1ee607988fe93db2606 Mon Sep 17 00:00:00 2001 From: Mark Date: Sat, 7 Sep 2019 21:28:11 -0400 Subject: [PATCH 6/8] really finish part 2 by fixing calculating on individual coordinate index to all 3 coordinate indices at the same time --- src/kernel.cu | 93 ++++++++++++++++++++++++++------------------------- 1 file changed, 48 insertions(+), 45 deletions(-) diff --git a/src/kernel.cu b/src/kernel.cu index 7a23618..da33fa0 100644 --- a/src/kernel.cu +++ b/src/kernel.cu @@ -267,8 +267,8 @@ __device__ glm::vec3 computeVelocityChange(int N, int iSelf, const glm::vec3 *po glm::vec3 avoidance_velocity(0.f, 0.f, 0.f); glm::vec3 perceived_velocity(0.f, 0.f, 0.f); glm::vec3 return_vel(0.f, 0.f, 0.f); - float neighbor_count_rule1 = 0; - float neighbor_count_rule3 = 0; + int neighbor_count_rule1 = 0; + int neighbor_count_rule3 = 0; //pre load all needed data glm::vec3 curr_boid_pos = pos[iSelf]; @@ -281,20 +281,20 @@ __device__ glm::vec3 computeVelocityChange(int N, int iSelf, const glm::vec3 *po glm::vec3 idx_boid_pos = pos[idx]; float dist = glm::distance(idx_boid_pos, curr_boid_pos); // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves - if (dist < rule1Distance) + if (dist <= rule1Distance) { perceived_center += idx_boid_pos; neighbor_count_rule1++; } // Rule 2: boids try to stay a distance d away from each other - if (dist < rule2Distance) + if (dist <= rule2Distance) { avoidance_velocity -= (idx_boid_pos - curr_boid_pos); } // Rule 3: boids try to match the speed of surrounding boids - if (dist < rule3Distance) + if (dist <= rule3Distance) { perceived_velocity += vel[idx]; neighbor_count_rule3++; @@ -307,12 +307,12 @@ __device__ glm::vec3 computeVelocityChange(int N, int iSelf, const glm::vec3 *po glm::vec3 rule3_component = glm::vec3(0.f, 0.f, 0.f); if (neighbor_count_rule1 > 0) { - rule1_component = (perceived_center / neighbor_count_rule1 - curr_boid_pos) * rule1Scale; + rule1_component = (perceived_center / (float)neighbor_count_rule1 - curr_boid_pos) * rule1Scale; } rule2_component = avoidance_velocity * rule2Scale; if (neighbor_count_rule3 > 0) { - rule3_component = (perceived_velocity / neighbor_count_rule3) * rule3Scale; + rule3_component = (perceived_velocity / (float)neighbor_count_rule3) * rule3Scale; } //helped by Hanna ReadMe Rule part sum all rules' and current velocity @@ -395,12 +395,6 @@ __global__ void kernComputeIndices(int N, int gridResolution, //pre store info glm::vec3 curr_pos = pos[idx]; - //test - //indices[idx] = idx; - //const glm::vec3 gridIndex = glm::floor((curr_pos - gridMin) * inverseCellWidth); - //gridIndices[idx] = gridIndex3Dto1D(gridIndex.x, gridIndex.y, gridIndex.z, gridResolution); - - //original int idx_x = (curr_pos.x - gridMin.x) * inverseCellWidth; int idx_y = (curr_pos.y - gridMin.y) * inverseCellWidth; int idx_z = (curr_pos.z - gridMin.z) * inverseCellWidth; @@ -486,8 +480,8 @@ __global__ void kernUpdateVelNeighborSearchScattered( int *particleArrayIndices, glm::vec3 *pos, glm::vec3 *vel1, glm::vec3 *vel2) { // TODO-2.1 - Update a boid's velocity using the uniform grid to reduce - // the number of boids that need to be checked. - // - Identify the grid cell that this particle is in +// the number of boids that need to be checked. +// - Identify the grid cell that this particle is in int idx = threadIdx.x + (blockIdx.x * blockDim.x); if (idx >= N) { return; @@ -504,19 +498,31 @@ __global__ void kernUpdateVelNeighborSearchScattered( //by calculating those cells that interact with the sphere with the neighbor_radius(std::max(std::max(rule1Distance, rule2Distance), rule3Distance)) -not allowed unless it is const float const float neighbor_radius = glm::max(glm::max(rule1Distance, rule2Distance), rule3Distance); - int max_x = glm::floor((curr_boid_pos.x + neighbor_radius - gridMin.x) * inverseCellWidth); - max_x = max_x > gridResolution - 1 ? gridResolution - 1 : max_x; - int min_x = glm::floor((curr_boid_pos.x - neighbor_radius - gridMin.x) * inverseCellWidth); - min_x = min_x < 0 ? 0 : min_x; - int max_y = glm::floor((curr_boid_pos.y + neighbor_radius - gridMin.y) * inverseCellWidth); - max_y = max_y > gridResolution - 1 ? gridResolution - 1 : max_y; - int min_y = glm::floor((curr_boid_pos.y + neighbor_radius - gridMin.y) * inverseCellWidth); - min_y = min_y < 0 ? 0 : min_y; - int max_z = glm::floor((curr_boid_pos.z + neighbor_radius - gridMin.z) * inverseCellWidth); - max_z = max_z > gridResolution - 1 ? gridResolution - 1 : max_z; - int min_z = glm::floor((curr_boid_pos.z + neighbor_radius - gridMin.z) * inverseCellWidth); - min_z = min_z < 0 ? 0 : min_z; - + //incorrect here, but why? + //int max_x = glm::floor((curr_boid_pos.x + neighbor_radius - gridMin.x) * inverseCellWidth); + //max_x = max_x > gridResolution - 1 ? gridResolution - 1 : max_x; + //int min_x = glm::floor((curr_boid_pos.x - neighbor_radius - gridMin.x) * inverseCellWidth); + //min_x = min_x < 0 ? 0 : min_x; + //int max_y = glm::floor((curr_boid_pos.y + neighbor_radius - gridMin.y) * inverseCellWidth); + //max_y = max_y > gridResolution - 1 ? gridResolution - 1 : max_y; + //int min_y = glm::floor((curr_boid_pos.y + neighbor_radius - gridMin.y) * inverseCellWidth); + //min_y = min_y < 0 ? 0 : min_y; + //int max_z = glm::floor((curr_boid_pos.z + neighbor_radius - gridMin.z) * inverseCellWidth); + //max_z = max_z > gridResolution - 1 ? gridResolution - 1 : max_z; + //int min_z = glm::floor((curr_boid_pos.z + neighbor_radius - gridMin.z) * inverseCellWidth); + //min_z = min_z < 0 ? 0 : min_z; + + //find the minimum index and max index of the cell we might need to use + glm::vec3 grid_index_min = glm::floor((curr_boid_pos - gridMin - glm::vec3(neighbor_radius)) * inverseCellWidth); + glm::vec3 grid_index_max = glm::floor((curr_boid_pos - gridMin + glm::vec3(neighbor_radius)) * inverseCellWidth); + + int max_x = imin(gridResolution - 1, grid_index_max.x); + int max_y = imin(gridResolution - 1, grid_index_max.y); + int max_z = imin(gridResolution - 1, grid_index_max.z); + + int min_x = imax(0, grid_index_min.x); + int min_y = imax(0, grid_index_min.y); + int min_z = imax(0, grid_index_min.z); // - For each cell, read the start/end indices in the boid pointer array. glm::vec3 perceived_center(0.f, 0.f, 0.f); glm::vec3 avoidance_velocity(0.f, 0.f, 0.f); @@ -596,7 +602,7 @@ __global__ void kernUpdateVelNeighborSearchScattered( { rule3_component = (perceived_velocity / (float)neighbor_count_rule3) * rule3Scale; } - new_velocity = vel1[idx] + rule1_component + rule2_component + rule3_component; + new_velocity += vel1[idx] + rule1_component + rule2_component + rule3_component; // - Clamp the speed change before putting the new speed in vel2 float curr_speed = glm::length(new_velocity); //if the total speed of vel is larger than maxSpeed, we normalize the vel and apply the maxSpeed we allow @@ -641,20 +647,17 @@ __global__ void kernUpdateVelNeighborSearchCoherent( int gridIndex = gridIndex3Dto1D(idx_x, idx_y, idx_z, gridResolution); // - Identify which cells may contain neighbors. This isn't always 8. //by calculating those cells that interact with the sphere with the neighbor_radius(std::max(std::max(rule1Distance, rule2Distance), rule3Distance)) - //float neighbor_radius = cellWidth / CELL_WIDTH_FACTOR; const float neighbor_radius = glm::max(glm::max(rule1Distance, rule2Distance), rule3Distance); - int max_x = (curr_boid_pos.x + neighbor_radius - gridMin.x) * inverseCellWidth; - max_x = max_x > gridResolution - 1 ? gridResolution - 1 : max_x; - int min_x = (curr_boid_pos.x - neighbor_radius - gridMin.x) * inverseCellWidth; - min_x = min_x < 0 ? 0 : min_x; - int max_y = (curr_boid_pos.y + neighbor_radius - gridMin.y) * inverseCellWidth; - max_y = max_y > gridResolution - 1 ? gridResolution - 1 : max_y; - int min_y = (curr_boid_pos.y + neighbor_radius - gridMin.y) * inverseCellWidth; - min_y = min_y < 0 ? 0 : min_y; - int max_z = (curr_boid_pos.z + neighbor_radius - gridMin.z) * inverseCellWidth; - max_z = max_z > gridResolution - 1 ? gridResolution - 1 : max_z; - int min_z = (curr_boid_pos.z + neighbor_radius - gridMin.z) * inverseCellWidth; - min_z = min_z < 0 ? 0 : min_z; + glm::vec3 grid_index_min = glm::floor((curr_boid_pos - gridMin - glm::vec3(neighbor_radius)) * inverseCellWidth); + glm::vec3 grid_index_max = glm::floor((curr_boid_pos - gridMin + glm::vec3(neighbor_radius)) * inverseCellWidth); + + int max_x = imin(gridResolution - 1, grid_index_max.x); + int max_y = imin(gridResolution - 1, grid_index_max.y); + int max_z = imin(gridResolution - 1, grid_index_max.z); + + int min_x = imax(0, grid_index_min.x); + int min_y = imax(0, grid_index_min.y); + int min_z = imax(0, grid_index_min.z); // - For each cell, read the start/end indices in the boid pointer array. glm::vec3 perceived_center(0.f, 0.f, 0.f); @@ -694,20 +697,20 @@ __global__ void kernUpdateVelNeighborSearchCoherent( glm::vec3 idx_boid_pos = pos[boid_idx]; float dist = glm::distance(idx_boid_pos, curr_boid_pos); // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves - if (dist < rule1Distance) + if (dist <= rule1Distance) { perceived_center += idx_boid_pos; neighbor_count_rule1++; } // Rule 2: boids try to stay a distance d away from each other - if (dist < rule2Distance) + if (dist <= rule2Distance) { avoidance_velocity -= (idx_boid_pos - curr_boid_pos); } // Rule 3: boids try to match the speed of surrounding boids - if (dist < rule3Distance) + if (dist <= rule3Distance) { perceived_velocity += vel1[boid_idx]; neighbor_count_rule3++; From cda89a242ff271b7e689169cdae7fd0ea371c800 Mon Sep 17 00:00:00 2001 From: Mark Date: Sun, 8 Sep 2019 03:17:58 -0400 Subject: [PATCH 7/8] add automation to generate performance data --- CMakeLists.txt | 1 + src/cxxopts.hpp | 2207 +++++++++++++++++++++++++++++++++++++++++++++++ src/main.cpp | 169 +++- src/main.hpp | 2 + 4 files changed, 2373 insertions(+), 6 deletions(-) create mode 100644 src/cxxopts.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 150664e..e250439 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -67,6 +67,7 @@ set(headers src/kernel.h src/main.hpp src/utilityCore.hpp + src/cxxopts.hpp ) set(sources diff --git a/src/cxxopts.hpp b/src/cxxopts.hpp new file mode 100644 index 0000000..ed3c6a2 --- /dev/null +++ b/src/cxxopts.hpp @@ -0,0 +1,2207 @@ +/* + +Copyright (c) 2014, 2015, 2016, 2017 Jarryd Beck + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. + +*/ + +#ifndef CXXOPTS_HPP_INCLUDED +#define CXXOPTS_HPP_INCLUDED + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef __cpp_lib_optional +#include +#define CXXOPTS_HAS_OPTIONAL +#endif + +#ifndef CXXOPTS_VECTOR_DELIMITER +#define CXXOPTS_VECTOR_DELIMITER ',' +#endif + +#define CXXOPTS__VERSION_MAJOR 2 +#define CXXOPTS__VERSION_MINOR 2 +#define CXXOPTS__VERSION_PATCH 0 + +namespace cxxopts +{ + static constexpr struct { + uint8_t major, minor, patch; + } version = { + CXXOPTS__VERSION_MAJOR, + CXXOPTS__VERSION_MINOR, + CXXOPTS__VERSION_PATCH + }; +} + +//when we ask cxxopts to use Unicode, help strings are processed using ICU, +//which results in the correct lengths being computed for strings when they +//are formatted for the help output +//it is necessary to make sure that can be found by the +//compiler, and that icu-uc is linked in to the binary. + +#ifdef CXXOPTS_USE_UNICODE +#include + +namespace cxxopts +{ + typedef icu::UnicodeString String; + + inline + String + toLocalString(std::string s) + { + return icu::UnicodeString::fromUTF8(std::move(s)); + } + + class UnicodeStringIterator : public + std::iterator + { + public: + + UnicodeStringIterator(const icu::UnicodeString* string, int32_t pos) + : s(string) + , i(pos) + { + } + + value_type + operator*() const + { + return s->char32At(i); + } + + bool + operator==(const UnicodeStringIterator& rhs) const + { + return s == rhs.s && i == rhs.i; + } + + bool + operator!=(const UnicodeStringIterator& rhs) const + { + return !(*this == rhs); + } + + UnicodeStringIterator& + operator++() + { + ++i; + return *this; + } + + UnicodeStringIterator + operator+(int32_t v) + { + return UnicodeStringIterator(s, i + v); + } + + private: + const icu::UnicodeString* s; + int32_t i; + }; + + inline + String& + stringAppend(String&s, String a) + { + return s.append(std::move(a)); + } + + inline + String& + stringAppend(String& s, int n, UChar32 c) + { + for (int i = 0; i != n; ++i) + { + s.append(c); + } + + return s; + } + + template + String& + stringAppend(String& s, Iterator begin, Iterator end) + { + while (begin != end) + { + s.append(*begin); + ++begin; + } + + return s; + } + + inline + size_t + stringLength(const String& s) + { + return s.length(); + } + + inline + std::string + toUTF8String(const String& s) + { + std::string result; + s.toUTF8String(result); + + return result; + } + + inline + bool + empty(const String& s) + { + return s.isEmpty(); + } +} + +namespace std +{ + inline + cxxopts::UnicodeStringIterator + begin(const icu::UnicodeString& s) + { + return cxxopts::UnicodeStringIterator(&s, 0); + } + + inline + cxxopts::UnicodeStringIterator + end(const icu::UnicodeString& s) + { + return cxxopts::UnicodeStringIterator(&s, s.length()); + } +} + +//ifdef CXXOPTS_USE_UNICODE +#else + +namespace cxxopts +{ + typedef std::string String; + + template + T + toLocalString(T&& t) + { + return std::forward(t); + } + + inline + size_t + stringLength(const String& s) + { + return s.length(); + } + + inline + String& + stringAppend(String&s, String a) + { + return s.append(std::move(a)); + } + + inline + String& + stringAppend(String& s, size_t n, char c) + { + return s.append(n, c); + } + + template + String& + stringAppend(String& s, Iterator begin, Iterator end) + { + return s.append(begin, end); + } + + template + std::string + toUTF8String(T&& t) + { + return std::forward(t); + } + + inline + bool + empty(const std::string& s) + { + return s.empty(); + } +} + +//ifdef CXXOPTS_USE_UNICODE +#endif + +namespace cxxopts +{ + namespace + { +#ifdef _WIN32 + const std::string LQUOTE("\'"); + const std::string RQUOTE("\'"); +#else + const std::string LQUOTE("‘"); + const std::string RQUOTE("’"); +#endif + } + + class Value : public std::enable_shared_from_this + { + public: + + virtual ~Value() = default; + + virtual + std::shared_ptr + clone() const = 0; + + virtual void + parse(const std::string& text) const = 0; + + virtual void + parse() const = 0; + + virtual bool + has_default() const = 0; + + virtual bool + is_container() const = 0; + + virtual bool + has_implicit() const = 0; + + virtual std::string + get_default_value() const = 0; + + virtual std::string + get_implicit_value() const = 0; + + virtual std::shared_ptr + default_value(const std::string& value) = 0; + + virtual std::shared_ptr + implicit_value(const std::string& value) = 0; + + virtual std::shared_ptr + no_implicit_value() = 0; + + virtual bool + is_boolean() const = 0; + }; + + class OptionException : public std::exception + { + public: + OptionException(const std::string& message) + : m_message(message) + { + } + + virtual const char* + what() const noexcept + { + return m_message.c_str(); + } + + private: + std::string m_message; + }; + + class OptionSpecException : public OptionException + { + public: + + OptionSpecException(const std::string& message) + : OptionException(message) + { + } + }; + + class OptionParseException : public OptionException + { + public: + OptionParseException(const std::string& message) + : OptionException(message) + { + } + }; + + class option_exists_error : public OptionSpecException + { + public: + option_exists_error(const std::string& option) + : OptionSpecException("Option " + LQUOTE + option + RQUOTE + " already exists") + { + } + }; + + class invalid_option_format_error : public OptionSpecException + { + public: + invalid_option_format_error(const std::string& format) + : OptionSpecException("Invalid option format " + LQUOTE + format + RQUOTE) + { + } + }; + + class option_syntax_exception : public OptionParseException { + public: + option_syntax_exception(const std::string& text) + : OptionParseException("Argument " + LQUOTE + text + RQUOTE + + " starts with a - but has incorrect syntax") + { + } + }; + + class option_not_exists_exception : public OptionParseException + { + public: + option_not_exists_exception(const std::string& option) + : OptionParseException("Option " + LQUOTE + option + RQUOTE + " does not exist") + { + } + }; + + class missing_argument_exception : public OptionParseException + { + public: + missing_argument_exception(const std::string& option) + : OptionParseException( + "Option " + LQUOTE + option + RQUOTE + " is missing an argument" + ) + { + } + }; + + class option_requires_argument_exception : public OptionParseException + { + public: + option_requires_argument_exception(const std::string& option) + : OptionParseException( + "Option " + LQUOTE + option + RQUOTE + " requires an argument" + ) + { + } + }; + + class option_not_has_argument_exception : public OptionParseException + { + public: + option_not_has_argument_exception + ( + const std::string& option, + const std::string& arg + ) + : OptionParseException( + "Option " + LQUOTE + option + RQUOTE + + " does not take an argument, but argument " + + LQUOTE + arg + RQUOTE + " given" + ) + { + } + }; + + class option_not_present_exception : public OptionParseException + { + public: + option_not_present_exception(const std::string& option) + : OptionParseException("Option " + LQUOTE + option + RQUOTE + " not present") + { + } + }; + + class argument_incorrect_type : public OptionParseException + { + public: + argument_incorrect_type + ( + const std::string& arg + ) + : OptionParseException( + "Argument " + LQUOTE + arg + RQUOTE + " failed to parse" + ) + { + } + }; + + class option_required_exception : public OptionParseException + { + public: + option_required_exception(const std::string& option) + : OptionParseException( + "Option " + LQUOTE + option + RQUOTE + " is required but not present" + ) + { + } + }; + + template + void throw_or_mimic(const std::string& text) + { + static_assert(std::is_base_of::value, + "throw_or_mimic only works on std::exception and " + "deriving classes"); + +#ifndef CXXOPTS_NO_EXCEPTIONS + // If CXXOPTS_NO_EXCEPTIONS is not defined, just throw + throw T{text}; +#else + // Otherwise manually instantiate the exception, print what() to stderr, + // and abort + T exception{text}; + std::cerr << exception.what() << std::endl; + std::cerr << "Aborting (exceptions disabled)..." << std::endl; + std::abort(); +#endif + } + + namespace values + { + namespace + { + std::basic_regex integer_pattern + ("(-)?(0x)?([0-9a-zA-Z]+)|((0x)?0)"); + std::basic_regex truthy_pattern + ("(t|T)(rue)?|1"); + std::basic_regex falsy_pattern + ("(f|F)(alse)?|0"); + } + + namespace detail + { + template + struct SignedCheck; + + template + struct SignedCheck + { + template + void + operator()(bool negative, U u, const std::string& text) + { + if (negative) + { + if (u > static_cast((std::numeric_limits::min)())) + { + throw_or_mimic(text); + } + } + else + { + if (u > static_cast((std::numeric_limits::max)())) + { + throw_or_mimic(text); + } + } + } + }; + + template + struct SignedCheck + { + template + void + operator()(bool, U, const std::string&) {} + }; + + template + void + check_signed_range(bool negative, U value, const std::string& text) + { + SignedCheck::is_signed>()(negative, value, text); + } + } + + template + R + checked_negate(T&& t, const std::string&, std::true_type) + { + // if we got to here, then `t` is a positive number that fits into + // `R`. So to avoid MSVC C4146, we first cast it to `R`. + // See https://github.com/jarro2783/cxxopts/issues/62 for more details. + return -static_cast(t-1)-1; + } + + template + T + checked_negate(T&& t, const std::string& text, std::false_type) + { + throw_or_mimic(text); + return t; + } + + template + void + integer_parser(const std::string& text, T& value) + { + std::smatch match; + std::regex_match(text, match, integer_pattern); + + if (match.length() == 0) + { + throw_or_mimic(text); + } + + if (match.length(4) > 0) + { + value = 0; + return; + } + + using US = typename std::make_unsigned::type; + + constexpr bool is_signed = std::numeric_limits::is_signed; + const bool negative = match.length(1) > 0; + const uint8_t base = match.length(2) > 0 ? 16 : 10; + + auto value_match = match[3]; + + US result = 0; + + for (auto iter = value_match.first; iter != value_match.second; ++iter) + { + US digit = 0; + + if (*iter >= '0' && *iter <= '9') + { + digit = static_cast(*iter - '0'); + } + else if (base == 16 && *iter >= 'a' && *iter <= 'f') + { + digit = static_cast(*iter - 'a' + 10); + } + else if (base == 16 && *iter >= 'A' && *iter <= 'F') + { + digit = static_cast(*iter - 'A' + 10); + } + else + { + throw_or_mimic(text); + } + + US next = result * base + digit; + if (result > next) + { + throw_or_mimic(text); + } + + result = next; + } + + detail::check_signed_range(negative, result, text); + + if (negative) + { + value = checked_negate(result, + text, + std::integral_constant()); + } + else + { + value = static_cast(result); + } + } + + template + void stringstream_parser(const std::string& text, T& value) + { + std::stringstream in(text); + in >> value; + if (!in) { + throw_or_mimic(text); + } + } + + inline + void + parse_value(const std::string& text, uint8_t& value) + { + integer_parser(text, value); + } + + inline + void + parse_value(const std::string& text, int8_t& value) + { + integer_parser(text, value); + } + + inline + void + parse_value(const std::string& text, uint16_t& value) + { + integer_parser(text, value); + } + + inline + void + parse_value(const std::string& text, int16_t& value) + { + integer_parser(text, value); + } + + inline + void + parse_value(const std::string& text, uint32_t& value) + { + integer_parser(text, value); + } + + inline + void + parse_value(const std::string& text, int32_t& value) + { + integer_parser(text, value); + } + + inline + void + parse_value(const std::string& text, uint64_t& value) + { + integer_parser(text, value); + } + + inline + void + parse_value(const std::string& text, int64_t& value) + { + integer_parser(text, value); + } + + inline + void + parse_value(const std::string& text, bool& value) + { + std::smatch result; + std::regex_match(text, result, truthy_pattern); + + if (!result.empty()) + { + value = true; + return; + } + + std::regex_match(text, result, falsy_pattern); + if (!result.empty()) + { + value = false; + return; + } + + throw_or_mimic(text); + } + + inline + void + parse_value(const std::string& text, std::string& value) + { + value = text; + } + + // The fallback parser. It uses the stringstream parser to parse all types + // that have not been overloaded explicitly. It has to be placed in the + // source code before all other more specialized templates. + template + void + parse_value(const std::string& text, T& value) { + stringstream_parser(text, value); + } + + template + void + parse_value(const std::string& text, std::vector& value) + { + std::stringstream in(text); + std::string token; + while(in.eof() == false && std::getline(in, token, CXXOPTS_VECTOR_DELIMITER)) { + T v; + parse_value(token, v); + value.emplace_back(std::move(v)); + } + } + +#ifdef CXXOPTS_HAS_OPTIONAL + template + void + parse_value(const std::string& text, std::optional& value) + { + T result; + parse_value(text, result); + value = std::move(result); + } +#endif + + inline + void parse_value(const std::string& text, char& c) + { + if (text.length() != 1) + { + throw_or_mimic(text); + } + + c = text[0]; + } + + template + struct type_is_container + { + static constexpr bool value = false; + }; + + template + struct type_is_container> + { + static constexpr bool value = true; + }; + + template + class abstract_value : public Value + { + using Self = abstract_value; + + public: + abstract_value() + : m_result(std::make_shared()) + , m_store(m_result.get()) + { + } + + abstract_value(T* t) + : m_store(t) + { + } + + virtual ~abstract_value() = default; + + abstract_value(const abstract_value& rhs) + { + if (rhs.m_result) + { + m_result = std::make_shared(); + m_store = m_result.get(); + } + else + { + m_store = rhs.m_store; + } + + m_default = rhs.m_default; + m_implicit = rhs.m_implicit; + m_default_value = rhs.m_default_value; + m_implicit_value = rhs.m_implicit_value; + } + + void + parse(const std::string& text) const + { + parse_value(text, *m_store); + } + + bool + is_container() const + { + return type_is_container::value; + } + + void + parse() const + { + parse_value(m_default_value, *m_store); + } + + bool + has_default() const + { + return m_default; + } + + bool + has_implicit() const + { + return m_implicit; + } + + std::shared_ptr + default_value(const std::string& value) + { + m_default = true; + m_default_value = value; + return shared_from_this(); + } + + std::shared_ptr + implicit_value(const std::string& value) + { + m_implicit = true; + m_implicit_value = value; + return shared_from_this(); + } + + std::shared_ptr + no_implicit_value() + { + m_implicit = false; + return shared_from_this(); + } + + std::string + get_default_value() const + { + return m_default_value; + } + + std::string + get_implicit_value() const + { + return m_implicit_value; + } + + bool + is_boolean() const + { + return std::is_same::value; + } + + const T& + get() const + { + if (m_store == nullptr) + { + return *m_result; + } + else + { + return *m_store; + } + } + + protected: + std::shared_ptr m_result; + T* m_store; + + bool m_default = false; + bool m_implicit = false; + + std::string m_default_value; + std::string m_implicit_value; + }; + + template + class standard_value : public abstract_value + { + public: + using abstract_value::abstract_value; + + std::shared_ptr + clone() const + { + return std::make_shared>(*this); + } + }; + + template <> + class standard_value : public abstract_value + { + public: + ~standard_value() = default; + + standard_value() + { + set_default_and_implicit(); + } + + standard_value(bool* b) + : abstract_value(b) + { + set_default_and_implicit(); + } + + std::shared_ptr + clone() const + { + return std::make_shared>(*this); + } + + private: + + void + set_default_and_implicit() + { + m_default = true; + m_default_value = "false"; + m_implicit = true; + m_implicit_value = "true"; + } + }; + } + + template + std::shared_ptr + value() + { + return std::make_shared>(); + } + + template + std::shared_ptr + value(T& t) + { + return std::make_shared>(&t); + } + + class OptionAdder; + + class OptionDetails + { + public: + OptionDetails + ( + const std::string& short_, + const std::string& long_, + const String& desc, + std::shared_ptr val + ) + : m_short(short_) + , m_long(long_) + , m_desc(desc) + , m_value(val) + , m_count(0) + { + } + + OptionDetails(const OptionDetails& rhs) + : m_desc(rhs.m_desc) + , m_count(rhs.m_count) + { + m_value = rhs.m_value->clone(); + } + + OptionDetails(OptionDetails&& rhs) = default; + + const String& + description() const + { + return m_desc; + } + + const Value& value() const { + return *m_value; + } + + std::shared_ptr + make_storage() const + { + return m_value->clone(); + } + + const std::string& + short_name() const + { + return m_short; + } + + const std::string& + long_name() const + { + return m_long; + } + + private: + std::string m_short; + std::string m_long; + String m_desc; + std::shared_ptr m_value; + int m_count; + }; + + struct HelpOptionDetails + { + std::string s; + std::string l; + String desc; + bool has_default; + std::string default_value; + bool has_implicit; + std::string implicit_value; + std::string arg_help; + bool is_container; + bool is_boolean; + }; + + struct HelpGroupDetails + { + std::string name; + std::string description; + std::vector options; + }; + + class OptionValue + { + public: + void + parse + ( + std::shared_ptr details, + const std::string& text + ) + { + ensure_value(details); + ++m_count; + m_value->parse(text); + } + + void + parse_default(std::shared_ptr details) + { + ensure_value(details); + m_default = true; + m_value->parse(); + } + + size_t + count() const noexcept + { + return m_count; + } + + // TODO: maybe default options should count towards the number of arguments + bool + has_default() const noexcept + { + return m_default; + } + + template + const T& + as() const + { + if (m_value == nullptr) { + throw_or_mimic("No value"); + } + +#ifdef CXXOPTS_NO_RTTI + return static_cast&>(*m_value).get(); +#else + return dynamic_cast&>(*m_value).get(); +#endif + } + + private: + void + ensure_value(std::shared_ptr details) + { + if (m_value == nullptr) + { + m_value = details->make_storage(); + } + } + + std::shared_ptr m_value; + size_t m_count = 0; + bool m_default = false; + }; + + class KeyValue + { + public: + KeyValue(std::string key_, std::string value_) + : m_key(std::move(key_)) + , m_value(std::move(value_)) + { + } + + const + std::string& + key() const + { + return m_key; + } + + const + std::string& + value() const + { + return m_value; + } + + template + T + as() const + { + T result; + values::parse_value(m_value, result); + return result; + } + + private: + std::string m_key; + std::string m_value; + }; + + class ParseResult + { + public: + + ParseResult( + const std::shared_ptr< + std::unordered_map> + >, + std::vector, + bool allow_unrecognised, + int&, char**&); + + size_t + count(const std::string& o) const + { + auto iter = m_options->find(o); + if (iter == m_options->end()) + { + return 0; + } + + auto riter = m_results.find(iter->second); + + return riter->second.count(); + } + + const OptionValue& + operator[](const std::string& option) const + { + auto iter = m_options->find(option); + + if (iter == m_options->end()) + { + throw_or_mimic(option); + } + + auto riter = m_results.find(iter->second); + + return riter->second; + } + + const std::vector& + arguments() const + { + return m_sequential; + } + + private: + + void + parse(int& argc, char**& argv); + + void + add_to_option(const std::string& option, const std::string& arg); + + bool + consume_positional(std::string a); + + void + parse_option + ( + std::shared_ptr value, + const std::string& name, + const std::string& arg = "" + ); + + void + parse_default(std::shared_ptr details); + + void + checked_parse_arg + ( + int argc, + char* argv[], + int& current, + std::shared_ptr value, + const std::string& name + ); + + const std::shared_ptr< + std::unordered_map> + > m_options; + std::vector m_positional; + std::vector::iterator m_next_positional; + std::unordered_set m_positional_set; + std::unordered_map, OptionValue> m_results; + + bool m_allow_unrecognised; + + std::vector m_sequential; + }; + + struct Option + { + Option + ( + const std::string& opts, + const std::string& desc, + const std::shared_ptr& value = ::cxxopts::value(), + const std::string& arg_help = "" + ) + : opts_(opts) + , desc_(desc) + , value_(value) + , arg_help_(arg_help) + { + } + + std::string opts_; + std::string desc_; + std::shared_ptr value_; + std::string arg_help_; + }; + + class Options + { + typedef std::unordered_map> + OptionMap; + public: + + Options(std::string program, std::string help_string = "") + : m_program(std::move(program)) + , m_help_string(toLocalString(std::move(help_string))) + , m_custom_help("[OPTION...]") + , m_positional_help("positional parameters") + , m_show_positional(false) + , m_allow_unrecognised(false) + , m_options(std::make_shared()) + , m_next_positional(m_positional.end()) + { + } + + Options& + positional_help(std::string help_text) + { + m_positional_help = std::move(help_text); + return *this; + } + + Options& + custom_help(std::string help_text) + { + m_custom_help = std::move(help_text); + return *this; + } + + Options& + show_positional_help() + { + m_show_positional = true; + return *this; + } + + Options& + allow_unrecognised_options() + { + m_allow_unrecognised = true; + return *this; + } + + ParseResult + parse(int& argc, char**& argv); + + OptionAdder + add_options(std::string group = ""); + + void + add_options + ( + const std::string& group, + std::initializer_list