Issue with equation of state "GS1" in XLALSimNeutronStarRadius()

While using lalsim.SimNeutronStarRadius (with lalsim version 3.1.0) I am getting a wrong shape in the plane m_ns vs r_ns for the EOS called "GS1". Here is the code I'm using and the plot I am getting:

import numpy as np
import random
import matplotlib
import matplotlib.pyplot as plt
import lalsimulation as lalsim
from lal import MSUN_SI

fig = plt.figure(num=None, figsize=(12, 7))
eos = lalsim.SimNeutronStarEOSByName("GS1")
fam = lalsim.CreateSimNeutronStarFamily(eos)
min_mass = lalsim.SimNeutronStarFamMinimumMass(fam)
max_mass = lalsim.SimNeutronStarMaximumMass(fam)
radius, mass = [], []
for i in range(0, 10000):
    m = random.uniform(min_mass, max_mass)
    mass.append(m/MSUN_SI)
    radius.append(lalsim.SimNeutronStarRadius(m, fam)*(0.001))
p0 = plt.scatter(radius, mass, s=1) 
plt.tight_layout()

I am comparing this trend with the one in Fig.2(on the right) in this paper and the two plots do not match.