code: def _gsa_force(grav, mass_i, mass_j, position_i, position_j):
      
  
      position_diff = numpy.subtract(position_j, position_i)
      distance = numpy.linalg.norm(position_diff)
  
      # The first 3 terms give the magnitude of the force
      # The last term is a vector that provides the direction
      # Epsilon prevents divide by zero errors
      return grav * (mass_i * mass_j) / (distance + EPSILON) * position_diff
docs: Gives the force of solution j on solution i.
  
      Variable name in GSA paper given in ()
  
      args:
          grav: The gravitational constant. (G)
          mass_i: The mass of solution i (derived from fitness). (M_i)
          mass_j: The mass of solution j (derived from fitness). (M_j)
          position_i: The position of solution i. (x_i)
          position_j: The position of solution j. (x_j)
  
      returns:
          numpy.array; The force vector of solution j on solution i.