docs: def _new_population_gsa(population, fitnesses, velocities, lower_bounds,
                          upper_bounds, grav_initial, grav_reduction_rate,
                          iteration, max_iterations):
      """Generate a new population as given by GSA algorithm.
  
      In GSA paper, grav_initial is G_i
      """
      # Update the gravitational constant, and the best and worst of the population
      # Calculate the mass and acceleration for each solution
      # Update the velocity and position of each solution
      population_size = len(population)
      solution_size = len(population[0])
  
      # In GSA paper, grav is G
      grav = _next_grav_gsa(grav_initial, grav_reduction_rate, iteration,
                            max_iterations)
      masses = _get_masses(fitnesses)
  
      # Create bundled solution with position and mass for the K best calculation
      # Also store index to later check if two solutions are the same
      # Sorted by solution fitness (mass)
      solutions = [{
          'pos': pos,
          'mass': mass,
          'index': i
      } for i, (pos, mass) in enumerate(zip(population, masses))]
      solutions.sort(key=lambda x: x['mass'], reverse=True)
  
      # Get the force on each solution
      # Only the best K solutions apply force
      # K linearly decreases to 1
      num_best = int(population_size - (population_size - 1) *
                     (iteration / float(max_iterations)))
      forces = []
      for i in range(population_size):
          force_vectors = []
          for j in range(num_best):
              # If it is not the same solution
              if i != solutions[j]['index']:
                  force_vectors.append(
                      _gsa_force(grav, masses[i], solutions[j]['mass'],
                                 population[i], solutions[j]['pos']))
          forces.append(_gsa_total_force(force_vectors, solution_size))
  
      # Get the acceleration of each solution
      accelerations = []
      for i in range(population_size):
          accelerations.append(_gsa_acceleration(forces[i], masses[i]))
  
      # Update the velocity of each solution
      new_velocities = []
      for i in range(population_size):
          new_velocities.append(
              _gsa_update_velocity(velocities[i], accelerations[i]))
  
      # Create the new population
      new_population = []
      for i in range(population_size):
          new_position = _gsa_update_position(population[i], new_velocities[i])
          # Constrain to bounds
          new_position = list(
              numpy.clip(new_position, lower_bounds, upper_bounds))
  
          new_population.append(new_position)
  
      return new_population, new_velocities