JKTEBOP backend RV semi-amplitudes not working with mass ratio != 1

[jsinkbaek]

Minimum working example included as a jupyter notebook in attachment. rv_jktebop_example.zip

Hi!

I ran across this issue when I wanted to make a forward model in PHOEBE using parameters from my (external) JKTEBOP fit. I tried using the JKTEBOP (internal) backend to compare, and it did not match very well at all.

I can reproduce this problem by starting with q=M2/M1=1, where the PHOEBE and JKTEBOP backend produces equivalent RV curves. Afterwards, when altering with q=0.5, the PHOEBE model is changed. But the JKTEBOP model produces the exact same result as before.

It looks like the relevant line is: https://github.com/phoebe-project/phoebe2/blob/c4e2e30d5222a01fc9a8bc2831145093841d09d4/phoebe/backend/backends.py#L2340

As far as I can see, calculation of the RV semi-amplitudes is using the system 'sma@binary' in the JKTEBOP backend, instead of 'sma@primary' and 'sma@secondary'?

[jsinkbaek]

I have jury-rigged my local package to correct this issue and tested it.

It looks like changes are only needed in the phoebe2/phoebe/backend/backends.py file.

Below, I have gone through the ones I made.

In class JktebopBackend, method _workersetup. Insert new lines 2148-2149: sma_A = b.get_value(qualifier='sma', component=starrefs[0], context='component', unit=u.solRad, **_skip_filter_checks) sma_B = b.get_value(qualifier='sma', component=starrefs[1], context='component', unit=u.solRad, **_skip_filter_checks)

Then, still in _workersetup, insert following new lines in the return dictionary: sma_A=sma_A, sma_B=sma_B

Now, in _run_single_dataset(), insert new lines (now ~2191-2192): sma_A = kwargs.get('sma_A') sma_B = kwargs.get('sma_B')

In the if statement if info['kind'] == 'rv': insert a new first line: sma_ = sma_A if info['component'] == starrefs[0] else sma_B

Finally, in the same if statement, modify the line calculating K to be: K = np.pi * 2*(sma_*u.solRad).to(u.km).value * np.sin((incl*u.deg).to(u.rad).value) / (period*u.d).to(u.s).value

The last plot in the notebook should now show the change in q for the JKTEBOP backend as well:

[jsinkbaek]

After performing this change, it still does not reproduce my results for my actual stellar observations compared with using a model generated with JKTEBOP outside of PHOEBE. See below figure that compares residuals.

To me, it looks from both figures (synthetic example and star) that the JKTEBOP backend in PHOEBE in general produces RV semi-amplitudes that are a tiny bit smaller than what the PHOEBE backend, and external JKTEBOP, gets.

I have tried illustrating this with the synthetic example through a similar plot, where it is much more clear:

[kecnry]

Thanks for the detailed report! Can you compare against phoebe's dynamical RVs (rv_method='dynamical') to see if that improves the remaining small discrepancy? Otherwise, I think it still might be worth getting the immediate fix out and looking into the discrepancy later (it could be a bug in phoebe, a bug in jktebop, a bug in the translation layer, or just a small difference between the two codes). Would you be interested in opening a PR for a bugfix (I'm happy to help get you setup to do that)?

[jsinkbaek]

I will also remark that this small discrepancy seems to be present already before my hotfix, see the notebook.

[jsinkbaek]

Hi @kecnry,

Is what you are requesting what I have in the edited last message https://github.com/phoebe-project/phoebe2/issues/699#issuecomment-1429872347 ?

[jsinkbaek]

To clarify, the PHOEBE model in the synthetic example was done using rv_method='dynamical. I will try and compare my external JKTEBOP model with the PHOEBE dynamical rv.

[jsinkbaek]

For the PHOEBE dynamical model, there is zero discrepancy with the external JKTEBOP model

[kecnry]

Your changes that you described to the code seem like a reasonable fix to me. If I understand your last comments correctly (and sorry - I haven't had a chance to actually read or run the script yet), it sounds like there is still some remaining bug in the translation error since you're seeing identical results in your own script. PHOEBE writes to a temporary file that is then passed to jktebop - you can try to open that file (its easier if you change the code to not use a temporary file) and compare to your own file you created and see what parameters are causing the difference.

If you're willing, you can open a pull request with those changes (and then get "credit" for the code changes) and that might also help us track down the remaining discrepancy. To do that, you'll need to fork the repository to your own account, create a branch off of the master branch, push it to your fork, and then go to the "Pull Request" tab here and open a new PR. That'll then run all the tests and eventually make its way into the next bugfix release. Or if you'd rather, we should be able to quickly reproduce the changes you made ourselves and get a bugfix release out in the next week or so.

[aprsa]

We also need to have a test unit added for the PR, to test against such discrepancies; I can help with that if you'd like.

[jsinkbaek]

So, I am working with two notebooks here. One, which I attached, is a synthetic example to illustrate the first point. The other is a jumble of code where I have tested different things for an actual stellar system, which includes measured RVs. For all the comparisons with "external JKTEBOP", that is my own JKTEBOP fit to that stellar system.

@kecnry what you say is correct. My fit using JKTEBOP outside of PHOEBE matches the "traditional" PHOEBE dynamical RV model much much better than the PHOEBE jktebop backend.

[jsinkbaek]

I might have found the issue with the smaller discrepancy, I will just test something and be right back.

I suspect it is because currently K = something something / period where it should be K = something something / (period * sqrt(1-e^2))

[jsinkbaek]

Yep, that solves it, see attached. This is the two synthetic data-sets, with q=1 and q=0.5, compared with the dynamical PHOEBE model. The difference is now on same order as the differences between PHOEBE and JKTEBOP (external).

[aprsa]

Perfect, thanks!! Are you comfortable with trying to issue a PR yourself, with our help of course, or would you prefer that we take care of it on our side?

[jsinkbaek]

I will try ;)

[aprsa]

Awesome! :) The first step is to write a pytest unit for this -- my suggestion is to do precisely what you did above, i.e. test RV residuals and make sure they're all "close" -- i.e., within a specified tolerance. That test should live in the phoebe/nosetests directory, and a pytest run needs to pass the test.

[jsinkbaek]

So a test of jktebop vs phoebe backend? I could write something that selects semi-random orbital parameters, since I believe this probably was not observed as e=0 and q=1 behaves nicely.

[kecnry]

We might need to have the output from jktebop stored somewhere in the tests (since installing jktebop on the CI machinery might not be simple). But if you want to get the code changes up in a PR, we can also take over and write the tests.

[aprsa]

Right: storing the jktebop output in a standalone file, and testing phoebe's backends against it, sounds like the best solution to me. If you can broaden the test to other parameters (rather than just q/e), that would be even better as it would provide us with more coverage.

[jsinkbaek]

Ehh, is it not the JKTEBOP backend this test would be for? Having a fixed model would then defeat the purpose if code changes are introduced without recalculating it. Unless you wanted it to verify the PHOEBE dynamical RV against it instead?

It is getting a bit technical for me. Would it be ok if I just submitted the pull request for the bugfix and asked you to handle the tests?

[aprsa]

Sure thing!

The problem with the test is that jktebop isn't a python package but, rather, a standalone executable. That makes it tricky (but not impossible) to install on github's continuous integration (CI). As of right now, we do not have it installed, so if we wanted to test against jktebop in real time, we'd need to do that as well. I'm not opposed to trying to do that as well, but probably as a separate ticket. In the mean time, I figured that testing against, say, 3 or 4 different circumstances of precomputed jktebop models and jktebop backend in phoebe (RVs, but perhaps also LCs) would be a good enough start.

[kecnry]

That's a good point - to test the wrapper itself we will either need to get jktebop installed on the CI system or test the temporary jktebop input files that phoebe produces. But we can discuss and take over and help write those tests in the PR itself if you can just get your code changes submitted there. Thanks again!

[jsinkbaek]

I see that finite time of integration support for JKTEBOP (numint) is also on a To-Do list. I would like to use this feature. If I can implement it and submit a PR as a separate branch, would that be ok?

[aprsa]

Sounds perfect!

[jsinkbaek]

For fti:

I don't see any reference to dataset in the JKTEBOP backend code. Is that handled neatly such that the bundle there only sees one dataset (the one it needs)? Or will I run into trouble there if it is defined like a regular parameter, in case there is two light curves defined?

[kecnry]

Awesome! Let's move the fti discussion over to #700 just so the conversations aren't all mangled together (or feel free to also message us on the old workshop slack if you still have access to that).

[jsinkbaek]

I just issued a pull request (#701) for this functionality