Issues when testing stardate

[cpiaulet]

I installed stardate on a new python=3.8 conda environment, following the instructions.

The following example MCMC fits from the readthedocs run fine (no error), but I am not getting the right results:

My code:

import stardate as sd

# Create a dictionary of observables
iso_params = {"teff": (4386, 50),     # Teff with uncertainty.
              "logg": (4.66, .05),    # logg with uncertainty.
              "feh": (0.0, .02),      # Metallicity with uncertainty.
              "parallax": (1.48, .1),  # Parallax in milliarcseconds.
              "maxAV": .1}            # Maximum extinction

prot, prot_err = 29, 3

# Set up the star object.
star = sd.Star(iso_params, prot=prot, prot_err=prot_err)  # Here's where you add a rotation period

# Run the MCMC
star.fit(max_n=1000)

# max_n is the maximum number of MCMC samples. I recommend setting this
# much higher when running for real, or using the default value of 100000.

# Print the median age with the 16th and 84th percentile uncertainties.
age, errp, errm, samples = star.age_results()
print("stellar age = {0:.2f} + {1:.2f} + {2:.2f}".format(age, errp, errm))

I get as an output 0.42 + 0.32 + 3.37 rather than 2.97 + 0.60 + 0.55, and the bias towards low stellar ages remains even for larger numbers of steps (10,000).

When directly calling the Praesepe model to get period from age or age from period however, I am getting the expected values.

Another issue is that I am getting 3 failed tests when running pytest (see below).

platform darwin -- Python 3.8.13, pytest-7.1.2, pluggy-1.0.0
rootdir: /Users/caroline/Research/GitHub/stardate
collected 13 items                                                                                    

test/isochrone_test.py F                                                                        [  7%]
test/lhf_test.py ..F..F......                                                                   [100%]

============================================== FAILURES ===============================================
___________________________________________ test_iso_lnlike ___________________________________________

    def test_iso_lnlike():
    
        iso_params = {"teff": (5770, 10),
                    "feh": (0, .01),
                    "logg": (4.44, .1),
                    "parallax": (1, 1)}
    
        mod = SingleStarModel(mist, **iso_params)  # StarModel isochrones obj
        params = [354, np.log10(4.56*1e9), 0., 1000, 0.]
        lnparams = [354, np.log10(4.56*1e9), 0., np.log(1000), 0.]
        lnpr = mod.lnprior(params)
        # lnparams = [350, 9, 0., 6, 0.]
        args1 = [mod, iso_params]
    
        # Calculate the lnprob above
        iso_lp = iso_lnprob(lnparams, *args1)
        start = time.time()
        for i in range(100):
            iso_lp = iso_lnprob(lnparams, *args1)
        end = time.time()
    
        # Calculate the stardate lnprob
        args2 = [mod, None, None, True, False, True, "praesepe"]
    
        start = time.time()
        for i in range(100):
            lp = lnprob(lnparams, *args2)[0]
        end = time.time()
        # print("time = ", end - start)
    
        ll = lnlike(lnparams, *args2)
        lnprior = lp - ll
    
        assert iso_lp == lp
        assert np.isclose(iso_lp - lnpr, ll)
>       assert lnpr == lnprior
E       assert -18.46720822861122 == -18.46720822861107

test/isochrone_test.py:62: AssertionError
____________________________________________ test_calc_bv _____________________________________________

    def test_calc_bv():
        sun = [355, np.log10(4.56*1e9), 0., np.log(1000), 0.]
        bv_sun = calc_bv(sun)
>       assert .6 < bv_sun, "Are you sure you're on the isochrones eep branch?"
E       AssertionError: Are you sure you're on the isochrones eep branch?
E       assert 0.6 < 0.42834431703148956

test/lhf_test.py:92: AssertionError
_______________________________________ test_gyro_model_rossby ________________________________________

    def test_gyro_model_rossby():
        """
        Make sure that the rossby model gives Solar values for the Sun.
        Make sure it gives a maximum rotation period of pmax.
        """
        age = np.log10(4.56*1e9)
        sun = [355, age, 0., np.log(1000), 0.]
    
        # iso_params = {"teff": (5777, 10),
        #               "logg": (4.44, .05),
        #               "feh": (0., .001),
        #               "parallax": (1., .01),  # milliarcseconds
        #               "B": (15.48, 0.02)}
    
        # # Set up the StarModel isochrones object.
        # mod = StarModel(mist, **iso_params)
        # # the lnprob arguments]
        # args = [mod, 26., 1., False, True, "praesepe"]
        # print(lnprob(sun, *args))
        # assert 0
    
        prot_sun = 10**gyro_model_rossby(sun, model="praesepe")[0]
>       assert 21 < prot_sun
E       assert 21 < 6.909888089400217

test/lhf_test.py:142: AssertionError