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
Generate a new population as given by GSA algorithm.

    In GSA paper, grav_initial is G_i