stevengiacalone / triceratops

Tool for Rating Interesting Candidate Exoplanets and Reliability Analysis of Transits Originating from Proximate Stars
https://triceratops.readthedocs.io
MIT License
20 stars 7 forks source link

TypingError in triceratops.target.calc_probs #15

Closed Christopher-Mann closed 1 year ago

Christopher-Mann commented 2 years ago

Hi there. I'm newly trying out triceratops to help with a TESS paper I'm finalizing (and need the kind of statistical validation that Triceratops can provide).

I finally got all the modules to load properly in a dedicated conda environemnt and I swapped the shell to bash so that lightkurve will run (Catalina seems to default to a zsh shell which causes issues). I'm working through the TESS tutorial notebook and hit the following issue in the 5th code cell, where we first call: target.calc_probs(time=lc.time.value, flux_0=lc.flux.value, flux_err_0=np.mean(lc.flux_err.value), P_orb=P_orb)

I get the following error: (I haven't edited the tutorial file at all)

---------------------------------------------------------------------------
TypingError                               Traceback (most recent call last)
File <timed exec>:9, in <module>

File /Applications/anaconda/envs/triceratops/lib/python3.8/site-packages/triceratops/triceratops.py:731, in target.calc_probs(self, time, flux_0, flux_err_0, P_orb, contrast_curve_file, filt, N, parallel, drop_scenario, verbose, flatpriors, exptime, nsamples, molusc_file, trilegal_fname)
    725 if verbose == 1:
    726     print(
    727         "Calculating TP scenario "
    728         + "probabilitiey for " + str(ID) + "."
    729         )
--> 731 res = lnZ_TTP(
    732     time, flux, flux_err, P_orb,
    733     M_s, R_s, Teff, Z,
    734     N, parallel, self.mission,
    735     flatpriors,
    736     exptime, nsamples
    737     )
    738 # self.res_TTP = res
    739 j = 0

File /Applications/anaconda/envs/triceratops/lib/python3.8/site-packages/triceratops/marginal_likelihoods.py:144, in lnZ_TTP(time, flux, sigma, P_orb, M_s, R_s, Teff, Z, N, parallel, mission, flatpriors, exptime, nsamples)
    142             continue
    143         if (incs[i] >= inc_min) & (coll[i] == False):
--> 144             lnL[i] = -0.5*ln2pi - lnsigma - lnL_TP(
    145                 time, flux, sigma, rps[i],
    146                 P_orb[i], incs[i], a, R_s, u1, u2,
    147                 eccs[i], argps[i],
    148                 exptime=exptime, nsamples=nsamples
    149                 )
    151 N_samples = 100
    152 idx = (-lnL).argsort()[:N_samples]

File /Applications/anaconda/envs/triceratops/lib/python3.8/site-packages/triceratops/likelihoods.py:186, in lnL_TP(time, flux, sigma, R_p, P_orb, inc, a, R_s, u1, u2, ecc, argp, companion_fluxratio, companion_is_host, exptime, nsamples)
    152 def lnL_TP(time: np.ndarray, flux: np.ndarray, sigma: float, R_p: float,
    153            P_orb: float, inc: float, a: float, R_s: float,
    154            u1: float, u2: float, ecc: float, argp: float,
   (...)
    157            exptime: float = 0.00139,
    158            nsamples: int = 20):
    159     """
    160     Calculates the log likelihood of a transiting planet scenario by
    161     comparing a simulated light curve and the TESS light curve.
   (...)
    184         Log likelihood (float).
    185     """
--> 186     model = simulate_TP_transit(
    187         time, R_p, P_orb, inc, a, R_s, u1, u2,
    188         ecc, argp,
    189         companion_fluxratio, companion_is_host,
    190         exptime, nsamples
    191         )
    192     return 0.5*(np.sum((flux-model)**2 / sigma**2))

File /Applications/anaconda/envs/triceratops/lib/python3.8/site-packages/triceratops/likelihoods.py:50, in simulate_TP_transit(time, R_p, P_orb, inc, a, R_s, u1, u2, ecc, argp, companion_fluxratio, companion_is_host, exptime, nsamples)
     48 # step 1: simulate light curve assuming only the host star exists
     49 tm.set_data(time, exptimes=exptime, nsamples=nsamples)
---> 50 flux = tm.evaluate_ps(
     51     k=R_p*Rearth/(R_s*Rsun),
     52     ldc=[float(u1), float(u2)],
     53     t0=0.0,
     54     p=P_orb,
     55     a=a/(R_s*Rsun),
     56     i=inc*(pi/180.),
     57     e=ecc,
     58     w=(90-argp)*(pi/180.)
     59     )
     60 # step 2: adjust the light curve to account for flux dilution
     61 # from non-host star
     62 if companion_is_host:

File /Applications/anaconda/envs/triceratops/lib/python3.8/site-packages/pytransit/models/ma_quadratic.py:183, in QuadraticModel.evaluate_ps(self, k, ldc, t0, p, a, i, e, w, copy)
    180 if ldc.size != 2 * self.npb:
    181     raise ValueError("The quadratic model needs two limb darkening coefficients per passband")
--> 183 flux = quadratic_model_s(self.time, k, t0, p, a, i, e, w, ldc,
    184                          self.lcids, self.pbids, self.epids, self.nsamples, self.exptimes, self.npb,
    185                          self.ed, self.ld, self.le, self.kt, self.zt, self.interpolate)
    186 return squeeze(flux)

File /Applications/anaconda/envs/triceratops/lib/python3.8/site-packages/numba/core/dispatcher.py:468, in _DispatcherBase._compile_for_args(self, *args, **kws)
    464         msg = (f"{str(e).rstrip()} \n\nThis error may have been caused "
    465                f"by the following argument(s):\n{args_str}\n")
    466         e.patch_message(msg)
--> 468     error_rewrite(e, 'typing')
    469 except errors.UnsupportedError as e:
    470     # Something unsupported is present in the user code, add help info
    471     error_rewrite(e, 'unsupported_error')

File /Applications/anaconda/envs/triceratops/lib/python3.8/site-packages/numba/core/dispatcher.py:409, in _DispatcherBase._compile_for_args.<locals>.error_rewrite(e, issue_type)
    407     raise e
    408 else:
--> 409     raise e.with_traceback(None)

TypingError: Failed in nopython mode pipeline (step: nopython frontend)
Failed in nopython mode pipeline (step: nopython frontend)
Failed in nopython mode pipeline (step: nopython frontend)
Cannot unify float64 and array(float64, 1d, C) for 't0.2', defined at /Applications/anaconda/envs/triceratops/lib/python3.8/site-packages/pytransit/orbits/taylor_z.py (320)

File "../../../../../../../../Applications/anaconda/envs/triceratops/lib/python3.8/site-packages/pytransit/orbits/taylor_z.py", line 320:
def find_contact_point(k: float, point: int, y0, vx, vy, ax, ay, jx, jy, sx, sy):
    <source elided>
    i = 0
    while abs(t2 - t0) > 1e-6 and i < 100:
    ^

During: typing of assignment at /Applications/anaconda/envs/triceratops/lib/python3.8/site-packages/pytransit/orbits/taylor_z.py (320)

File "../../../../../../../../Applications/anaconda/envs/triceratops/lib/python3.8/site-packages/pytransit/orbits/taylor_z.py", line 320:
def find_contact_point(k: float, point: int, y0, vx, vy, ax, ay, jx, jy, sx, sy):
    <source elided>
    i = 0
    while abs(t2 - t0) > 1e-6 and i < 100:
    ^

During: resolving callee type: type(CPUDispatcher(<function find_contact_point at 0x7fac899713a0>))
During: typing of call at /Applications/anaconda/envs/triceratops/lib/python3.8/site-packages/pytransit/orbits/taylor_z.py (360)

During: resolving callee type: type(CPUDispatcher(<function find_contact_point at 0x7fac899713a0>))
During: typing of call at /Applications/anaconda/envs/triceratops/lib/python3.8/site-packages/pytransit/orbits/taylor_z.py (360)

File "../../../../../../../../Applications/anaconda/envs/triceratops/lib/python3.8/site-packages/pytransit/orbits/taylor_z.py", line 360:
def t14(k: float, y0, vx, vy, ax, ay, jx, jy, sx, sy):
    t1 = find_contact_point(k, 1, y0, vx, vy, ax, ay, jx, jy, sx, sy)
    ^

During: resolving callee type: type(CPUDispatcher(<function t14 at 0x7fac8995cdc0>))
During: typing of call at /Applications/anaconda/envs/triceratops/lib/python3.8/site-packages/pytransit/models/numba/ma_quadratic_nb.py (611)

During: resolving callee type: type(CPUDispatcher(<function t14 at 0x7fac8995cdc0>))
During: typing of call at /Applications/anaconda/envs/triceratops/lib/python3.8/site-packages/pytransit/models/numba/ma_quadratic_nb.py (611)

File "../../../../../../../../Applications/anaconda/envs/triceratops/lib/python3.8/site-packages/pytransit/models/numba/ma_quadratic_nb.py", line 611:
def quadratic_model_s(t, k, t0, p, a, i, e, w, ldc, lcids, pbids, epids, nsamples, exptimes, npb, edt, ldt, let, kt, zt, interpolate):
    <source elided>
    y0, vx, vy, ax, ay, jx, jy, sx, sy = vajs_from_paiew(p, a, i, e, w)
    half_window_width = 0.025 + 0.5 * t14(k[0], y0, vx, vy, ax, ay, jx, jy, sx, sy)
    ^

Can't say I really know what's going on, but it seems pretty deep inside the triceratops source code. Any insight on how to fix this?

Christopher-Mann commented 2 years ago

Update on my previous post:

I pressed on with the TESS tutorial and, for whatever reason, Examples 2 and 3 run just fine. Whatever the issue is, seems to be specific to Example 1.

I was able to apply the structure of examples 2 and 3 to my own analysis and things seem to be working fine.

stevengiacalone commented 2 years ago

The error seems to be an issue with pytransit, which is used to model the light curves. Specifically it seems that the non-parallelized LC calculation is failing, but the parallelized LC calculation is not. I'll need to do some troubleshooting to fix this, but for now setting "parallel = True" in calc_probs should be a safe bet.

supremeKAI40 commented 2 years ago

Hello there, I happened to face the same issue. I think this is being caused by recent version of numba (v0.56.x) released on 2022-07-25. It seems to raise data type conflicts somewhere inside pytransit.

Although downgraded version of triceratops (1.0.15) works just fine..

Not sure on the cause but I hope it helps

hpparvi commented 2 years ago

Hi, I'll see if I can fix this in PyTransit asap (1-2 days).

hpparvi commented 2 years ago

Hi,

Sorry for taking a bit longer before I could get to this. This is a problem in Triceratops marginal_likelihoods.py:

The orbital periord P_orb is not a float but an array in marginal_likelihoods.py lines 66-71, and because of this, also the semi-major axis a in line 74 is an array. The crash comes when the transit model is tried to evaluate for an array a instead of float.

I tested that changing the lines 147 and 568 (the numbers might be off by two or three since I've added some debugging breakpoints to my version) in marginal_likelihoods.py

P_orb[i], incs[i], a, R_s, u1, u2,

into

P_orb[i], incs[i], a[i], R_s, u1, u2,

seems to fix the issue. However, I think @stevengiacalone needs to take a closer look at this to make sure the bug is not repeated elsewhere in the code. I believe a is expected to be a float while a_arr is expected to be an ndarray, but at the moment both are arrays.

The parallelised code works because there PyTransit expects to get the parameters as arrays, but the non-parallelized code breaks because then all the orbital parameters should be floats.

stevengiacalone commented 2 years ago

Thanks for finding the root of the error, @hpparvi ! This is due to an update I made a few months back and should be easily fixable. I'll upload a new patched version of the code as soon as I'm able to.

stevengiacalone commented 1 year ago

Fixed in latest update.